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

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) <- "Mean Sq" typeIII <- typeIII["Sum of Sq"]/typeIII["Df"] # seq entry has originally 0 df # must equalize names else rbind() dosn't function names(typeIII) <- "Sum Sq"; names(typeI) <- "Pr(>F)" names(typeIII) <- "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") <- "(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, format="f", digits=2) ,"\n")   cat("Upper confidence limit.(%)...............:",       formatC(100*CI,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 Ing. Helmut Schütz 