# # Epage 193 # # 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) } data(lakes) attach(lakes) summary(lakes) lakes <- na.omit( lakes ) # remove NA values: lakes$Elev <- abs(lakes$Elev) # there is only one negative value in lakes$Elev and that seems to be very different than the others pairs(~MaxDepth+MeanDepth+Cond+Elev+Lat+Long+Dist+NLakes+Photo+Area+Species ) myFrame <- data.frame( lakes[,c(1,2,3,4,5,6,7,8,9,10,11)] ) cor(myFrame)[1,] # based on the correlation matrix we'll drop MeanDepth (0.03408282) which also should be very correlated to MaxDepth #pairs(~MaxDepth+Cond+Elev+Lat+Long+Dist+NLakes+Photo+Area+Species ) # alternativly consider the three variables that are most highly correlated with Species sort(abs(cor(myFrame)[1,])) # gives: Dist Elev Area # 0.35328247 0.40441612 0.74810087 # pairs(~Elev+Dist+Area+Species ) # pick values such that the independent variables are transformed to a set of predictors as close to # linearly related as possible # ans <- bctrans(Species~Elev+Dist+Area, data=lakes) summary(ans) # transform Dist and Area with logarithmic transformations: # postscript("../../WriteUp/Graphics/Chapter8/prob_12_pairs_plot.eps", onefile=FALSE, horizontal=FALSE) pairs( ~log(Elev)+log(Dist)+log(Area)+Species ) dev.off() lakesTrans <- powtran( ans, lambda=c(0,0,0) ) m0 <- lm( Species ~ logElev+logDist+logArea, data=lakesTrans ) inv.tran.plot( lakesTrans$Species, predict(m0) ) # this indicates no transformation of the response is needed. summary(m0) # based on this and the pairs plot delete the variable log(Elev) m1 <- lm( Species ~ logDist+logArea, data=lakesTrans ) summary(m1) # lets plot a residual plot (constant variance does not seem to be violated): plot( predict(m1), residuals(m1,type="pearson") ) postscript("../../WriteUp/Graphics/Chapter8/prob_12_res_plots.eps", onefile=FALSE, horizontal=FALSE) residual.plots( m1 ) dev.off() # lets do a visual assessment: mmps( m1, sd=TRUE )