Bioequivalence and Bioavailability Forum • Oops.

Oops. [Regulatives / Guidelines]

posted by Helmut – Vienna, Austria, 2011-11-04 17:05 (5350 d 15:11 ago) – Posting: # 7607
Views: 11,125

Dear simulants!

❝ or what I would do

❝ FullModel <- lm(log(y1) ~ 0 + seq1 + sub1 %in% seq1 + per1 + trt1)

Oh THX! Affects also the T-R. Latest version (flexible sample sizes, corrected T/R1 and T/R2):

set.seed(12041958) # keep this line to compare results only

sims   <- 10000    # <- no of sims

size   <- 12       # <- no of subjects



# T/R1 = 95%, T/R2 = 105.26%

mean1 <- 0.95;   cv1 <- mean1/100*30 # <- means & (total!) CVs here

mean2 <- 1.00;   cv2 <- mean2/100*30 # <-

mean3 <- 0.95^2; cv3 <- mean3/100*30 # <-



seq1 <- as.factor(rep(rep(1:6, each = 3), size/6))

sub1 <- as.factor(rep(1:size, each = 3))

per1 <- as.factor(rep(1:3, size))

trt1 <- as.factor(rep(c(1,2,3, 2,3,1, 3,1,2, 1,3,2, 2,1,3, 3,2,1), size/6))



# massacrating out everything related to R2

# means that when T was before R1 we now have a seq 1 otherwise a seq 2

seq2 <- as.factor(rep(c(1,1, 2,2, 1,1, 1,1, 2,2, 2,2), size/6))

sub2 <- as.factor(rep(1:size, each = 2))

per2 <- as.factor(rep(1:2, size))

trt2 <- as.factor(rep(c(1,2, 2,1, 1,2, 1,2, 2,1, 2,1), size/6))



# see Martin's post: http://forum.bebac.at/mix_entry.php?id=5272

meanl1 <- log(mean1)-0.5*log(cv1^2+1); sdl1 <- sqrt(log(cv1^2+1))

meanl2 <- log(mean2)-0.5*log(cv2^2+1); sdl2 <- sqrt(log(cv2^2+1))

meanl3 <- log(mean3)-0.5*log(cv3^2+1); sdl3 <- sqrt(log(cv3^2+1))



RM.be  <- 0; FM.be <- 0

FMCVintra <- NULL; RMCVintra <- NULL

for (iter in 1:sims) # T is 1, R1 is 2 and R2 is 3

  {

  y1 = rlnorm(n=size*3, meanlog=meanl1, sdlog=sdl1)

  for (i in 1:size*3) if (trt1[i ]==2) y1[i ]=rlnorm(n=1,meanlog=meanl2,sdlog=sdl2)

  for (i in 1:size*3) if (trt1[i ]==3) y1[i ]=rlnorm(n=1,meanlog=meanl3,sdlog=sdl3)

  FullModel <- lm(log(y1) ~ 0 + seq1 + sub1 %in% seq1 + per1 + trt1)

  FMDelta   <- mean(log(y1)[trt1==1])-mean(log(y1)[trt1==2])

  FMdf      <- aov(FullModel)$df.residual

  FMMSE     <- summary(FullModel)$sigma^2/FMdf

  FMlo      <- 100*exp(FMDelta - qt(1-0.05,FMdf)*sqrt(2*FMMSE/size))

  FMhi      <- 100*exp(FMDelta + qt(1-0.05,FMdf)*sqrt(2*FMMSE/size))

  FMCI      <- FMhi - FMlo

  FMCVintra <- c(FMCVintra, sqrt(exp(FMMSE)-1))

  if ((FMlo>=80) && (FMhi<=125)) FM.be <- FM.be + 1

  # do the pruning!

  y2 <- y1

  length(y2) <- size*2

  j <- 0

  for (i in 1:size*3) if (trt1[i ] != 3)

    { j     <- j+1;

      y2[j] <- y1[i ]
    }

  ReducedModel <- lm(log(y2) ~ 0 + seq2 + sub2 %in% seq2 + per2 + trt2) # EMA?!

  RMDelta   <- mean(log(y2)[trt2==1])-mean(log(y2)[trt2==2])

  RMdf      <- aov(ReducedModel)$df.residual

  RMMSE     <- summary(ReducedModel)$sigma^2/RMdf

  RMlo      <- 100*exp(RMDelta - qt(1-0.05,RMdf)*sqrt(2*RMMSE/size))

  RMhi      <- 100*exp(RMDelta + qt(1-0.05,RMdf)*sqrt(2*RMMSE/size))

  RMCI      <- RMhi - RMlo

  RMCVintra <- c(RMCVintra, sqrt(exp(RMMSE)-1))

  if ((RMlo>=80) && (RMhi<=125)) RM.be <- RM.be + 1

}

if (RM.be>FM.be){Outcome <- "liberal."}else{Outcome <- "conservative."}

result <- paste(paste(

          " Reduced model (EMA: R2 dropped) compared to full model (Williams’ design)\n"),

          paste("Full model passed BE:    ",

          format(100*FM.be/sims,digits=4,nsmall=2,zero.print=TRUE),

          "%, CVintra: ",

          format(100*median(FMCVintra),digits=2,nsmall=2,

          zero.print=TRUE),"%\n"),

          paste("Reduced model passed BE: ",

          format(100*RM.be/sims,digits=4,nsmall=2,zero.print=TRUE),

          "%, CVintra: ",

          format(100*median(RMCVintra),digits=2,nsmall=2,

          zero.print=TRUE),"%\n"),

          paste("Assessment: Reduced model is",Outcome,"\n"),

          paste("Based on",sims,"simulated BE studies; n =",size,"each.","\n"))

cat(result)

Results:

       T  R1  R2   Full BE (CV_intra)  Red. BE (CV_intra)

Sc.1  30  30  30   87.40% (6.17%)    82.22% ( 8.69%) Red. = conservative

Sc.2  30  50  30   81.63% (6.71%)    70.28% (10.31%) Red. = conservative

Sc.3  30  30  50   86.61% (6.65%)    78.43% ( 9.19%) Red. = conservative

I have no idea how to simulate CV_intra; CV_total is not what I really want (to see if reference’s different variances have an influence on the result).

—
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:

Anticonservativism?! ElMaestro 2010-02-06 17:10 [Regulatives / Guidelines]
- Conservativism? Helmut 2011-11-03 21:12
  - In hindsight... ElMaestro 2011-11-03 21:43
    - In hindsight... Helmut 2011-11-04 01:56
      - In hindsight... ElMaestro 2011-11-04 08:41
      - rlnorm simulates what? d_labes 2011-11-04 09:40
        
        Oops.Helmut 2011-11-04 16:05
        
        Oops. ElMaestro 2011-11-04 16:22
        
        Simulation of intra-subject variability d_labes 2011-11-07 11:16
        
        Simulation of intra-subject variability ElMaestro 2011-11-08 11:02
        
        Oops. Oops. d_labes 2011-11-25 13:54
        
        Another Oops. Helmut 2011-11-25 14:23
- Simul Ants questions d_labes 2011-11-04 15:46
  - Simul Ants questions ElMaestro 2011-11-04 16:06
  - Simul Ants questions Helmut 2011-11-04 16:09
- Liberal Conservatives d_labes 2011-11-08 11:31
  - Liberal Conservatives martin 2011-11-08 20:50
    - Liberal Conservatives Helmut 2011-11-08 22:56
      - intra-subject correlation martin 2011-11-09 09:01
        
        intra-subject correlation Helmut 2011-11-25 17:16