Power Analysis: Non-inferiority Cluster Randomised Trial

Sample Size and Power Calculation for Non-Inferiority Test

Null Hypothesis: The efficacy of the intervention arm (\(\pi_2\)) is inferior (worse) than the SOC arm (\(\pi_1\)) by margin of \(\delta\). \[H_0: \pi_{1} - \pi_{2} \ge \delta\]

Alternative Hypothesis: The efficacy of the intervention arm (\(\pi_2\)) is no worse than the SOC arm (\(\pi_1\)) by margin of \(\delta\). \[H_1: \pi_{1} - \pi_{2} < \delta\]

Standard Error for the Difference of Two Proportions:

\[SE = \sqrt{\hat{p_1}\cdot(1-\hat{p_1})/n_1 + \hat{p_2}\cdot(1-\hat{p_2})/n_2}\]

Standard Error for the Difference of Two Means: \[SE = \sigma\sqrt{1/n_1 + 1/n_2}\]

Test Statistics:

\[z = \frac{\hat{p_1}- \hat{p_2} - \delta}{SE}\]

Single-sided \((100-\alpha)\%\) Confidence Interval:

\[[-1,\hat{p_1} - \hat{p_2} + z_{1-\alpha} \cdot SE]\]

Sample Size for Individual-Randomization:

\[n_I=\frac{(z_{\alpha/2}+z_\beta)^2(\pi_0(1-\pi_0)+\pi_1(1-\pi_1))}{(\pi_0-\pi_1)^2}\]

Design Effect:

\[Deff=1+(m-1)\rho\]

Reference:

  • Richard J. Hayes, Lawrence H. Moulton (2017) Cluster Randomised Trials, 2nd Edition. Chapman & Hall/CRC Biostatistics Series
  • Blackwelder WC (1992) “Proving the null hypothesis” in clinical trials. Control Clin Trials 3:345-353.
  • Jennifer Schumi and Janet T Wittes (2011) Through the looking glass: understanding non-inferiority. Trials 2011, 12:106

Non-inferiority Cluster Randomised Trial

Richard J. Hayes & Lawrence H. Moulton. Cluster Randomised Trials, Second Edition, 2017, P. 156

  • This formula uses alpha value of 0.05 not 0.025, to ensure a standard 95% confidence interval on delta would excluded the -0.05 margin with 80% probability.
  • C: number of clusters
  • m: number of subjects per cluster
  • delta = pi_1 - pi_0: true difference between two arms
  • pi_1: proportion on the treatment arm
  • pi_0: proportion on the control arm
  • beta:
  • alpha:
  • k: coefficient of variance

[1] 0.7786671

Power with ICC

Assuming the proportion of outcome in the treatment arm is 0.3, then by enrolling 292 people into our study (approximately 124 in the treatment arm and 168 in the control arm), we anticipate having 80% power to detect a 60% increase in outcome in the treatment arm (i.e., relative risk [RR]=1.6, proportion outcome in treatment arm=0.48). These and all other power calculations assume a type I error rate of 0.05. These calculations treat people who are lost to follow-up (LTFU) as non-adherent, which is how they will be handled in the analyses. The sensitivity of these calculations to varying intraclass correlations, effects sizes, and levels of outcome in the treatment arm is shown in the table. Again, it should be noted that these are approximate calculations that will be made more precise during the pilot phase of the funding period.

With n=292 where we compute the power, we vary the effect size from 1.5, 1.6, 1.7; vary the intraclass correlation in the treatment group from 0.2, 0.25, 0.3; and vary the levels of outcome in the control arm from 0.2, 0.3, 0.4.

  • n_total: Total number of participants
  • n_0 & n_1: Number of participants on individual PrEP arm and group PrEP arm
  • rho_0 & rho_1: ICC assumption on individual PrEP arm and group PrEP arm
  • deff_0 & deff_1: Design effect
  • n_0_s & n_1_s: Sample size under individual randomization
  • pi_0: Proportion assumption for the individual PrEP arm
  • RR: Relative risk
  • Power
