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

posted by mittyri – Russia, 2016-10-10 18:33  – Posting: # 16710
Views: 26,951

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. MSEsubject:

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

Activity
 Admin contact
20,145 posts in 4,248 threads, 1,385 registered users;
online 5 (0 registered, 5 guests [including 2 identified bots]).
Forum time (Europe/Vienna): 14:51 CET

A little Learning is a dang’rous Thing;
Drink deep, or taste not the Pierian Spring:
There shallow Draughts intoxicate the Brain,
And drinking largely sobers us again.    Alexander Pope

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