Bioequivalence and Bioavailability Forum • FDA group model in R

FDA group model in R [Two-Stage / GS Designs]

posted by mittyri – Russia, 2016-10-11 15:43 (3103 d 02:08 ago) – Posting: # 16723
Views: 39,516
(edited on 2016-10-11 16:53)

Dear zizou,

Thank you for suggestion!
The FDA group model for PHX is described by Helmut here:

fixed:  Group+Sequence+Sequence(Group)+Period(Group)+Treatment+Treatment×Group

random: Subject(Sequence×Group)

I think in this case PHX easily switches from GLM to MIXED procedure and does not produce SS tables.
We all know that the GLM procedure should be used for conventional FDA model. Could we use the same procedure (GLM) here?

The relationship between GLM/MIXED is described very well in 'must-read' post by ElMaestro

I executed the code with unequal sequences (removing Subj18 from the reference dataset):
PHX Sequential tests:

  Dependent Hypothesis Numer_DF Denom_DF        F_stat     P_value

1   Ln(Var)        int        1       13 2830.51538137 0.000000000

2   Ln(Var)      Group        1       13    1.01359006 0.332416756

3   Ln(Var)        Seq        1       13    0.69165428 0.420616500

4   Ln(Var)  Group*Seq        1       13    0.12190952 0.732563224

5   Ln(Var)  Group*Per        2       13    2.99449184 0.085186463

6   Ln(Var)        Trt        1       13    3.42927592 0.086890932

7   Ln(Var)  Trt*Group        1       13    0.34285051 0.568214408

PHX Partial tests:

  Dependent Hypothesis Numer_DF Denom_DF        F_stat     P_value

1   Ln(Var)        int        1       13 2803.64675463 0.000000000

2   Ln(Var)      Group        1       13    1.14134985 0.304804606

3   Ln(Var)        Seq        1       13    0.65955233 0.431338627

4   Ln(Var)  Group*Seq        1       13    0.12190952 0.732563224

5   Ln(Var)  Group*Per        2       13    3.18289277 0.074972760

6   Ln(Var)        Trt        1       13    3.53470601 0.082691107

7   Ln(Var)  Trt*Group        1       13    0.34285051 0.568214408

muddle <- lm(log(Var)~Group+Seq+Seq*Group+Subj%in%(Seq)+Per%in%Group+Trt+Group*Trt, data=data)

anova(muddle)

              Df  Sum Sq   Mean Sq  F value     Pr(>F)    

Group           1 0.30128 0.3012782 40.66580 2.4299e-05 ***

Seq             1 0.20559 0.2055864 27.74955 0.00015211 ***

Trt             1 0.02101 0.0210135  2.83635 0.11599640    

Group:Seq       1 0.03624 0.0362362  4.89108 0.04551863 *  

Group:Per       2 0.04876 0.0243815  3.29096 0.06975471 .  

Group:Trt       1 0.00254 0.0025401  0.34285 0.56821441    

Group:Seq:Subj 13 3.86410 0.2972387 40.12055 2.5513e-08 ***

Residuals      13 0.09631 0.0074086

drop1(muddle,test="F")


Single term deletions

Model:

log(Var) ~ Group + Seq + Seq * Group + Subj %in% (Seq * Group) + 

    Per %in% Group + Trt + Group * Trt

               Df Sum of Sq     RSS       AIC  F value     Pr(>F)    

<none>                      0.09631 -157.4617                        

Group:Per       2   0.04716 0.14347 -147.9107  3.18289   0.074973 .  

Group:Trt       1   0.00254 0.09885 -158.5766  0.34285   0.568214    

Group:Seq:Subj 13   3.86410 3.96042  -57.1004 40.12055 2.5513e-08 ***

and CI's are far from equality... :-(

PS: I tried to switch to the mixed model

library(nlme)

Group222FDA<- function(data, alpha=0.05)

