# load data
data(LakeAcidity)

# code converting latitude and longitude to x-y coordinates
convert <- function(lat, lon) {
    lat <- lat/180*pi
    lon <- lon/180*pi
    m.lat <- (max(lat)+min(lat))/2
    m.lon <- (max(lon)+min(lon))/2
    x <- cos(m.lat)*sin(m.lon-lon)
    y <- sin(lat-m.lat)
    cbind(x,y)
}
sum(abs(convert(LakeAcidity$lat,LakeAcidity$lon)-LakeAcidity$geog))

# fit model with interaction and obtain diagnostics
fit.lake <- ssanova(ph~log(cal)*geog,data=LakeAcidity,
                    type="tp",order=2,method="m")
sum.lake <- summary(fit.lake,diag=TRUE)
round(sum.lake$kappa,2)
round(sum.lake$cos,2)

# fit additive model
fit.lake.add <- update(fit.lake,.~.-log(cal):geog)
# r-square and pi
sum.lake.add <- summary(fit.lake.add,diag=TRUE)
sum.lake.add$r.squared
sum.lake.add$pi

# evaluate the additive fit
est.cal <- predict(fit.lake.add,se=TRUE,inc="log(cal)",
                   new=model.frame(~log(cal),data=LakeAcidity))
grid <- seq(-.04,.04,len=31)
new <- model.frame(~geog,list(geog=cbind(rep(grid,31),
                   rep(grid,rep(31,31)))))
est.geog <- predict(fit.lake.add,new,se=TRUE,inc="geog")

# plot figure 4.5 top frame
ind<-order(LakeAcidity$cal)
plot(LakeAcidity$cal[ind],est.cal$fit[ind],log="x",type="l",lty=1,ylim=c(-.7,.9),
xlab="Calcium concentration",ylab="")
lines(LakeAcidity$cal[ind],(est.cal$fit+1.96*est.cal$se)[ind],lty=2)
lines(LakeAcidity$cal[ind],(est.cal$fit-1.96*est.cal$se)[ind],lty=2)
# convert x-y grid to longitude and latitude
m.lon <- (max(LakeAcidity$lon)+min(LakeAcidity$lon))/2
m.lat <- (max(LakeAcidity$lat)+min(LakeAcidity$lat))/2
grid.lon <- grid.lat <- (0:30)/375-.04
grid.lon <- -(m.lon-180*asin(grid.lon/cos(m.lat/180*pi))/pi)
grid.lat <- 180*asin(grid.lat)/pi+m.lat
# plot figure 4.5 left frame
par(pty="s")
contour(grid.lon,grid.lat,matrix(est.geog$fit,31,31),lty=3,labcex=.6,
        levels=c(-.25,-.15,-.05,.05,.15,.25,.35),
        xlab="longitude",ylab="latitude")
contour(grid.lon,grid.lat,matrix(est.geog$se,31,31),lty=2,labcex=.01,
        levels=c(.15),add=TRUE)
# plot figure 4.5 right frame
contour(grid.lon,grid.lat,matrix(est.geog$se,31,31),lty=3,labcex=.6,
        levels=c(.1,.15,.2,.3,.4,.5),xlab="longitude",ylab="latitude")
points(-LakeAcidity$lon,LakeAcidity$lat,cex=.6,col=5)


# use filled.contour
grid.lon <- grid.lat <- (0:30)/375-.04
m.lon <- (max(LakeAcidity$lon)+min(LakeAcidity$lon))/2
m.lat <- (max(LakeAcidity$lat)+min(LakeAcidity$lat))/2
tick.lon <- cos(m.lat/180*pi)*sin((m.lon-84:80)/180*pi)
tick.lat <- sin((33:37-m.lat)/180*pi)
# figure 4.5 left frame
filled.contour(grid.lon,grid.lat,matrix(est.geog$fit,31,31),asp=1,
               plot.title=title(xlab="longitude",ylab="latitude"),
               plot.axes={points(LakeAcidity$geog);
                          axis(1,at=tick.lon,lab=84:80);
                          axis(2,at=tick.lat,lab=33:37)})
# figure 4.5 right frame
filled.contour(grid.lon,grid.lat,matrix(est.geog$se,31,31),asp=1,
               plot.title=title(xlab="longitude",ylab="latitude"),
               plot.axes={points(LakeAcidity$geog);
                          axis(1,at=tick.lon,lab=84:80);
                          axis(2,at=tick.lat,lab=33:37)})