RowNumber n_total RR pi_0 n_0 rho_0 deff_0 n_0_s n_1 rho_1 deff_1 n_1_s Power
1 292 1.5 0.2 124 0.05 1.1 112 168 0.20 1.4 120 0.43
2 292 1.5 0.2 124 0.05 1.1 112 168 0.25 1.5 112 0.41
3 292 1.5 0.2 124 0.05 1.1 112 168 0.30 1.6 105 0.40
4 292 1.5 0.3 124 0.05 1.1 112 168 0.20 1.4 120 0.67
5 292 1.5 0.3 124 0.05 1.1 112 168 0.25 1.5 112 0.65
6 292 1.5 0.3 124 0.05 1.1 112 168 0.30 1.6 105 0.64
7 292 1.5 0.4 124 0.05 1.1 112 168 0.20 1.4 120 0.87
8 292 1.5 0.4 124 0.05 1.1 112 168 0.25 1.5 112 0.86
9 292 1.5 0.4 124 0.05 1.1 112 168 0.30 1.6 105 0.85
10 292 1.6 0.2 124 0.05 1.1 112 168 0.20 1.4 120 0.56
11 292 1.6 0.2 124 0.05 1.1 112 168 0.25 1.5 112 0.54
12 292 1.6 0.2 124 0.05 1.1 112 168 0.30 1.6 105 0.53
13 292 1.6 0.3 124 0.05 1.1 112 168 0.20 1.4 120 0.82
14 292 1.6 0.3 124 0.05 1.1 112 168 0.25 1.5 112 0.80
15 292 1.6 0.3 124 0.05 1.1 112 168 0.30 1.6 105 0.79
16 292 1.6 0.4 124 0.05 1.1 112 168 0.20 1.4 120 0.96
17 292 1.6 0.4 124 0.05 1.1 112 168 0.25 1.5 112 0.96
18 292 1.6 0.4 124 0.05 1.1 112 168 0.30 1.6 105 0.95
19 292 1.7 0.2 124 0.05 1.1 112 168 0.20 1.4 120 0.68
20 292 1.7 0.2 124 0.05 1.1 112 168 0.25 1.5 112 0.67
21 292 1.7 0.2 124 0.05 1.1 112 168 0.30 1.6 105 0.65
22 292 1.7 0.3 124 0.05 1.1 112 168 0.20 1.4 120 0.92
23 292 1.7 0.3 124 0.05 1.1 112 168 0.25 1.5 112 0.91
24 292 1.7 0.3 124 0.05 1.1 112 168 0.30 1.6 105 0.90
25 292 1.7 0.4 124 0.05 1.1 112 168 0.20 1.4 120 0.99
26 292 1.7 0.4 124 0.05 1.1 112 168 0.25 1.5 112 0.99
27 292 1.7 0.4 124 0.05 1.1 112 168 0.30 1.6 105 0.99

Matched Pair Design

Individual-Randomization \[n_I=\frac{(z_{\alpha/2}+z_\beta)^2(\pi_0(1-\pi_0)+\pi_1(1-\pi_1))}{(\pi_0-\pi_1)^2}\]

Design Effect \[Deff=1+(m-1)\rho\]

Relationship between ICC and Coefficient of Variance (k) \[\rho = \frac{k^2\pi}{1-\pi}\]

Matched-Pair Cluster-Randomization \[c_m=2 + (z_{\alpha/2} + z_{\beta})^2\frac{\pi_0(1-\pi_0)/m + \pi_1(1-\pi_1)/m + k_m^2(\pi_0^2 + \pi_1^2)}{(\pi_0-\pi_1)^2}\]

Reference: Richard J. Hayes, Lawrence H. Moulton. Cluster Randomised Trials, 2nd Edition. Chapman & Hall/CRC Biostatistics Series

