R-code for finding simultaneous upper confidence bounds for relative risks in multiple comparisons to control:

You can automatically load all required R code by typing:

source(“http://sites.williams.edu/bklingen/files/2012/10/simRR_Oct2012.r”)

Alternatively, you can source the following file, which contains all necessary functions, into R: simRR_Oct2012

Attention: For the computation of Dunnett and Berger-Boos adjusted simultaneous confidence bounds, you need to download the “mvtnorm” package for multivariate normal calculations. (simRR.R will call the mvtnorm package automatically, so you just need to have this package installed)

Here is the example that’s analyzed in the manuscript:

> source(“http://lanfiles.williams.edu/~bklingen/simRR/simRR.R”)
> y <- c(70,26,24) # successes
> n <- c(3900,3900,3900) #sample sizes
> RRci.S2.Unadj(y,n,conflev=0.975) # Unadjusted score bounds
UB1 UB2
0.5790503 0.5416130
> RRci.S2.Sid(y,n,conflev=0.975) # Sidak adj. simult. score bounds
UB1 UB2
0.6160457 0.5772260
> RRci.S2.Dun(y,n,conflev=0.975) # Dunnett adj. simult. score bounds
UB1 UB2
0.6149397 0.5761532
> RRci.S2.BB(y,n,conflev=0.975) # Berger-Boos adj. simult. score bounds
UB1 UB2
0.6138859 0.5751453
> RRci.T1.Sid(y,n,conflev=0.975) # Sidak adj. simult. Wald bounds
UB1 UB2
0.6195966 0.5808495
> RRci.T1.Dun(y,n,conflev=0.975) # Dunnett adj. simult. Wald bounds
UB1 UB2
0.6184456 0.5797378
> RRci.T1.BB(y,n,conflev=0.975) # Berger-Boos adj. simult. Wald bounds
UB1 UB2
0.6173531 0.5786832

Evaluating coverage and power through simulation:

The R code above (simRR.R) contains a function (cov.pow.calc) that estimates the simultaneous coverage rate of the simultaneous upper bounds for various methods, for given value of p0 (the successs probability in the control group), given values of the true relative risks (RR) and given sample size (n). It also evaluates the power, i.e., the probability of the upper bounds being below a given treshold (RR.null).

For example, it is expected that the success probability in the control group is 10% (p0=0.1), that the relative risk with the first and second treatment is RR=c(0.7,0.6), respectively and that the sample sizes in the three groups (first group is control) is n=c(250,500,500). Additionally, we want to find the power for the upper bounds of the three relative risks being below RR.null=c(1,1), when the simultaneous coverage is controlled at conflev=0.975. For the Berger-Boos approach, the prec=0.001 determines how fine the grid for the numerical search should be (the smaller, the more precise). By default, unadjusted, Dunnet-adjusted and Berger-Boos adjusted simultaneous upper confidence bounds based on the score method are computed:

> result <- sim.cov.pow(p0=0.1, RR=c(0.7,0.6), n=c(250,500,500), conflev=0.975, RR.null=c(1,1), sims=1000, prec=0.001)
> result$mean.coverage
S2.Unadj S2.Dun S2.BB
0.953 0.981 0.980
> result$mean.pow.least
S2.Unadj S2.Dun S2.BB
0.555 0.454 0.459
> result$mean.pow.all
S2.Unadj S2.Dun S2.BB
0.239 0.169 0.169

This returns the estimated coverage probability for the unadjusted, Dunnett-adjusted and Berger-Boos adjusted simultaneous upper bounds with the score method (test statistics S2). E.g., the estimated coverage probabilities under the asumed values for the true parameters when inverting the score statistic using the Berger-Boos approach for the multiplicity adjustment is 98.0% (compare to a nominal coverage of 97.5%).

Additional statistics can be requested with the which.stat option. Available statitics are (the first three are the default): “S2.Unadj”, “S2.Sid”, “S2.Dun”, “S2.BB”, “T1.Sid”, “T1.Dun”, “T1.BB”. (E.g., which.stat=c(“S2.Unadj”, “S2.Sid”, “T1.Sid”) ).

In addition, by specifying RR.null=c(1,1) we estimated the power that at least one of the upper bounds is below its superiority margin (power.least) or that all upper bounds are below their respective superiority margin (1 in our case, power.all). For instance, under the assumed values for the true parameters, the power to see that the true relative risk is below 1 for at least one of the comparisons to control is 45.9% with the score confidence interval and the Berger-Boos multiplicity adjustment method.

 

Coming sometime in the future: R code for computing simultaneous lower bounds (please email me if there is a need)