## R-Code for Power in 2×2 Cross-over [R for BE/BA]

Dear all,

following this post and just out of curiosity I implemented SAS-code by Jones and Kenward (2000) in R (v2.4.1).
To quote the beginning of Jones'/Kenward's SAS-Code (which is copyrighted material and therefore not given here; but you can afford buying a book for sure):
/*** WARNING : PROGRAM OFFERED FOR USE WITHOUT ANY GUARANTEES    ***/ /*** NO LIABILITY IS ACCEPTED FOR ANY LOSS RESULTING FROM USE OF ***/ /*** THIS SET OF SAS INTRUCTIONS                                 ***/

Their example (pp 336-337) yields from a within-subject standard deviation of 0.355 in 58 subjects and a point estimate of 1.00 power=90.4%.
Following R-code gives 90.3988%.
  # significance level   a      <- 0.05   # variance of difference of two observations on log scale   # sigmaW = within-subjects standard deviation   sigmaW <- 0.355   s      <- sqrt(2)*sigmaW   # total number of subjects (needs to be a multiple of 2)   # for power calculation   n      <- 58   # error degrees of freedom for AB/BA crossover with   # n subjects in total   n2     <- n-2   # ratio = mu_T/mu_R   ratio  <- 1.00   # calculate power   t1     <- qt(1-a,n2)   t2     <- -t1   nc1    <- sqrt(n)*(log(ratio)-log(0.8))/s   nc2    <- sqrt(n)*(log(ratio)-log(1.25))/s   df     <- sqrt(n2)*(nc1-nc2)/(2*t1)   prob1  <- pt(t1,df,nc1)   prob2  <- pt(t2,df,nc2)   answer <- prob2-prob1   power  <- answer*100   power

The following code reproduces the power-curves for n=8/12/18/24/30/36/40/60 and CV=20% given in Fig.1c by Diletti et al. (1991)
 a       <- 0.05       # alpha CV      <- 0.20       # intra-subject coefficient of variation # calculate within-subjects standard deviation # from CV sigmaW  <- sqrt(log(CV^2+1)) s       <- sqrt(2)*sigmaW Theta1  <- 0.80       # lower acceptance limit Theta2  <- 1/Theta1   # upper acceptance limit Step    <- 0.005      # Interval ssx     <- 0.65       # x-coordinate of sample size labels in % CVLabel <- paste("CV =",as.character(CV*100),"%" ) # vector of sample sizes for (n in c(8,12,18,24,30,36,48,60))   {   n2    <- n-2   # vector of ratios from 'Theta1' to 'Theta2', interval 'Step'   assign("ratio",seq(Theta1,Theta2,Step))   t1    <- qt(1-a,n2)   t2    <- -t1   nc1   <- sqrt(n)*(log(ratio)-log(Theta1))/s   nc2   <- sqrt(n)*(log(ratio)-log(Theta2))/s   df    <- sqrt(n2)*(nc1-nc2)/(2*t1)   prob1 <- pt(t1,df,nc1)   prob2 <- pt(t2,df,nc2)   power <- prob2-prob1   if (n == 8)     # set-up of plot, first curve     {       plot(ratio, power, axes=F, type="l", lwd=2, col=2,         xlim=c(Theta1,Theta2), ylim=c(0,1),         main="2×2 Cross-over", xlab="µT/µR", ylab="Power" )       text(1, 0.4, CVLabel, cex=1.25, adj=0.5 )       axis(1, at=seq(Theta1,Theta2,Step*10),         labels=seq(Theta1,Theta2,Step*10), pos=0 )       axis(2, at=seq(0,1,0.1), labels=seq(0,1,0.1), las=1 )       abline(h=0.05, col = "black", lty = "dotted" )       abline(h=seq(0.6,1,0.1), col = "black", lty = "dotted" )       abline(v=seq(Theta1,Theta2,Step*10), col = "black", lty = "dotted" )       box()     }   else            # add curves to plot     lines(ratio, power, lwd=2, col=2 )     if (ratio[as.integer(length(ratio)*ssx)] )       text(ratio[as.integer(length(ratio)*ssx)],         power[as.integer(length(ratio)*ssx)], n, adj=0.5 )   } Values calculated by the R-code for n=24 match Diletti et al.'s Fig.1c:
+------+-------+ |  PE  | Power | +------+-------+ | 0.85 | 28.16 | | 0.90 | 64.08 | | 0.95 | 89.19 | | 1.00 | 96.28 | | 1.05 | 89.90 | | 1.10 | 70.02 | | 1.15 | 41.92 | | 1.20 | 18.30 | +------+-------+

Please note, this code is not validated, but only a toy you can play around with!

With some ambition you can modify the code (essentially varying the point estimate instead of the sample size) in order to match also Diletti et al.'s Fig.2c. Now it should be only a small step to re-calculate their table of sample sizes...
1. Jones B and MG Kenward
Design and Analysis of Cross-Over Trials
Chapman & Hall/CRC, Boca Raton, (2nd Edition 2000)
2. Diletti E, Hauschke D and VW Steinijans
Sample size determination for bioequivalence assessment by means of confidence intervals
Int J Clin Pharm Ther Toxicol 29/1, 1-8 (1991)

Dif-tor heh smusma 🖖
Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes

### Complete thread:

Admin contact
21,540 posts in 4,502 threads, 1,525 registered users;
online 27 (0 registered, 27 guests [including 6 identified bots]).
Forum time: Thursday 21:18 UTC (Europe/Vienna)

Laws are like sausages – it is better
not to see them being made.    Otto von Bismarck

The Bioequivalence and Bioavailability Forum is hosted by Ing. Helmut Schütz 