# load data
data(ozone)

# fit model 0 and obtain diagnostics
fit.ozone0<-ssanova(log10(upo3)~(ibtp+hmdt+ibht+dgpg+vsty)^2,
                    data=ozone,type="linear")
sum.ozone0<-summary(fit.ozone0,TRUE)
round(sum.ozone0$kappa,2)
round(sum.ozone0$cos,2)

# fit model 1 and obtain diagnostics
fit.ozone1<-ssanova(log10(upo3)~ibtp+hmdt+dgpg+vsty
                    +ibtp:(hmdt+vsty)+hmdt:dgpg,data=ozone)
sum.ozone1<-summary(fit.ozone1,TRUE)
round(sum.ozone1$kappa,2)
round(sum.ozone1$cos,2)

# fit model 2 and obtain diagnostics
fit.ozone2<-update(fit.ozone1,.~.+vdht+sbtp+wdsp-(ibtp+dgpg):hmdt)
sum.ozone2<-summary(fit.ozone2,TRUE)
round(sum.ozone2$kappa,2)
round(sum.ozone2$cos,2)

# fit model 3 and obtain diagnostics
fit.ozone3<-update(fit.ozone2,.~.-hmdt-vdht-wdsp)
sum.ozone3<-summary(fit.ozone3,TRUE)
round(sum.ozone3$kappa,2)
round(sum.ozone3$cos,2)
round(sum.ozone3$pi,2)

# fit model 4 and obtain diagnostics
fit.ozone4<-update(fit.ozone3,.~.-sbtp)
sum.ozone4<-summary(fit.ozone4,TRUE)
round(sum.ozone4$kappa,2)
round(sum.ozone4$cos,2)
round(sum.ozone4$pi,2)

# evaluate model 3 and model 4 for figure 3.7
est3.ibtp<-predict(fit.ozone3,ozone,se=TRUE,include="ibtp")
est4.ibtp<-predict(fit.ozone4,ozone,se=TRUE,include="ibtp")
est3.dgpg<-predict(fit.ozone3,ozone,se=TRUE,include="dgpg")
est4.dgpg<-predict(fit.ozone4,ozone,se=TRUE,include="dgpg")
grid<-seq(0,350,len=51)
est3.vsty<-predict(fit.ozone3,data.frame(vsty=grid),se=TRUE,include="vsty")
est4.vsty<-predict(fit.ozone4,data.frame(vsty=grid),se=TRUE,include="vsty")

# draw figure 3.7
par(mar=c(5,5,2,2)+.1,mfcol=c(2,3),mex=.6,pty="s")

ind<-order(ozone$ibtp)
plot(ozone$ibtp[ind],est3.ibtp$fit[ind],type="l",lty=1,ylim=c(-.65,.65),
xlab="inversion base temperature",ylab="",cex.lab=1.2)
lines(ozone$ibtp[ind],est3.ibtp$fit[ind]+1.96*est3.ibtp$se[ind],lty=2)
lines(ozone$ibtp[ind],est3.ibtp$fit[ind]-1.96*est3.ibtp$se[ind],lty=2)
rug(ozone$ibtp,ticksize=.05,lwd=.005)

plot(ozone$ibtp[ind],est4.ibtp$fit[ind],type="l",lty=1,ylim=c(-.65,.65),
xlab="inversion base temperature",ylab="",cex.lab=1.2)
lines(ozone$ibtp[ind],est4.ibtp$fit[ind]+1.96*est4.ibtp$se[ind],lty=2)
lines(ozone$ibtp[ind],est4.ibtp$fit[ind]-1.96*est4.ibtp$se[ind],lty=2)
rug(ozone$ibtp,ticksize=.05,lwd=.005)

ind<-order(ozone$dgpg)
plot(ozone$dgpg[ind],est3.dgpg$fit[ind],type="l",lty=1,ylim=c(-.6,.6),
xlab="Dagget pressure gradient",ylab="",cex.lab=1.2)
lines(ozone$dgpg[ind],est3.dgpg$fit[ind]+1.96*est3.dgpg$se[ind],lty=2)
lines(ozone$dgpg[ind],est3.dgpg$fit[ind]-1.96*est3.dgpg$se[ind],lty=2)
rug(ozone$dgpg,ticksize=.05,lwd=.005)

plot(ozone$dgpg[ind],est4.dgpg$fit[ind],type="l",lty=1,ylim=c(-.6,.6),
xlab="Dagget pressure gradient",ylab="",cex.lab=1.2)
lines(ozone$dgpg[ind],est4.dgpg$fit[ind]+1.96*est4.dgpg$se[ind],lty=2)
lines(ozone$dgpg[ind],est4.dgpg$fit[ind]-1.96*est4.dgpg$se[ind],lty=2)
rug(ozone$dgpg,ticksize=.05,lwd=.005)

grid<-seq(0,350,len=51)
plot(grid,est3.vsty$fit,type="l",lty=1,ylim=c(-.6,.6),
xlab="visibility",ylab="",cex.lab=1.2)
lines(grid,est3.vsty$fit+1.96*est3.vsty$se,lty=2)
lines(grid,est3.vsty$fit-1.96*est3.vsty$se,lty=2)
rug(jitter(ozone$vsty,factor=.4,amount=0),ticksize=.05,lwd=.005)

plot(grid,est4.vsty$fit,type="l",lty=1,ylim=c(-.6,.6),
xlab="visibility",ylab="",cex.lab=1.2)
lines(grid,est4.vsty$fit+1.96*est4.vsty$se,lty=2)
lines(grid,est4.vsty$fit-1.96*est4.vsty$se,lty=2)
rug(jitter(ozone$vsty,factor=.4,amount=0),ticksize=.05,lwd=.005)