###################################################
### chunk number 1: libs
###################################################
options(width = 65) 
library(spBayes)
library(MBA)
library(geoR)
library(fields)
library(rgl)
library(sp)
library(maptools)
library(rgdal)
library(classInt)
library(gstat)
library(lattice)
library(xtable)


###################################################
### chunk number 2: coords
###################################################
data(BEF.dat)

BEF.dat <- BEF.dat[BEF.dat$ALLBIO02_KGH>0,]
bio <- BEF.dat$ALLBIO02_KGH*0.001;
log.bio <- log(bio)

coords <- as.matrix(BEF.dat[,c("XUTM","YUTM")])

plot(coords, pch=19, cex=0.5, xlab="Easting (m)", ylab="Northing (m)")


###################################################
### chunk number 3: MBASurf
###################################################
x.res <- 100; y.res <- 100

surf <- mba.surf(cbind(coords, bio), no.X=x.res, no.Y=y.res, h=5, m=2, extend=FALSE)$xyz.est
image.plot(surf, xaxs = "r", yaxs = "r", xlab="Easting (m)", ylab="Northing (m)")
points(coords)


###################################################
### chunk number 4: rglMBASurf eval=FALSE
###################################################
## col <- rbind(0, cbind(matrix(drape.color(surf[[3]]), x.res-1, y.res-1), 0))
## 
## surface3d(surf[[1]], surf[[2]], surf[[3]], col=col)
## 
## axes3d(); title3d(main="Biomass", xlab="Easting (m)", ylab="Northing (m)", zlab="Log metric tons of biomass")


###################################################
### chunk number 5: MBACovarSurf
###################################################
lm.bio <- lm(log.bio~ELEV+SLOPE+SUM_02_TC1+SUM_02_TC2+SUM_02_TC3, data=BEF.dat)
summary(lm.bio)

bio.resid <- resid(lm.bio)

par(mfrow=c(2,3))

surf <- mba.surf(cbind(coords, bio.resid), no.X=x.res, no.Y=y.res, h=5, m=2, extend=FALSE)$xyz.est
image.plot(surf, xaxs = "r", yaxs = "r", main = "Log metric tons of biomass residuals", xlab="Easting (m)", ylab="Northing (m)")

covars <- c("ELEV","SLOPE","SUM_02_TC1","SUM_02_TC2","SUM_02_TC3")

for(i in 1:5){
  surf <- mba.surf(cbind(coords, BEF.dat[,covars[i]]), no.X=x.res, no.Y=y.res, extend=FALSE)$xyz.est
  image.plot(surf, xaxs = "r", yaxs = "r", main=covars[i])
}


###################################################
### chunk number 6: geoRVariog
###################################################
max.dist <- 0.5*max(iDist(coords))
bins <- 20

vario.bio <- variog(coords=coords, data=log.bio, uvec=(seq(0, max.dist, length=bins)))

fit.bio <-variofit(vario.bio, ini.cov.pars=c(0.08, 500/-log(0.05)), ##sigma^2 and 1/phi 
                         cov.model="exponential", minimisation.function="nls", weights="equal")

vario.bio.resid <- variog(coords=coords, data=bio.resid, uvec=(seq(0, max.dist, length=bins)))

fit.bio.resid <-variofit(vario.bio.resid, ini.cov.pars=c(0.08, 500/-log(0.05)),
                         cov.model="exponential", minimisation.function="nls", weights="equal")

par(mfrow=c(1,2))
plot(vario.bio, main="Log metric tons of biomass")
lines(fit.bio)
abline(h=fit.bio$nugget, col="blue")##nugget
abline(h=fit.bio$cov.pars[1]+fit.bio$nugget, col="green")##sill
abline(v=-log(0.05)*fit.bio$cov.pars[2], col="red3")##effective range

plot(vario.bio.resid, main="Log metric tons of biomass residuals")
lines(fit.bio.resid)
abline(h=fit.bio.resid$nugget, col="blue")
abline(h=fit.bio.resid$cov.pars[1]+fit.bio.resid$nugget, col="green")
abline(v=-log(0.05)*fit.bio.resid$cov.pars[2], col="red3")


###################################################
### chunk number 7: spExact
###################################################
p <- 6
beta.prior.mean <- as.matrix(rep(0, times=p))
beta.prior.precision <- matrix(0, nrow=p, ncol=p)

phi <- 0.014

alpha <- 0.016/0.08

sigma.sq.prior.shape <- 2.0
sigma.sq.prior.rate <- 0.08

sp.exact <- bayesGeostatExact(log.bio~ELEV+SLOPE+SUM_02_TC1+SUM_02_TC2+SUM_02_TC3, data=BEF.dat, coords=coords, n.samples=1000,
                              beta.prior.mean=beta.prior.mean,
                              beta.prior.precision=beta.prior.precision,
                              phi=phi, alpha=alpha,
                              sigma.sq.prior.shape=sigma.sq.prior.shape,
                              sigma.sq.prior.rate=sigma.sq.prior.rate, sp.effects=FALSE)
