Bioequivalence and Bioavailability Forum • Fixed effects model: changing the F, p values

Fixed effects model: changing the F, p values [Two-Stage / GS Designs]

posted by mittyri – Russia, 2016-10-10 20:33 (3188 d 10:07 ago) – Posting: # 16710
Views: 40,779

Dear All,

I cannot change the previous post, so please find my corrected code
I've found my mistake regarding Fs
As ElMaestro suggested, all beween-subject factors should be assessed vs. MSE_subject:

# http://link.springer.com/article/10.1208/s12248-014-9661-0

Analyse222BE <- function (data, alpha=0.05, Group = FALSE) { 

  

  

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

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

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

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

  # be explicite

  ow <- options()

  options(contrasts=c("contr.treatment","contr.poly"))

  if(!Group) {

    muddle <- lm(log(Var)~Trt+Per+Seq+Subj, data=data) # conventional model

  }  else {

    data$Group  <- factor(data$Group)  # Group

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

  }

  

  # in the standard contrasts option "contr.treatment" 

  # the coefficient TrtT is the difference T-R since TrtR is set to zero

  lnPE     <- coef(muddle)["TrtT"]


  lnCI     <- confint(muddle,c("TrtT"), level=1-2*alpha)

  typeI    <- anova(muddle)

  

  #--- Type III sum of squares like SAS

  typeIII  <- drop1(muddle,test="F")

  typeIII$AIC <- NULL # null out this collumn

  # typeIII 3rd column is RSS - residual SS of model fit without term 

  # we make Mean Sq in that column

  names(typeIII)[3] <- "Mean Sq"

  typeIII[3]     <- typeIII["Sum of Sq"]/typeIII["Df"]
  # seq entry has originally 0 df

  # must equalize names else rbind() dosn't function

  names(typeIII)[2] <- "Sum Sq";

  names(typeI)[5]   <- "Pr(>F)"

  names(typeIII)[5] <- "Pr(>F)"

  if(!Group) {

    typeIII    <- rbind(typeIII[2:3,], # tmt, period

                        typeI["Seq",],

                        typeIII["Subj",],

                        typeI["Residuals",])

    # correct the sequence test, denominator is subject ms

    MSden <- typeIII["Subj","Mean Sq"]


    df2   <- typeIII["Subj","Df"]


    fval  <- typeIII["Seq","Mean Sq"]/MSden

    df1   <- typeIII["Seq","Df"]
    typeIII["Seq",4] <- fval

    typeIII["Seq",5] <- 1-pf(fval, df1, df2)

  } else {

    typeIII    <- rbind(typeIII[2:3,], # tmt, period

                        typeI["Group",],

                        typeI["Seq",],

                        typeI["Group:Seq",],

                        typeIII["Group:Seq:Subj",],

                        typeI["Residuals",])

    # correct the sequence test, denominator is subject ms

    MSden <- typeIII["Group:Seq:Subj","Mean Sq"]


    df2   <- typeIII["Group:Seq:Subj","Df"]
    

    fval_Seq  <- typeIII["Seq","Mean Sq"]/MSden

    df1_Seq   <- typeIII["Seq","Df"]    

    typeIII["Seq",4] <- fval_Seq

    typeIII["Seq",5] <- 1-pf(fval_Seq, df1_Seq, df2)    

    fval_Group  <- typeIII["Group","Mean Sq"]/MSden    

    df1_Group   <- typeIII["Group","Df"]


    typeIII["Group",4] <- fval_Group

    typeIII["Group",5] <- 1-pf(fval_Group, df1_Group, df2)    

    fval_Group_Seq  <- typeIII["Group:Seq","Mean Sq"]/MSden     

    df1_Group_Seq   <- typeIII["Group:Seq","Df"]


    typeIII["Group:Seq",4] <- fval_Group_Seq

    typeIII["Group:Seq",5] <- 1-pf(fval_Group_Seq, df1_Group_Seq, df2)      

  }

  

  attr(typeIII,"heading")<- attr(typeIII,"heading")[2:3] # suppress "single term deletion"

  attr(typeIII,"heading")[1] <- "(mimicing SAS/SPSS but with Seq tested against Subj)\n"

  # no need for the next

  mse      <- summary(muddle)$sigma^2

  

  # back transformation to the original domain

  CV <- 100*sqrt(exp(mse)-1)

  PE <- exp(lnPE)

  CI <- exp(lnCI)

  # output

  

  options(digits=8)

  cat("Type I sum of squares: ")

  print(typeI)

  cat("\nType III sum of squares: ")

  print(typeIII)

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

  cat("% confidence interval\n")

  cat("CV (%) ..................................:", 

      formatC(CV, format="f", digits=2),"\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")

  

  #reset options

  options(ow)

}



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)

Group2$Group <- 2

data<-rbind(Group1, Group2)



Analyse222BE(data, alpha=0.05, Group = FALSE)

Analyse222BE(data, alpha=0.05, Group = TRUE)

Now Group factor seems to be more suitable.
What bothers me is that Group effect is not equal in R and PHX, cannot figure out the reason

—
Kind regards,
Mittyri

Complete thread:

Potvin C in the EU Helmut 2013-04-16 17:40 [Two-Stage / GS Designs]
- 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 valuesmittyri 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 R mittyri 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