# # Epage 137 # # Written by: # -- # John L. Weatherwax 2009-04-21 # # email: wax@alum.mit.edu # # Please send comments and especially bug reports to the # above email address. # #----- # version$language == "R" for R version$language == NULL for SPlus if(is.null(version$language) == FALSE) require(alr3) else library(alr3) if(is.null(version$language) == FALSE) data(mile) attach(mile) # 6.18.1: # indsM <- mile$Gender== "Male" indsF <- mile$Gender== "Female" postscript("../../WriteUp/Graphics/Chapter6/prob_18_sp.eps", onefile=FALSE, horizontal=FALSE) plot( Year[indsM], Time[indsM], pch=c("m"), col=c("red"), ylim=range(Time[indsM],Time[indsF]), xlim=range(Year[indsM],Year[indsF]), xlab="Year", ylab="Time") par(new=T) plot( Year[indsF], Time[indsF], pch=c("f"), col=c("black"), ylim=range(Time[indsM],Time[indsF]), xlim=range(Year[indsM],Year[indsF]), xlab="Year", ylab="Time" ) par(new=F) dev.off() # 6.18.2: set this code up to do the model specifcation in general: # Xvar = mile$Year Fvar = factor(mile$Gender,ordered=FALSE) Yvar = mile$Time plot(Fvar,Yvar,xlab="factor index" ,ylab="Response") # fit the various four cases of models to this dateset: # m1 <- lm(Yvar~ Xvar + Fvar + Fvar:Xvar, na.action=na.omit) # this is model 1 the most general m2 <- lm(Yvar~ Xvar + Fvar , na.action=na.omit) # this is model 2 parallel m3 <- lm(Yvar~ Xvar + Fvar:Xvar , na.action=na.omit) # this is model 3 common intercept m4 <- lm(Yvar~ Xvar , na.action=na.omit) # this is model 4 the least general (all the same) # get some statistics on m1: # print( c(m1$df.residual, sum( m1$residuals^2 )) ) # compare all models to the most general model (model 1): # a2<-anova(m2,m1) print( c( a2$Res.Df[1], a2$RSS[1], a2$F[2], a2$'Pr(>F)'[2] ), digits=3 ) a3<-anova(m3,m1) print( c( a3$Res.Df[1], a3$RSS[1], a3$F[2], a3$'Pr(>F)'[2] ), digits=3 ) a4<-anova(m4,m1) print( c( a4$Res.Df[1], a4$RSS[1], a4$F[2], a4$'Pr(>F)'[2] ), digits=3 ) # compute various linear models depending: # mile$D <- Fvar # this is also the most general model but specified so that each factor has a different slope and itercept n1 <- lm(Yvar~ -1+D+D:Xvar, mile,na.action=na.omit) # the coefficients from this model "n1" are given by: # #> coef(n1) # b0 b1 b2 b3 # DFemale DMale DFemale:Xvar DMale:Xvar #2309.4247477 953.7469611 -1.0336960 -0.3661867 delta.method(n1, "(240 - b0)/b2") # lets check that I got the coefficients correct by computing the point estimate of \hat{y} by hand: # yHat = ( 240 - coef(n1)[1] ) / coef(n1)[3] # Yes! # lets just fit the data directy: indsF <- mile$Gender == "Female" Xvar <- mile$Year[indsF] Yvar <- mile$Time[indsF] m0 <- lm( Yvar ~ Xvar ) summary( m0 ) # use this new linear model to determine the standard error ... it is the same! delta.method(m0, "(240 - b0)/b1") # 6.18.4: # delta.method(n1, "(b0 - b1)/(b3 - b2)") # lets check that I got the coefficients correct by computing the point estimate of \hat{y} by hand: # yHat = ( coef(n1)[1] - coef(n1)[2] )/( coef(n1)[4] - coef(n1)[3] ) print(yHat)