summary(sp.exact$p.samples)


###################################################
### chunk number 8: spLM
###################################################
bef.sp <- spLM(log.bio~ELEV+SLOPE+SUM_02_TC1+SUM_02_TC2+SUM_02_TC3, data=BEF.dat, coords=coords,
               starting=list("phi"=3/200,"sigma.sq"=0.08, "tau.sq"=0.02),
               sp.tuning=list("phi"=0.1, "sigma.sq"=0.05, "tau.sq"=0.05),
               priors=list("phi.Unif"=c(3/1500, 3/50), "sigma.sq.IG"=c(2, 0.08),
                 "tau.sq.IG"=c(2, 0.02)),
               cov.model="exponential",
               n.samples=1, verbose=TRUE, n.report=1)


##save(bef.sp, file = "R-data/bef.splm")
load(file = "R-data/bef.splm")


###################################################
### chunk number 9: spLMBeta
###################################################
par(mfrow=c(3,2))
plot(bef.sp$p.samples[,1:6], auto.layout=TRUE, density=FALSE)


###################################################
### chunk number 10: spLMSpatial
###################################################
par(mfrow=c(2,2))
effective.range <- 3/bef.sp$p.samples[,9]
plot(mcmc(cbind(bef.sp$p.samples[,7:9],effective.range)), auto.layout=TRUE, density=FALSE)


###################################################
### chunk number 11: spSummary
###################################################
n.samples <- nrow(bef.sp$p.samples)

burn.in <- 0.25*n.samples
summary(mcmc(bef.sp$p.samples[burn.in:n.samples,]))


###################################################
### chunk number 12: spLMW
###################################################
w.hat.mu <- apply(bef.sp$sp.effects[,burn.in:n.samples],1,mean)
w.hat.sd <- apply(bef.sp$sp.effects[,burn.in:n.samples],1,sd)

par(mfrow=c(1,2))
surf <- mba.surf(cbind(coords, bio.resid), no.X=x.res, no.Y=y.res, extend=FALSE)$xyz.est
z.lim <- range(surf[[3]], na.rm=TRUE)
image.plot(surf, xaxs = "r", yaxs = "r", zlim=z.lim, main="LM residuals")

surf <- mba.surf(cbind(coords, w.hat.mu), no.X=x.res, no.Y=y.res, extend=FALSE)$xyz.est
image.plot(surf, xaxs = "r", yaxs = "r", zlim=z.lim, main="Mean spatial effects")


###################################################
### chunk number 13: spLMSDW
###################################################
par(mfrow=c(1,2))
surf <- mba.surf(cbind(coords, bio.resid), no.X=x.res, no.Y=y.res, extend=FALSE)$xyz.est
image.plot(surf, xaxs = "r", yaxs = "r", main="LM residuals")

surf <- mba.surf(cbind(coords, w.hat.sd), no.X=x.res, no.Y=y.res, extend=FALSE)$xyz.est
image.plot(surf, xaxs = "r", yaxs = "r", main="SD spatial effects")


###################################################
### chunk number 14: spPredict
###################################################
BEF.shp <- readShapePoly("BEF-data/BEF_bound.shp")
BEF.poly <- as.matrix(BEF.shp@polygons[[1]]@Polygons[[1]]@coords)

BEF.grids <- readGDAL("BEF-data/dem_slope_lolosptc_clip_60.img")

pred.covars <- cbind(BEF.grids[["band1"]],BEF.grids[["band2"]],
                     BEF.grids[["band3"]],BEF.grids[["band4"]],BEF.grids[["band5"]])

pred.covars <- cbind(rep(1, nrow(pred.covars)), pred.covars)
pred.coords <- SpatialPoints(BEF.grids)@coords

pred.covars <- pred.covars[pointsInPoly(BEF.poly, pred.coords),]
pred.coords <- pred.coords[pointsInPoly(BEF.poly, pred.coords),]

burn.in <- 0.5*n.samples
bef.bio.pred <- spPredict(bef.sp, pred.coords=pred.coords, pred.covars=pred.covars, start=burn.in, 
                          end=burn.in+1, verbose=FALSE)

##save(bef.bio.pred, file = "R-data/bef.splm.predict")
load(file = "R-data/bef.splm.predict")


###################################################
### chunk number 15: spPredMean
###################################################
bef.bio.pred.mu <- apply(bef.bio.pred$y.pred,1,mean)
bef.bio.pred.sd <- apply(bef.bio.pred$y.pred,1,sd)

