# this is the file script for Cashmore water % data from Lark et al (005)
# first we just do some playing around with the data to look at the sampling points and trend
# read it in from .csv file
#
# cols are  x y soilwater xsca ysca x2 y2 xy
#
cashmore <- read.table('cashmore.csv',header=T,sep=',')
# three potential plots
#figures 1 and 2
#
#
trellis.device(color=F)
xyplot(ysca~xsca,data=cashmore,ylab='y (m) /100',xlab='x (m) /100')
par(mfrow=c(1,2))
plot(cashmore$ysca,cashmore$soilwater,ylab='Water Content %',xlab='y (m)/100')
plot(cashmore$xsca,cashmore$soilwater,ylab='Water Content %',xlab='x (m)/100')

#
# clear linear trend for both x and y
# fit some prelimiary models and do variograms for them
models.out <- list()
models.out[[1]] <- samm(soilwater~1,data=cashmore)
models.out[[2]] <- samm(soilwater~1+xsca+ysca,data=cashmore)
models.out[[3]] <- samm(soilwater~1+xsca+ysca+x2+y2+xy,data=cashmore)
#> summary(models.out[[1]])$coef
#            solution std error  z ratio
#(Intercept)   22.275 0.3494769 63.73812
#> summary(models.out[[2]])$coef
#             solution std error    z ratio
#ysca        -5.933036 0.3145672 -18.860949
#xsca        -1.730776 0.2004775  -8.633267
#(Intercept) 64.550098 2.1159049  30.507088
#> summary(models.out[[3]])$coef
#               solution  std error    z ratio
#xy            1.3074717  0.4327023  3.0216427
#y2            3.9587296  0.5817520  6.8048410
#x2            0.1813361  0.2243806  0.8081633
##xsca        -11.1076369  3.4730860 -3.1982038
#(Intercept) 220.2671939 23.7761887  9.2641927
#>
dvrgs <- list()
vrgs <-list()
dvrgs[[1]] <- samm.variogram(x=cashmore$x, y=cashmore$y,
        z=residuals(models.out[[1]]),angle=c(0,45,90,135),angle.tol=45,grid=F)
vrgs[[1]] <- samm.variogram(x=cashmore$x, y=cashmore$y,z=residuals(models.out[[1]]),grid=F)
dvrgs[[2]] <- samm.variogram(x=cashmore$x, y=cashmore$y,
        z=residuals(models.out[[2]]),angle=c(0,45,90,135),angle.tol=45,grid=F)
vrgs[[2]] <- samm.variogram(x=cashmore$x, y=cashmore$y,z=residuals(models.out[[2]]),grid=F)
dvrgs[[3]] <- samm.variogram(x=cashmore$x, y=cashmore$y,
        z=residuals(models.out[[3]]),angle=c(0,45,90,135),angle.tol=45,grid=F)
vrgs[[3]] <- samm.variogram(x=cashmore$x, y=cashmore$y,z=residuals(models.out[[3]]),grid=F)

vario.df <- rbind(dvrgs[[1]],vrgs[[1]],dvrgs[[2]],vrgs[[2]],dvrgs[[3]],vrgs[[3]])
vario.df$type <- factor(rep(rep(c('0','45','90','135','omni'),c(table(dvrgs[[1]]$angle),dim(vrgs[[1]])[1])),3),
    levels=c('0','45','90','135','omni'))
vario.df$model <- factor(rep(c('NULL','LIN','QUAD'),each=(dim(dvrgs[[1]])[1]+dim(vrgs[[1]])[1])),levels=c('NULL','LIN','QUAD'))

####################
# figure 3 variograms for each of three trend models
###################
#

