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

posted by Helmut Homepage – Vienna, Austria, 2006-12-28 22:33 (6590 d 02:23 ago) – Posting: # 422
Views: 56,198

Dear all,

following this post and just out of curiosity I implemented SAS-code by Jones and Kenward (2000)[1] 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)[2]
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 )
  }

[image]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.

[image]

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 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes

Complete thread:

UA Flag
Activity
 Admin contact
23,359 posts in 4,905 threads, 1,671 registered users;
102 visitors (0 registered, 102 guests [including 22 identified bots]).
Forum time: 00:57 CET (Europe/Vienna)

Trust, but verify.    Anonymous (Russian proverb)

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