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 o=lm(y~time+temp+time:temp,data=d) model.matrix(o) coef(o) vcov(o) #Cell means are # # temp 20 temp 30 # ------------------------- # time 3 mu mu+temp30 # time 6 mu+time6 mu+time6+temp30+time6:temp30 # #Time main effects? # # (mu+mu+temp30)/2 # -(mu+time6+mu+time6+temp30+time6:temp30)/2 # ---------------------------- # -time6-time6:temp30/2 # # H0:time6+time6:temp30/2=0 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) } o\$coe Ctime=matrix(c( 0,1,0,.5 ),nrow=1,byrow=T) Ctemp=matrix(c( 0,0,1,.5 ),nrow=1,byrow=T) Ctimetempint=matrix(c( 0,0,0,1 ),nrow=1,byrow=T) test(o,Ctime) test(o,Ctemp) test(o,Ctimetempint) #The R function anova will produce tests for #the presence of time main effects, #temp main effects, and time-by-temp interaction. #However, these are "Type I Tests" that do not #in general match the "Type III Tests" above. anova(o) #Any difference among the four treatment means? # # temp 20 temp 30 # ------------------------- # time 3 mu mu+temp30 # time 6 mu+time6 mu+time6+temp30+time6:temp30 Coverall=matrix(c( 0,1,0,0, 0,0,1,0, 0,0,0,1 ),nrow=3,byrow=T) test(o,Coverall) #Other choices for C can provide the same test. Coverall=matrix(c( 0,0,1,0, 0,0,1,1, 0,1,0,0 ),nrow=3,byrow=T) test(o,Coverall) #Any difference between the 3 month, 20 degree mean #and the average of the other three treatment means? # # temp 20 temp 30 # ------------------------- # time 3 mu mu+temp30 # time 6 mu+time6 mu+time6+temp30+time6:temp30 # # # 2*time6/3+2*temp30/3+time6:temp30/3 C=matrix(c( 0,2/3,2/3,1/3 ),nrow=1,byrow=T) test(o,C) #The following function can be used to obtain #confidence intervals for each element of an #estimable C*beta. ci=function(lmout,C,a=0.05) { b=coef(lmout) V=vcov(lmout) df=lmout\$df Cb=C%*%b se=sqrt(diag(C%*%V%*%t(C))) tval=qt(1-a/2,df) low=Cb-tval*se up=Cb+tval*se m=cbind(C,Cb,se,low,up) dimnames(m)[[2]]=c(paste("c",1:ncol(C),sep=""), "estimate","se", paste(100*(1-a),"% Conf.",sep=""), "limits") m } ci(o,Ctime) ci(o,Coverall) summary(o) ################################################## # temp 20 temp 30 # ------------------------- # time 3 mu mu+temp30 # time 6 mu+time6 mu+time6+temp30+time6:temp30 ##################################################