trellis.device(color=F)
sps <- trellis.par.get("superpose.symbol")
sps$pch <- 1:7
spl <- trellis.par.get("plot.line")
spl$lty <- 1:7
spl$type <- "b"
spl$pch <- 1:5
#spl$col <- sps$col
spl$col <- 1
sps$col <- 1
trellis.par.set("superpose.symbol",sps)
trellis.par.get("superpose.symbol")
trellis.par.set("superpose.line",spl)
trellis.par.get("superpose.line")
xyplot(gamma~distance|model, data=vario.df,groups=type,
    xlab="Distance",ylab='Semi-variance',panel=panel.superpose,type="b",
    scales=list(  x=list(relation="free",limits=c(0,max(vario.df$distance))) , 
        y=list( relation="free", limits=list(
        c(0,max(dvrgs[[1]]$gamma)),c(0,max(dvrgs[[2]]$gamma)),
        c(0,max(dvrgs[[3]]$gamma)) )  )  ),
    as.table=TRUE,
    key=list(columns=1,
        text=list(paste(c("   ", "  ", "  ", " ", ""),
                    unique(vario.df$type) )),
        lines=Rows(spl,1:5),    # next line previously had y=0.95
                                # changed from y=0.93 to 0.85 with main title
        x=0.75,y=0.40,corner=c(0,1),title="Angle (degrees)",cex.title=1,divide=1.5  )
    )

######################
# figure 4 correlation contours for ani- mat
######################
par(mfrow=c(1,1))
rx <- range(cashmore$xsca)
rx <- (rx[2]-rx[1])/2
ry <- range(cashmore$ysca)
ry <- (ry[2] - ry[1])/2

h1 <- seq(-rx,rx,2*rx/100)
h2 <- seq(-ry,ry,2*ry/100)

res <- outer(h1,h2,function(x,y,...) my.matern(x,y,...),zeta=.187,nu=.5,alpha=.699,delta=4.04)
contour(h1,h2,res,xlim=c(-1,1)/1.6,ylim=c(-1,1)/1.6,levels=seq(.1,.9,.1))
mtext("Anisotropic - (4.04, 40deg)",cex=1.0)
#
######################################
# now the modelling in samm. I have done model selection in ASReml and the model is exponential
######################################
models.out[[4]] <- samm(soilwater~1+xsca+ysca+x2+y2+xy,random=~units,rcov=~iexp(xsca,ysca),data=cashmore)
models.out[[5]] <- samm(soilwater~1+xsca+ysca+x2+y2+xy,random=~units,rcov=~ieuc(xsca,ysca),data=cashmore)
models.out[[6]] <- samm(soilwater~1+xsca+ysca+x2+y2+xy,random=~units,rcov=~aexp(xsca,ysca),data=cashmore)

h1 <- outer(cashmore$xsca,cashmore$xsca,function(x,y) (x-y)^2)
h2 <- outer(cashmore$ysca,cashmore$ysca,function(x,y) (x-y)^2)
deuc <- sqrt(h1+h2)
h1 <- outer(cashmore$xsca,cashmore$xsca,function(x,y) abs(x-y))
h2 <- outer(cashmore$ysca,cashmore$ysca,function(x,y) abs(x-y))
dman <- h1+h2
dman <- as.vector(dman)
deuc <- as.vector(deuc)
vveuc <- .473 + 1.426*(1- .00487^deuc)
vvman <- .694 + 13.39*(1-.763^dman)

#
##################
# plots to show x/y effects
#################
temp <- cashmore
temp$res <- residuals(models.out[[5]])
temp$res1 <- residuals(models.out[[6]])
temp <- temp[order(temp$ysca,temp$xsca),]
xyplot(res~xsca,data=temp[temp$ysca==5.5,],type='b')
xyplot(res~ysca,data=temp[temp$xsca==6.5,],type='b')
xyplot(soilwater~ysca,data=temp[temp$xsca==6.5,],type='b')
xyplot(res~xsca|ysca,data=temp,type='b')


###############################
#
#
#
# now for the plots from the analysed data. I have two runs, one both with nu=1, but isotropic and
# one anisotropic
# yht files contain the hat matrix and the fitted values for plotting
#
##############################
xuniq<-sort(unique(cashmore$xsca))
yuniq<-sort(unique(cashmore$ysca))
datagrid<-expand.grid(xsca=xuniq,ysca=yuniq)
# first a check of my.krige
# it gives the right answers
#
X <- cbind(cashmore$xy,cashmore$y2,cashmore$x2,cashmore$ysca,cashmore$xsca,1)
Xp <- X
y <- cashmore$soilwater
locations <- cbind(cashmore$xsca,cashmore$ysca)
locations <- rbind(locations,locations)
out.pred <- my.krige(out=models.out[[6]],Xp=Xp, X=X, y=y, locations = locations, phi=models.out[[6]]$gammas[3:4])
Xp <- cbind(datagrid$xsca*datagrid$ysca,datagrid$ysca^2, datagrid$xsca^2, datagrid$ysca, datagrid$xsca,1) 
locations <- rbind(cbind(cashmore$xsca,cashmore$ysca),cbind(datagrid$xsca,datagrid$ysca))

