library(spBayes)

set.seed(1234)



##Portions of this example requires MBA package to make surfaces
library(MBA)

##This is for pretty image maps, but no good for B&W
jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", 
"#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))


###########################################
##  Prediction for bayes.geostat.exact
###########################################
data(FBC07.dat)
Y <- FBC07.dat[1:150,"Y.2"]
coords <- as.matrix(FBC07.dat[1:150,c("coord.X", "coord.Y")])

n.samples <-200
n = length(Y)
p = 1

phi <- 0.15
nu <- 0.5

beta.prior.mean <- as.matrix(rep(0, times=p))
beta.prior.precision <- matrix(0, nrow=p, ncol=p)

alpha <- 5/5

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

##############################
##Simple linear model with
##the default exponential
##spatial decay function
##############################
m.1 <- bayes.geostat.exact(Y~1, n.samples=n.samples,
                           beta.prior.mean=beta.prior.mean,
                           beta.prior.precision=beta.prior.precision,
                           coords=coords, phi=phi, alpha=alpha,
                           sigma.sq.prior.shape=sigma.sq.prior.shape,
                           sigma.sq.prior.rate=sigma.sq.prior.rate,
                           sp.effects=TRUE)


##Now prediction

pred.coords <- expand.grid(seq(0,100,length=10),seq(0,100,length=10))
pred.covars <- as.matrix(rep(1,nrow(pred.coords)))

m.1.pred <- sp.predict(m.1, pred.coords=pred.coords,
                       pred.covars=pred.covars, thin=5)

par(mfrow=c(2,2))
obs.surf <-
  mba.surf(cbind(coords, Y), no.X=100, no.Y=100, extend=T)$xyz.est

image(obs.surf, xaxs = "r", yaxs = "r", main="Observed response", col=jet.colors(200))
points(coords, pch=19, cex=1, col="green")
contour(obs.surf, add=T)

w.hat <- rowMeans(m.1$sp.effects)
w.surf <-
  mba.surf(cbind(coords, w.hat), no.X=100, no.Y=100, extend=T)$xyz.est

image(w.surf, xaxs = "r", yaxs = "r", main="Random effects",col=jet.colors(200))
points(coords, pch=19, cex=1, col="green")
contour(w.surf, add=T)

y.hat <- rowMeans(m.1.pred)
y.surf <-
  mba.surf(cbind(pred.coords, y.hat), no.X=100, no.Y=100, extend=T)$xyz.est

image(y.surf, xaxs = "r", yaxs = "r", main="Predicted response", col=jet.colors(200))
points(pred.coords, pch=19, cex=1, col="black")
rect(0, 0, 50, 50, col=NA, border="green")
contour(y.surf, add=T)

y.var <- apply(m.1.pred, 1, var)
y.surf <-
  mba.surf(cbind(pred.coords, y.var), no.X=100, no.Y=100, extend=T)$xyz.est
image(y.surf, xaxs = "r", yaxs = "r", main="Predicted response\nvariance", col=jet.colors(200))
points(coords, pch=19, cex=1, col="green")
points(pred.coords, pch=19, cex=1, col="black")
rect(0, 0, 50, 50, col=NA, border="green")
contour(y.surf, add=T)



###########################################
##       Prediction for sp.lm
###########################################
data(rf.n200.dat)

Y <- rf.n200.dat$Y
coords <- as.matrix(rf.n200.dat[,c("x.coords","y.coords")])
w <- rf.n200.dat$w

pred.coords <- expand.grid(seq(1,10,1), seq(1,10,1))
n.pred <- nrow(pred.coords)

###############################
##Prediction with a sp.lm model
###############################
m.2 <- sp.lm(Y~1, coords=coords,
             starting=list("phi"=0.6,"sigma.sq"=1, "tau.sq"=1),
             sp.tuning=list("phi"=0.01, "sigma.sq"=0.05, "tau.sq"=0.05),
             priors=list("phi.Unif"=c(0.3, 3), "sigma.sq.IG"=c(2, 1),
               "tau.sq.IG"=c(2, 1)),
             cov.model="exponential",
             n.samples=1000, verbose=TRUE, n.report=100, sp.effects=TRUE)

pred <- sp.predict(m.2, pred.coords,
                   pred.covars=as.matrix(rep(1,n.pred)))

par(mfrow=c(1,2))
obs.surf <-
  mba.surf(cbind(coords, Y), no.X=100, no.Y=100, extend=TRUE)$xyz.est
image(obs.surf, xaxs = "r", yaxs = "r", main="Observed response", col=jet.colors(200))
points(coords)
contour(obs.surf, add=T)

y.hat <- rowMeans(pred$y.pred)
y.pred.surf <-
  mba.surf(cbind(pred.coords, y.hat), no.X=100, no.Y=100, extend=TRUE)$xyz.est

image(y.pred.surf, xaxs = "r", yaxs = "r", main="Predicted response", col=jet.colors(200))
points(coords, pch=1, cex=1)
points(pred.coords, pch=19, cex=1)
contour(y.pred.surf, add=T)
legend(1.5,2.5, legend=c("Obs.", "Pred."), pch=c(1,19),
       cex=c(1,1), bg="white")