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

posted by mittyri – Russia, 2016-10-10 18:33 (1396 d 23:27 ago) – Posting: # 16710
Views: 29,190

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,977 posts in 4,374 threads, 1,460 registered users;
online 17 (0 registered, 17 guests [including 9 identified bots]).
Forum time: Friday 18:01 CEST (Europe/Vienna)

When someone says his conclusions are objective,
he means that they are based on prejudices
which many other people share.    Celia Green

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