out.aexp <- my.krige(out=models.out[[6]],Xp=Xp, X=X, y=y, locations = locations)
out.ieuc <- my.krige(out=models.out[[5]],Xp=Xp, X=X, y=y, locations = locations, phi=models.out[[5]]$gammas[3],model='ieuc')

##########################################
##########################################
# figure 5 image plots of the predictions and the mseps
# one for the notes and one for the book
#
##########################################

par(mfrow=c(2,2))
image(xuniq,yuniq,matrix(as.vector(out.aexp$ypred),34,23),col=terrain.colors(50),xlim=c(5.0,7.5),ylim=c(4.2,6.00),ylab='y',
xlab='x')
mtext('Water Content: Ani-CB')
contour(xuniq,yuniq,matrix(as.vector(out.aexp$ypred),34,23),levels=seq(5,40,1),
        add=TRUE,col="peru")

image(xuniq,yuniq,matrix(as.vector(out.aexp$pev),34,23),col=terrain.colors(50),xlim=c(5.0,7.5),ylim=c(4.2,6.00),ylab='y',
xlab='x')
mtext("MSEP: Ani-CB")
contour(xuniq,yuniq,matrix(as.vector(out.aexp$pev),34,23),levels=seq(0,4,.1),
        add=TRUE,col="peru")

image(xuniq,yuniq,matrix(as.vector(out.ieuc$ypred),34,23),col=terrain.colors(50),xlim=c(5.0,7.5),ylim=c(4.2,6.00),ylab='y',
xlab='x')
mtext('Water Content: Iso-Euc')
contour(xuniq,yuniq,matrix(as.vector(out.ieuc$ypred),34,23),levels=seq(5,40,1),
        add=TRUE,col="peru")

image(xuniq,yuniq,matrix(as.vector(out.ieuc$pev),34,23),col=terrain.colors(100),xlim=c(5.0,7.5),ylim=c(4.2,6.00),ylab='y',
xlab='x')
mtext("MSEP: Iso-Euc")
contour(xuniq,yuniq,matrix(as.vector(out.ieuc$pev),34,23),levels=seq(0,4,.1),
        add=TRUE,col="peru")
        
###############
# figure 5 for the book
#
##############

par(mfrow=c(2,2))
contour(xuniq,yuniq,matrix(as.vector(out.aexp$ypred),34,23),levels=seq(5,40,1),
        add=FALSE,col="black")
mtext('Water Content: Ani-CB')

contour(xuniq,yuniq,matrix(as.vector(out.aexp$pev),34,23),levels=seq(0,4,.1),
        add=FALSE,col="black")
mtext("MSEP: Ani-CB")

contour(xuniq,yuniq,matrix(as.vector(out.ieuc$ypred),34,23),levels=seq(5,40,1),
        add=FALSE,col="black")
mtext('Water Content: Iso-Euc')

contour(xuniq,yuniq,matrix(as.vector(out.ieuc$pev),34,23),levels=seq(0,4,.1),
        add=FALSE,col="black")
mtext("MSEP: Iso-Euc")

############
# comparisons of the predictions and pevs from the two models
#
# figure 6
###########
par(mfrow=c(1,2))
plot(out.aexp$ypred,out.ieuc$ypred,ylab='Iso-Euc',xlab='Ani-CB')
mtext('Predictions')
identify(out.aexp$ypred,out.ieuc$ypred,datagrid$ysca)
plot(out.aexp$pev,out.ieuc$pev,ylab='Iso-Euc',xlab='Ani-CB')
mtext('MSEPs')
identify(out.aexp$pev,out.ieuc$pev,datagrid$ysca)