Sample Size and Power
# Matched Pairs Patients per Clinic Rate on SOC RR Rate on TRT ICC Power
1 14 27 0.65 1.20 0.78 0.03 0.719
2 14 27 0.65 1.20 0.78 0.05 0.588
3 14 27 0.65 1.22 0.79 0.03 0.797
4 14 27 0.65 1.22 0.79 0.05 0.668
5 14 27 0.65 1.22 0.79 0.07 0.567
6 14 27 0.65 1.24 0.81 0.03 0.861
7 14 27 0.65 1.24 0.81 0.05 0.740
8 14 27 0.65 1.24 0.81 0.07 0.638
9 14 27 0.65 1.26 0.82 0.03 0.909
10 14 27 0.65 1.26 0.82 0.05 0.803
11 14 27 0.65 1.26 0.82 0.07 0.703
12 14 27 0.65 1.28 0.83 0.03 0.944
13 14 27 0.65 1.28 0.83 0.05 0.855
14 14 27 0.65 1.28 0.83 0.07 0.762
15 14 27 0.65 1.30 0.85 0.03 0.967
16 14 27 0.65 1.30 0.85 0.05 0.897
17 14 27 0.65 1.30 0.85 0.07 0.814
18 14 27 0.65 1.32 0.86 0.03 0.982
19 14 27 0.65 1.32 0.86 0.05 0.929
20 14 27 0.65 1.32 0.86 0.07 0.857
21 14 27 0.65 1.34 0.87 0.03 0.990
22 14 27 0.65 1.34 0.87 0.05 0.953
23 14 27 0.65 1.34 0.87 0.07 0.893
24 14 27 0.65 1.36 0.88 0.03 0.995
25 14 27 0.65 1.36 0.88 0.05 0.970
26 14 27 0.65 1.36 0.88 0.07 0.922
27 14 27 0.65 1.38 0.90 0.03 0.998
28 14 27 0.65 1.38 0.90 0.05 0.981
29 14 27 0.65 1.38 0.90 0.07 0.944
30 14 27 0.65 1.40 0.91 0.03 0.999
31 14 27 0.65 1.40 0.91 0.05 0.989
32 14 27 0.65 1.40 0.91 0.07 0.961 
p1 <- 0.65
rr <- seq(1.1,1.4,0.01)
c <- c(14)
n <- 27
icc <- c(0.03,0.05,0.07)
alpha <- c(0.05)
df <- expand.grid(n,p1,rr,c,icc,alpha)
df <- as.data.table(as.data.frame(df))
colnames(df) <- c("n","p1","rr","c","icc","alpha")

# k is calculated on p1
df <- df[
  ,p2 := p1*rr
][p2 < 0.95
][,k := sqrt(icc*(1-p1)/p1)
][,z_alpha := qnorm(1-alpha/2)
][,z_beta := sqrt((c - 2)*(p1-p2)^2/(p1*(1-p1)/n + p2*(1-p2)/n + k^2*(p1^2 + p2^2))) - z_alpha
# ][z_beta > 0
][,power:= pnorm(z_beta)
]
# dim(df)
# [,ni := (qnorm(1-0.05/2) + qnorm(0.8))^2*(p1*(1-p1) + p2*(1-p2))/(p1-p2)^2
# ][,deff := n*c/ni
# ][,deff2 := 2 + icc*(n-1)
# ][power >= 0.8
# ][,diff := power - 0.8
# ][order(c,n,p1,rr,icc,abs(diff))
# ][,rnb := 1:.N,by=list(c,n,p1,rr,icc)
# ][rnb == 1
# ][order(c,p1,rr,icc,n)
# ][,rnb2 :=1:.N,by=list(c,p1,rr,icc)
# ][rnb2 == 1
# ][order(p1,rr,icc,n,c)
# ][,rn:=1:.N]

# View(df) dim(df)
df$alpha <- factor(df$alpha)
df$icc <- factor(df$icc)
df_p80 <- as.data.frame(cbind(x1=0,y1=0.8,x2=1.26,y2=0.8))
# colnames(df_p80)

Color.Line <- "#666666"

p <- ggplot(data=df, aes(x=rr, y=power, group = icc, colour = icc))
p <- p + geom_line(size=2)
p <- p + geom_segment(aes(x = 1.26
                          , y = 0
                          , xend = 1.26
                          , yend = 0.90),linetype=2,colour= Color.Line)
p <- p + geom_segment(aes(x = 1.1
                          , y = 0.7
                          , xend = 1.26
                          , yend = 0.7),linetype=2,colour= Color.Line)