{

  data$Per  <- factor(data$Per)  # Period

  data$Seq  <- factor(data$Seq)  # Sequence

  data$Trt  <- factor(data$Trt)  # Treatment

  data$Subj <- factor(data$Subj) # Subject

  data$Group  <- factor(data$Group)

  

  muddlelme <- lme(log(Var) ~ Group + Seq + Seq %in% Group

                   + Per %in% Group + Trt +  Trt * Group , 

                   random=~1|Subj, data=data, method="REML")

  lnPE     <- coef(muddlelme)[1, "TrtT"]


  lnCI     <- confint(muddlelme,c("TrtT"), level=1-2*alpha)[[1]][1:2]


  PE <- exp(lnPE)

  CI <- exp(lnCI)

  typeI <- anova(muddlelme)

  # output

  options(digits=8)

  cat("Type I sum of squares: \n")

  print(typeI)

  cat("\nBack-transformed PE and ",100*(1-2*alpha))

  cat("% confidence interval\n")

  cat("Point estimate (GMR).(%).................:", 

      formatC(100*PE, format="f", digits=2),"\n")

  cat("Lower confidence limit.(%)...............:", 

      formatC(100*CI[1], format="f", digits=2) ,"\n")

  cat("Upper confidence limit.(%)...............:", 

      formatC(100*CI[2],format="f", digits=2) ,"\n")

}

data <- read.delim("https://static-content.springer.com/esm/art%3A10.1208%2Fs12248-014-9661-0/MediaObjects/12248_2014_9661_MOESM1_ESM.txt")



Group1 <- subset(data, Subj <= 9)

Group1$Group <- 1

Group2 <- subset(data, Subj > 9& Subj < 18)

Group2$Group <- 2

data<-rbind(Group1, Group2)

Group222FDA(data, alpha=0.05)

The ANOVA type I is the same as in PHX, but PE and CIs are not. I'm in stuck....

—
Kind regards,
Mittyri

Complete thread:

Potvin C in the EU Helmut 2013-04-16 17:40
- Potvin C in the EU ElMaestro 2013-04-17 02:31
  - Potvin C in the EU Helmut 2013-04-17 12:23
    - 2 Groups model FDA d_labes 2013-04-17 13:10
      - 2 Groups model FDA Helmut 2013-04-17 16:37
        
        Group effects obsolete? d_labes 2013-04-18 09:55
        
        Group effects FDA/EMA Helmut 2013-04-19 14:34
        
        Group effects FDA/EMA mittyri 2014-08-21 15:34
        
        Group effects FDA/EMA ElMaestro 2014-08-21 17:31
        
        Group effects EMA mittyri 2014-10-01 11:25
        
        Group effects EMA ElMaestro 2014-10-01 11:31
        Group effects EMA Helmut 2014-10-01 13:24
        
        Group effects EMA VStus 2016-10-06 22:06
        
        Group effects EMA ElMaestro 2016-10-07 12:45
        
        Group effects EMA mittyri 2016-10-07 15:06
        
        Fixed effects model with Group term mittyri 2016-10-09 13:27
        
        Fixed effects model: changing the F, p values mittyri 2016-10-10 18:33
        
        Fixed effects model: changing the F, p values ElMaestro 2016-10-10 20:10
        
        Fixed effects model: changing the F, p values Helmut 2016-10-11 00:17
        Holy War of type III d_labes 2016-10-11 13:17
        
        Significance of Group effect in Russia: why the type III is so 'important' mittyri 2016-10-11 15:07
        
        Fixed effects model: changing the F, p values zizou 2016-10-11 01:06
        
        FDA group model in Rmittyri 2016-10-11 13:43
        
        FDA group model in R ElMaestro 2016-10-12 01:07
        FDA group model in R VStus 2016-10-17 15:32
        
        FDA group model in R ElMaestro 2016-10-17 18:58
        
        FDA group model in R VStus 2016-10-20 13:57
        
        Fixed effects model with Group term VStus 2016-10-20 13:54
        
        Group effects FDA/EMA zizou 2016-12-31 02:54
        
        Group effects FDA/EMA Helmut 2017-06-03 14:57
        
        2 Groups model FDA hiren379 2013-08-22 14:32
        
        Homework Helmut 2013-08-23 11:53