# NOT RUN {
str(case1402)
attach(case1402)
## EXPLORATION AND MODEL DEVELOPMENT; FORREST CULTIVAR
logForrest <- log(Forrest)
# Fit model without interactions first--to examine partial residual plots.
myLm1 <- lm(logForrest ~ Stress + SO2 + O3)
if(require(car)){ # Use the car library
crPlots(myLm1) # Partial res plots => linear effects of SO2 and O3 look ok.
}
myLm2 <- lm(logForrest ~ (Stress + SO2 + O3)^2) # all 2-factor interactions.
plot(myLm1, which=1) # Residual plot looks ok.
anova(myLm1,myLm2) # Test for interactive effects.
## INFERENCE AND INTERPRETATION; FORREST CULTIVAR
summary(myLm1)
betaF <- myLm1$coef
# Effect of 0.01 increase in SO2 (note coeff is negative):
100*(1 - exp(0.01*betaF[3]))
#1.655701; a 1.65% decrease in median yield
100*(1-exp(0.01*confint(myLm1,"SO2")))
#3.772277 -0.5074294: 95% onfidence interval for effect of 0.01 increase in SO2
# Effect of 0.01 increase in O3 (note coeff is negative):
100*(1 - exp(0.01*betaF[4]))
# 5.585979 I.e. a 5.6% decrease in median yield
100*(1-exp(0.01*confint(myLm1,"O3")))
#7.445966 3.688613: 95% confidence interval for effect of 0.01 increase in O3
# Effect of water stress (note coeff is positive for effect of well-watered):
100*(1 - exp(-betaF[2])) # Get estimate for negative of this beta
#3.220556: a 3.2% decrease in median yield due to water stress
100*(1-exp(-confint(myLm1,2)))
#-7.875521 13.17529: 95% confidence interval
## DISPLAY FOR PRESENTATION; FORREST CULTIVAR
jO3 <- jitter(O3,factor=.25) # Jitter for plotting; jittering factor 0.25.
jS <- jitter(SO2,factor=.25) # Jitter SO2 for plotting.
cexfac <- 1.75 # Use character expansion factor of 1.75 for plotting symbols.
opar <- par(no.readonly=TRUE) # Store current graphics parameters settings
par(mfrow=c(1,2)) # Make a panel of 2 plots in 1 row.
myPointCode <- ifelse(Stress=="Well-watered",21,24)
myPointColor <- ifelse(Stress=="Well-watered","green","orange")
par(mar=c(4.1,4.1,2.1,0.1))
plot(Forrest ~ jO3, log="y", ylab="Yield of Forrest Cultivar (kg/ha)",
xlab=expression(paste(italic("Ozone ("),mu,"L/L), jittered")),
pch=myPointCode, lwd=2, bg=myPointColor, cex=cexfac)
legend(.02,2400, c("Well-watered","Water Stressed"), pch=c(21,24),
pt.cex=cexfac, pt.bg=c("green","orange"), pt.lwd=2, lty=c(3,1), lwd=c(2,2))
dummyO <- seq(min(O3), max(O3), length=2)
curve1 <- exp(betaF[1] + betaF[3]*mean(SO2) + betaF[4]*dummyO)
curve2 <- exp(betaF[1] + betaF[2] + betaF[3]*mean(SO2)+ betaF[4]*dummyO)
lines(curve1 ~ dummyO,lwd=2)
lines(curve2 ~ dummyO,lwd=2,lty=3)
# PLOT FORREST VS SO2
par(mar=c(4.1,2.1,2.1,2.1)) # Change margins
plot(Forrest ~ jS, log="y", ylab="",
xlab=expression(paste(italic("Sulfur Dioxide ("),mu,"L/L), jittered")),
yaxt="n", pch=myPointCode, lwd=2, bg=myPointColor, cex=cexfac)
dummyS <- seq(min(SO2),max(SO2),length=2)
curve1 <- exp(betaF[1] + betaF[3]*dummyS + betaF[4]*mean(O3))
curve2 <- exp(betaF[1] + betaF[2] + betaF[3]*dummyS + betaF[4]*mean(O3))
lines(curve1 ~ dummyS,lwd=2)
lines(curve2 ~ dummyS,lwd=2,lty=3)
par(opar) # Restore previous graphics parameter settings
detach(case1402)
# }
Run the code above in your browser using DataLab