time=factor(rep(c(3,6),each=5))
temp=factor(rep(c(20,30,20,30),c(2,3,4,1)))
y=c(2,5,9,12,15,6,6,7,7,16)
d=data.frame(time,temp,y)
d
full=lm(y~time+temp+time:temp,data=d)
model.matrix(full)
coef(full)
##################################################
# temp 20 temp 30
# -------------------------
# time 3 mu mu+temp30
#
# time 6 mu+time6 mu+time6+temp30+time6:temp30
#
##################################################
test=function(lmout,C,d=0){
b=coef(lmout)
V=vcov(lmout)
dfn=nrow(C)
dfd=lmout$df
Cb.d=C%*%b-d
Fstat=drop(t(Cb.d)%*%solve(C%*%V%*%t(C))%*%Cb.d/dfn)
pvalue=1-pf(Fstat,dfn,dfd)
list(Fstat=Fstat,pvalue=pvalue)
}
Coverall=matrix(c(
0,1,0,0,
0,0,1,0,
0,0,0,1
),nrow=3,byrow=T)
test(full,Coverall)
reduced=lm(y~1,data=d)
model.matrix(reduced)
rvsf=function(reduced,full)
{
sser=deviance(reduced)
ssef=deviance(full)
dfer=reduced$df
dfef=full$df
dfn=dfer-dfef
Fstat=(sser-ssef)/dfn/
(ssef/dfef)
pvalue=1-pf(Fstat,dfn,dfef)
list(Fstat=Fstat,dfn=dfn,dfd=dfef,pvalue=pvalue)
}
rvsf(reduced,full)
anova(reduced,full)
Cinteraction=matrix(c(
0,0,0,1
),nrow=1,byrow=T)
test(full,Cinteraction)
anova(full)
summary(full)
additive=lm(y~time+temp,data=d)
model.matrix(additive)
coef(additive)
####################################
# temp 20 temp 30
# -------------------------
# time 3 mu mu+temp30
#
# time 6 mu+time6 mu+time6+temp30
#
####################################
rvsf(additive,full)
anova(additive,full)
drop1(full,test="F")
anova(additive)
drop1(additive,test="F")
anova(lm(y~temp+time,data=d))
#Previously we saw how to test for time main effects
#and temp main effects in the full model by testing
#H0: C beta = d.
#
#It is possible but not as easy to test for these main
#effects using the reduced versus full model approach.
#
#We will use the test for time main effects as an
#example.
#
#We need to find a matrix whose column space allows
#for temp main effects and time-by-temp interaction
#but no time main effects.
#
#It is natural to try the following model
#specification.
o=lm(y~temp+time:temp,data=d)
model.matrix(o)
#Examination of this design matrix shows that the
#cell means are modeled as
#
# temp 20 temp 30
# --------------------------------
# time 3 mu mu+temp30
#
# time 6 mu+temp20:time6 mu+temp30+temp30:time6
#
#It is easy to see that this is just the full model
#in which each treatment group is allowed to have
#its own mean. Thus, we can't use this code to
#obtain our reduced model fit.
#One way to obtain the test for time main effects
#by comparing a reduced and full model is as follows.
full=lm(y~time+temp+time:temp,data=d)
C=matrix(c(
0,1,0,.5
),nrow=1,byrow=T)
B=matrix(c(
1,0,0,0,
0,0,1,0,
0,0,0,1,
0,1,0,.5
),nrow=4,byrow=T)
#Note that X beta = X B^{-1} B beta.
#
#Let W = X B^{-1} and alpha = B beta.
#
#Then C beta = 0 is equivalent to alpha_4=0.
W=model.matrix(full)%*%solve(B)
W0=W[,1:3]
newfull=lm(y~W-1)
newreduced=lm(y~W0-1)
anova(newreduced,newfull)
rvsf(newreduced,newfull)
test(full,C)