surf <- mba.surf(cbind(coords, log.bio), no.X=x.res, no.Y=x.res, extend=TRUE, sp=TRUE)$xyz.est
surf  <- surf [!is.na(overlay(surf, BEF.shp)),]
surf  <- as.image.SpatialGridDataFrame(surf)
z.lim <- range(surf[["z"]], na.rm=TRUE)

pred.grid <- as.data.frame(list(pred.coords, pred.mu=bef.bio.pred.mu, pred.sd=bef.bio.pred.sd))
coordinates(pred.grid) = c("x", "y")
gridded(pred.grid) <- TRUE

pred.mu.image <- as.image.SpatialGridDataFrame(pred.grid["pred.mu"])

par(mfrow=c(1,2))
image.plot(surf, axes=TRUE, zlim=z.lim, col=tim.colors(25), xaxs = "r", yaxs = "r", main="Log metric tons of biomass")
plot(BEF.shp, add=TRUE)

image.plot(pred.mu.image, zlim=z.lim, col=tim.colors(25), xaxs = "r", yaxs = "r", main="Mean predicted log metric tons of biomass")
plot(BEF.shp, add=TRUE)


###################################################
### chunk number 16: spPredSD
###################################################
pred.sd.image <- as.image.SpatialGridDataFrame(pred.grid["pred.sd"])
image.plot(pred.sd.image, axes=TRUE, col=tim.colors(25), xaxs = "r", yaxs = "r")
plot(BEF.shp, add=TRUE)
points(coords, cex=1)


###################################################
### chunk number 17: predExport
###################################################
writeGDAL(pred.grid["pred.mu"], "BEF_Pred_mu_biomass.tif")
writeGDAL(pred.grid["pred.sd"], "BEF_Pred_sd_biomass.tif")


###################################################
### chunk number 18: modelComp
###################################################
set.seed(1)
n.samples <- 1

##Candidate model 1
lm.obj <- lm(log.bio~ELEV+SLOPE+SUM_02_TC1+SUM_02_TC2+SUM_02_TC3, data=BEF.dat)

cand.1 <- bayesLMRef(lm.obj, n.samples)

##save(cand.1, file = "R-data/dic.cand.1")
load(file = "R-data/dic.cand.1")

##Candidate model 2
cand.2 <- spLM(log.bio~1, data=BEF.dat, coords=coords,
               starting=list("phi"=3/200,"sigma.sq"=0.08, "tau.sq"=0.02),
               sp.tuning=list("phi"=0.1, "sigma.sq"=0.05, "tau.sq"=0.05),
               priors=list("phi.Unif"=c(3/1500, 3/50), "sigma.sq.IG"=c(2, 0.08),
                 "tau.sq.IG"=c(2, 0.02)),
               cov.model="exponential",
               n.samples=n.samples, verbose=FALSE)

##save(cand.2, file = "R-data/dic.cand.2")
load(file = "R-data/dic.cand.2")


##Candidate model 3
cand.3 <- spLM(log.bio~ELEV+SLOPE+SUM_02_TC1+SUM_02_TC2+SUM_02_TC3, data=BEF.dat, coords=coords,
               starting=list("phi"=3/200,"sigma.sq"=0.08, "tau.sq"=0.02),
               sp.tuning=list("phi"=0.1, "sigma.sq"=0.05, "tau.sq"=0.05),
               priors=list("phi.Unif"=c(3/1500, 3/50), "sigma.sq.IG"=c(2, 0.08),
                 "tau.sq.IG"=c(2, 0.02)),
               cov.model="exponential",
               n.samples=n.samples, verbose=FALSE)

##save(cand.3, file = "R-data/dic.cand.3")
load(file = "R-data/dic.cand.3")

cand.1.DIC <- spDIC(cand.1, start=2500, thin=10)
cand.2.DIC <- spDIC(cand.2, start=2500, thin=10)
cand.3.DIC <- spDIC(cand.3, start=2500, thin=10)



###################################################
### chunk number 19: tabMargDIC
###################################################
marg <- t(round(cbind(cand.1.DIC, cand.2.DIC$DIC.marg, cand.3.DIC$DIC.marg),1))
rownames(marg) <- c("Non-spatial", "Spatial intercept only", "Spatial with predictors")

xtable(marg, caption = "Candidate model comparison using marginalized DIC", label = "tab:margDIC",
table.placement = "tbp", caption.placement = "top")


###################################################
### chunk number 20: tabUnmargDIC
###################################################
marg <- t(round(cbind(cand.2.DIC$DIC.unmarg, cand.3.DIC$DIC.unmarg),1))
rownames(marg) <- c("Spatial intercept only", "Spatial with predictors")

xtable(marg, caption = "Candidate model comparison using unmarginalized DIC", label = "tab:margDIC",
table.placement = "tbp", caption.placement = "top")