p <- p + geom_segment(aes(x = 1.1
                          , y = 0.8
                          , xend = 1.26
                          , yend = 0.8),linetype=2,colour= Color.Line)
p <- p + geom_segment(aes(x = 1.1
                          , y = 0.9
                          , xend = 1.26
                          , yend = 0.9),linetype=2,colour= Color.Line)
p <- p + labs(x = "Relative Risk"
              ,y="Power"
              ,title="Power of Study Design"
              )
p <- p + scale_colour_manual(values=c("#15274F","#115F11","#764115"))
p <- p + scale_x_continuous(breaks=c(1.1,1.2,1.26,1.3,1.4))
p <- p + scale_y_continuous(breaks=c(0.2,0.4,0.7,0.8,0.9)
                            ,labels = c("20%","40%","70%","80%","90%"))
p <- p + theme(legend.position="top")
p <- p + theme(panel.background = element_rect(fill = '#FFFFFF'
                                               , colour = '#EEEEEE'))
p

Power on Harzard

With design effect 1.1 on the control arm and design effect 1.4, 1.5, and 1.6 on the treatment arm, we calucalte power for the allocation of 124 on the control arm and 168 on the treatment arm, assuming 20% - 50% reduction on hazard (HR 0.8 to 0.5) for the treatment arm, and the proportion of LTFU for the control arm 20%, 30%, or 40%.

  • HR: Hazard Ratio
  • n_0, n_1, & n_total: Number of participants on the individual arm, group arm, and total
  • L_0, L_1, & L: Number of Events from individual PrEP arm, group arm, and total
  • pi_0: Proportion of LTFU on the individual PrEP arm
  • rho_0 and rho_1: ICC
  • deff_0 and deff_1: Design effect
  • Power
RowNumber HR pi_0 n_0 n_0_s L_0 n_1 rho_1 deff_1 n_1_s L_1 L Power
1 0.5 0.2 124 112 22 168 0.20 1.4 120 12 34 0.52
2 0.5 0.2 124 112 22 168 0.25 1.5 112 11 33 0.51
3 0.5 0.2 124 112 22 168 0.30 1.6 105 10 32 0.50
4 0.5 0.3 124 112 33 168 0.20 1.4 120 18 51 0.70
5 0.5 0.3 124 112 33 168 0.25 1.5 112 16 49 0.68
6 0.5 0.3 124 112 33 168 0.30 1.6 105 15 48 0.67
7 0.5 0.4 124 112 44 168 0.20 1.4 120 24 68 0.81
8 0.5 0.4 124 112 44 168 0.25 1.5 112 22 66 0.80
9 0.5 0.4 124 112 44 168 0.30 1.6 105 21 65 0.80
10 0.6 0.2 124 112 22 168 0.20 1.4 120 14 36 0.33
11 0.6 0.2 124 112 22 168 0.25 1.5 112 13 35 0.33
12 0.6 0.2 124 112 22 168 0.30 1.6 105 12 34 0.32
13 0.6 0.3 124 112 33 168 0.20 1.4 120 21 54 0.47
14 0.6 0.3 124 112 33 168 0.25 1.5 112 20 53 0.46
15 0.6 0.3 124 112 33 168 0.30 1.6 105 18 51 0.45
16 0.6 0.4 124 112 44 168 0.20 1.4 120 28 72 0.58
17 0.6 0.4 124 112 44 168 0.25 1.5 112 26 70 0.57
18 0.6 0.4 124 112 44 168 0.30 1.6 105 25 69 0.56
19 0.7 0.2 124 112 22 168 0.20 1.4 120 16 38 0.19
20 0.7 0.2 124 112 22 168 0.25 1.5 112 15 37 0.19
21 0.7 0.2 124 112 22 168 0.30 1.6 105 14 36 0.19
22 0.7 0.3 124 112 33 168 0.20 1.4 120 25 58 0.27
23 0.7 0.3 124 112 33 168 0.25 1.5 112 23 56 0.27
24 0.7 0.3 124 112 33 168 0.30 1.6 105 22 55 0.26
25 0.7 0.4 124 112 44 168 0.20 1.4 120 33 77 0.35
26 0.7 0.4 124 112 44 168 0.25 1.5 112 31 75 0.34
27 0.7 0.4 124 112 44 168 0.30 1.6 105 29 73 0.33
28 0.8 0.2 124 112 22 168 0.20 1.4 120 19 41 0.11
29 0.8 0.2 124 112 22 168 0.25 1.5 112 17 39 0.10
30 0.8 0.2 124 112 22 168 0.30 1.6 105 16 38 0.10
31 0.8 0.3 124 112 33 168 0.20 1.4 120 28 61 0.14
32 0.8 0.3 124 112 33 168 0.25 1.5 112 26 59 0.14
33 0.8 0.3 124 112 33 168 0.30 1.6 105 25 58 0.13
34 0.8 0.4 124 112 44 168 0.20 1.4 120 38 82 0.17
35 0.8 0.4 124 112 44 168 0.25 1.5 112 35 79 0.17
36 0.8 0.4 124 112 44 168 0.30 1.6 105 33 77 0.16

