# # 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) } # load the data and append decision variables to each data frame: data(forbes) forbes$D <- rep(1,length(forbes$Pressure)) data(hooker) hooker$Lpres <- 100 * log( hooker$Pressure, 10 ) hooker$D <- rep(2,length(hooker$Temp)) # create a total data frame: allTemps <- c(forbes$Temp, hooker$Temp) allLpres <- c(forbes$Lpres, hooker$Lpres) allClass <- c(forbes$D, hooker$D) tempd <- data.frame( Temp=allTemps, Lpres=allLpres, D=allClass ) # create a factor from D: tempd$DF <- factor( tempd$D, ordered=FALSE ) # set this code up to do the model specifcation in general: # Xvar = tempd$Temp Fvar = tempd$DF Yvar = tempd$Lpres #postscript("../../WriteUp/Graphics/Chapter6/prob_5_boxplot.eps", onefile=FALSE, horizontal=FALSE) plot(Fvar,Yvar,xlab="factor index" ,ylab="Response") #dev.off() # 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 ) # Make a figure with lattice fitting all four models. # library(lattice) # this is needed for the following lattice plot ... NC <- length(unique(Fvar)) # the number of classes n <- length(Fvar) # the total number of observations # create a data frame with the needed information for the plot: # tmpDataFrame <- tempd tmpDataFrame$D <- Fvar tmpDataFrame$Xvar <- Xvar tmpDataFrame$Yvar <- Yvar # compute various linear models depending: # n1 <- lm(Yvar~ -1+D+D:Xvar, tmpDataFrame,na.action=na.omit) # model #1: this is the most general model n2 <- lm(Yvar~ -1+D+Xvar, tmpDataFrame,na.action=na.omit) # model #2: n3 <- lm(Yvar~ +1+D:Xvar, tmpDataFrame,na.action=na.omit) # model #3 n4 <- lm(Yvar~ Xvar, tmpDataFrame,na.action=na.omit) # model #4 postscript("../../WriteUp/Graphics/Chapter6/prob_5_model_comparison.eps", onefile=FALSE, horizontal=FALSE) if(is.null(version$language) == FALSE) { d <- rbind(tmpDataFrame,tmpDataFrame,tmpDataFrame,tmpDataFrame) labels <- c("(a) General", "(b) Parallel","(c) Common intercept", "(d) Common regression") d$f <- rep(c(1,2,3,4), rep(n,4)) d$fig <- rep(labels, rep(n,4)) xyplot(Yvar~Xvar|fig,group=D, d, subscripts=TRUE,ID=d$f, as.table=TRUE, key=simpleKey(as.character(1:NC),FALSE,FALSE,T), xlab=expression(Xvar), cex=.75, between=list(x=.32,y=.32), panel = function(x,y,groups,subscripts,ID){ panel.superpose(x,y,groups=groups,subscripts=subscripts, pch=as.character(1:NC)) g <- unique(ID[subscripts]) # which panel if(g==1) {panel.abline(coef(n1)[c(1,3)],lty=1) panel.abline(coef(n1)[c(2,4)],lty=2)} if(g==3) {panel.abline(coef(n3)[c(1,2)],lty=1) panel.abline(coef(n3)[c(1,3)],lty=2)} if(g==2) {panel.abline(coef(n2)[c(1,3)],lty=1) panel.abline(coef(n2)[c(2,3)],lty=2)} if(g==4) {panel.abline(coef(n4)[c(1,2)],lty=1)} }) } dev.off()