library(spBayes)

set.seed(1234)

##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"))

###########################################################
##Fitting univariate spatial models with 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 <- 500
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
##############################
set.seed(1)
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)


print(summary(m.1$p.samples))

##Requires MBA package to
##make surfaces
library(MBA)
par(mfrow=c(1,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)
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="Estimated random effects", col=jet.colors(200))
points(coords)
contour(w.surf, add=T)

##############################
##Simple linear model with
##the matern spatial decay
##function. Note, nu=0.5 so
##should produce the same
##estimates as m.1
##############################
set.seed(1)
m.2 <- bayes.geostat.exact(Y~1, n.samples=n.samples,
                           beta.prior.mean=beta.prior.mean,
                           beta.prior.precision=beta.prior.precision,
                           coords=coords, cov.model="matern",
                           phi=phi, nu=nu, alpha=alpha,
                           sigma.sq.prior.shape=sigma.sq.prior.shape,
                           sigma.sq.prior.rate=sigma.sq.prior.rate)

print(summary(m.2$p.samples))


##############################
##Another example but this
##time with covariates
##############################
data(FORMGMT.dat)

n = nrow(FORMGMT.dat)
p = 5 ##an intercept and four covariates

n.samples <- 50

phi <- 0.0012

coords <- cbind(FORMGMT.dat$Longi, FORMGMT.dat$Lat)
coords <- coords*(pi/180)*6378

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

alpha <- 1/1.5

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

m.4 <-
  bayes.geostat.exact(Y~X1+X2+X3+X4, data=FORMGMT.dat, 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)

print(summary(m.4$p.samples))

##Requires MBA package to
##make surfaces
library(MBA)
par(mfrow=c(1,2))
obs.surf <-
  mba.surf(cbind(coords, resid(lm(Y~X1+X2+X3+X4, data=FORMGMT.dat))),
           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)

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

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


#####################################################################
#####################################################################

###########################################################
##Fitting univariate spatial models with 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

##############################
##Estimating all parameters
##############################

m.1 <- 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)

print(summary(m.1$p.samples))
plot(m.1$p.samples, trace=TRUE, density=FALSE)

##Requires MBA package to
##make surfaces
library(MBA)
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)

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

image(w.surf, xaxs = "r", yaxs = "r", main="Estimated random effects", col=jet.colors(200))
points(coords)
points(m.1$knot.coords, pch=19, cex=1)
contour(w.surf, add=T)

##############################
##Fixing spatial and variance
##parameters
##############################

m.2 <- sp.lm(Y~1, coords=coords,
             fixed=list("phi"=0.6,"sigma.sq"=10, "tau.sq"=1),
             cov.model="exponential",
             n.samples=1000, verbose=TRUE, n.report=100, sp.effects=TRUE)

print(summary(m.2$p.samples))
plot(m.2$p.samples, trace=TRUE, density = FALSE)