Count Negative Binomial

Because of the fact that the variance is much larger than the mean, we believe there exists overdispertion in the distribution of the count data – the number of times for events. Therefore, we assume it follows a negative binomial distribution. By fitting a negative binomial model on the previous survey results, we have attained an estimate for the dispersion parameter of .

\[ P(X=x) = \binom{x+r-1}{x}p^r(1-p)^x\]

\[\mu = \frac{r(1-p)}{p}\]

\[ p = \frac{r}{u+r}\]

\[\sigma^2 = \frac{r(1-p)}{p^2}\] \[\sigma^2 = \mu + \frac{\mu^2}{r}\] \[ r = \frac{\mu^2}{\sigma^2-\mu}\]

\[ Deff = 1 + (m-1)\rho\]

Given the dispersion parameter of , in case we prolong the obersavation period to one month, the negative binomial distribution will give us a monthly dispersion parameter of which is 4.36 times of the weekly dispersion parameter.

All participants will be followed up for 15 months, and their monthly exposures will be recorded. The way negative binomial works is that there exists significant large difference across individual healers. we assume that even though most of variance on exposure is due to individual practice pattern, some of variance should be due to seasonal and other monthly factors. Therefore conservatively we assume 60% to 80% (70% in average) of variance is due to individual practice pattern, with 15 repeated measures on each participant, we will expect a design effect of 10.8.

With 100 participants enrolled on each arm and 10% lost of follow up, 15 repeated measures for each participant with design effect of 10.8, we calculate power under 1% and 5% type I error based on negative binomial regression with 1000 simulation for each level of treatment effect size, from 40% to 80% rate ratio on blood exposure. The simulation is conducted in R.

library(MASS)
ng_power <- function(eff){
  set.seed(654321)
  for(rep in 1:1000) {
    # print(rep)
    n <- 90*1.388889
    r_week <- 0.075
    mu_week <- 0.88
    p <- r_week/(r_week+mu_week)
    mu_month <- mu_week*4.36
    r_month <- mu_month*p/(1-p)
    size <- r_month
    
    # rnbinom(n=103,size=0.06852744,prob=0.07206736)
    S_c <- rnbinom(
      n = n
      , size = size
      # , prob = p
      , mu = mu_month
      )
    S_t <- rnbinom(
      n = n
      , size = size
      # , prob = p*eff
      , mu = mu_month*eff
      )
    
    d1 <- rbind(
      data.frame(count=S_c,treatment=0)
      ,data.frame(count=S_t,treatment=1))
    M1 <- glm.nb(
      count ~ treatment
      ,data = d1
    )
    SM1 <- summary(M1)
    SM1 <- SM1[["coefficients"]]
    p1 <- SM1[2,4]
    p2 <- p1 < 0.01
    p1 <- p1 < 0.05
    if(rep==1){pv <- p1}else{pv <- c(pv,p1)}
    
    if(rep==1){pv2 <- p2}else{pv2 <- c(pv2,p2)}
    # return(pv)
  }
  return(c(sum(pv)/length(pv),sum(pv2)/length(pv)))
}

seq_eff <- seq(0.4,0.8,by=0.025)
# x <- ng_power(0.5)
## length(seq_eff)

