Bioequivalence and Bioavailability Forum • In hindsight...

In hindsight... [Regulatives / Guidelines]

posted by Helmut – Vienna, Austria, 2011-11-04 02:56 (4930 d 19:04 ago) – Posting: # 7601
Views: 9,238

Salute!

❝ I am not sure why ElMaestro used a mixed model here; doesn't this specification just give the same as a linear model with seq1, per1, trt1, and sub1 as fixed factors :confused: ?

In terms of the residual variance – yes. Good ol’ days of subjects as random. Met last week another statistical pro who was upset upon the ‘all fixed’ story. ;-)

❝ I think we need to take degrees of freedom into consideration too. […] I think the df's differ between the two scenarios.

Right.

❝ I hope this does not lead into a discussion about how df's are calculated for a mixed model (R will not do Satterthwaite etc.),

No, no. Everything is balanced here. 20 dfs in the Williams’ design and 10 after R2 is thrown out.

Another try:

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

sims   <- 10000    # <- no of sims here



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

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

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

mean3 <- 1/0.95; cv3 <- mean3/100*30 # <-



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

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

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

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



# massacrating out everything related to R2

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

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

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

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

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



# 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))



Liberal <- 0; Conserv <- 0

FMCVintra <- NULL; RMCVintra <- NULL

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

  {

  y1 = rlnorm(n=36, meanlog=meanl1, sdlog=sdl1)

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

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

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

  FMDelta   <- mean(y1[trt1==1])-mean(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/12))

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

  FMCI      <- FMhi - FMlo

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

  # do the pruning!

  y2 <- y1

  length(y2) <- 24

  j <- 0

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

    { j     <- j+1;

      y2[j] <- y1[i]
    }

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

  RMDelta   <- mean(y2[trt2==1])-mean(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/12))

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

  RMCI      <- RMhi - RMlo

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

  if (RMCI  < FMCI){Liberal <- Liberal + 1}else{Conserv <- Conserv + 1}

}

result     <- paste(paste(" Eliminated R2 (EMA) compared to full model (Williams’ design)\n"),

              paste("Conservative: ", 100*Conserv/sims,

              "%, CVintra: ",round(100*median(FMCVintra), digits=2),"%\n"),

              paste("Liberal:      ", 100*Liberal/sims,

              "%, CVintra: ",round(100*median(RMCVintra), digits=2),"%\n"),

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

cat(result)

Sorry for the lengthy code; couldn’t do better.

       T  R1  R2   Cons. (CV_intra)  Liberal (CV_intra)

Sc.1  30  30  30   91.57% (6.56%)   8.43% ( 8.58%)

Sc.2  30  50  30   96.90% (7.83%)   3.10% (11.27%)

Sc.3  30  30  50   62.87% (8.26%)  37.13% ( 8.58%)

Hhm. Kannitverstan.

—
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