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:
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
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
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
- FDA group model in R ElMaestro 2016-10-17 18:58
- FDA group model in R mittyri 2016-10-11 13:43
- Fixed effects model: changing the F, p values ElMaestro 2016-10-10 20:10
- Fixed effects model with Group term VStus 2016-10-20 13:54
- Fixed effects model: changing the F, p valuesmittyri 2016-10-10 18:33
- Fixed effects model with Group term mittyri 2016-10-09 13:27
- Group effects EMA mittyri 2016-10-07 15:06
- Group effects EMA ElMaestro 2016-10-07 12:45
- Group effects EMA VStus 2016-10-06 22:06
- Group effects FDA/EMA zizou 2016-12-31 02:54
- Group effects FDA/EMA Helmut 2017-06-03 14:57
- Group effects EMA mittyri 2014-10-01 11:25
- Group effects FDA/EMA ElMaestro 2014-08-21 17:31
- Group effects FDA/EMA mittyri 2014-08-21 15:34
- Group effects FDA/EMA Helmut 2013-04-19 14:34
- 2 Groups model FDA hiren379 2013-08-22 14:32
- Homework Helmut 2013-08-23 11:53
- Group effects obsolete? d_labes 2013-04-18 09:55
- 2 Groups model FDA Helmut 2013-04-17 16:37
- 2 Groups model FDA d_labes 2013-04-17 13:10
- Potvin C in the EU Helmut 2013-04-17 12:23
- Potvin C in the EU ElMaestro 2013-04-17 02:31