power_exp <- data.frame(
  reduction = seq_eff
  ,power = rbind(
    ng_power(seq_eff[1])
    ,ng_power(seq_eff[2])
    ,ng_power(seq_eff[3])
    ,ng_power(seq_eff[4])
    ,ng_power(seq_eff[5])
    ,ng_power(seq_eff[6])
    ,ng_power(seq_eff[7])
    ,ng_power(seq_eff[8])
    ,ng_power(seq_eff[9])
    ,ng_power(seq_eff[10])
    ,ng_power(seq_eff[11])
    ,ng_power(seq_eff[12])
    ,ng_power(seq_eff[13])
    ,ng_power(seq_eff[14])
    ,ng_power(seq_eff[15])
    ,ng_power(seq_eff[16])
    ,ng_power(seq_eff[17])
  )
)


power_exp$reduction <- WuPercent(power_exp$reduction,1)
power_exp$power.1 <- WuPercent(power_exp$power.1,1)
power_exp$power.2 <- WuPercent(power_exp$power.2,1)

knitr::kable(
  power_exp
  ,caption = "Power Calculation for Blood Exposure, 90 per arm"
  ,col.names = c("Rate Ratio","Power with 5% Type I Error", "Power with 1% Type I Error")
  ,align = c("r","r","r")
)

R sessionInfo

R version 4.2.0 (2022-04-22) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Ubuntu 20.04.3 LTS

Matrix products: default BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0 LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale: [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
[4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
[7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C

attached base packages: [1] stats graphics grDevices utils datasets methods base

other attached packages: [1] powerSurvEpi_0.1.3 Wu_0.0.0.9000 flexdashboard_0.5.2 [4] lme4_1.1-29 Matrix_1.4-0 mgcv_1.8-38
[7] nlme_3.1-152 png_0.1-7 scales_1.2.0
[10] nnet_7.3-16 labelled_2.9.1 kableExtra_1.3.4
[13] plotly_4.10.0 gridExtra_2.3 ggplot2_3.3.6
[16] DT_0.23 tableone_0.13.2 magrittr_2.0.3
[19] lubridate_1.8.0 dplyr_1.0.9 plyr_1.8.7
[22] data.table_1.14.2 rmdformats_1.0.4 knitr_1.39

loaded via a namespace (and not attached): [1] httr_1.4.3 sass_0.4.1 tidyr_1.2.0 jsonlite_1.8.0
[5] viridisLite_0.4.0 splines_4.2.0 bslib_0.3.1 assertthat_0.2.1 [9] highr_0.9 yaml_2.3.5 pillar_1.8.0 lattice_0.20-45
[13] glue_1.6.2 digest_0.6.29 rvest_1.0.2 minqa_1.2.4
[17] colorspace_2.0-3 htmltools_0.5.3 survey_4.1-1 pkgconfig_2.0.3
[21] haven_2.5.0 bookdown_0.27 purrr_0.3.4 webshot_0.5.3
[25] svglite_2.1.0 pracma_2.3.8 tibble_3.1.8 farver_2.1.1
[29] generics_0.1.3 ellipsis_0.3.2 withr_2.5.0 klippy_0.0.0.9500 [33] lazyeval_0.2.2 cli_3.3.0 survival_3.2-13 evaluate_0.15
[37] fansi_1.0.3 MASS_7.3-54 forcats_0.5.1 xml2_1.3.3
[41] tools_4.2.0 hms_1.1.1 mitools_2.4 lifecycle_1.0.1
[45] stringr_1.4.0 munsell_0.5.0 compiler_4.2.0 jquerylib_0.1.4
[49] systemfonts_1.0.4 rlang_1.0.4 grid_4.2.0 nloptr_2.0.3
[53] rstudioapi_0.13 htmlwidgets_1.5.4 rmarkdown_2.14 boot_1.3-28
[57] gtable_0.3.0 DBI_1.1.2 R6_2.5.1 fastmap_1.1.0
[61] utf8_1.2.2 stringi_1.7.8 Rcpp_1.0.9 vctrs_0.4.1
[65] tidyselect_1.1.2 xfun_0.31