##Load libraries
library(spBayes)


###############################################
##subset data
###############################################
set.seed(1234)

data(BEF)

n.subset <- 100
BEF.subset <- BEF[sample(1:nrow(BEF),n.subset),]

##############################################
##general ggt.sp setup and model fit
##############################################

##Specify the priors, hyperparameters, and variance parameter starting values.
sigmasq.prior <- prior(dist="IG", shape=2, scale=30)
tausq.prior <- prior(dist="IG", shape=2, scale=30)
##the prior on phi corresponds to a prior of 500-2000 meters
##for the effective range (i.e., -log(0.05)/0.0015, -log(0.05)/0.006, when
##using the exponential covariance function).
phi.prior <- prior(dist="LOGUNIF", a=0.0015, b=0.006)

var.update.control <- list("sigmasq"=list(sample.order=0, starting=30, tuning=0.5, prior=sigmasq.prior),
                           "tausq"=list(sample.order=0, starting=30, tuning=0.5, prior=tausq.prior),
                           "phi"=list(sample.order=0, starting=0.006, tuning=0.8, prior=phi.prior))

##specify the number of samples
run.control <- list("n.samples" = 1500)

##specify some model, assign the prior and starting values for the regressors
model <- BE_BA_AREA~ELEV+SPR_02_TC1+SPR_02_TC2+SPT_02_TC3

##note, precision of 0 is same as flat prior.
beta.prior <- prior(dist="NORMAL", mu=rep(0, 5), precision=diag(0, 5))

beta.control <- list(beta.update.method="gibbs", beta.starting=rep(0, 5), beta.prior=beta.prior)

model.coords <- cbind(BEF.subset$XUTM, BEF.subset$YUTM)

ggt.sp.out <-ggt.sp(formula=model, data=BEF.subset,
                    coords=model.coords,
                    var.update.control=var.update.control,
                    beta.control=beta.control,
                    cov.model="exponential",
                    run.control=run.control, DIC=TRUE, DIC.start=501, verbose=TRUE)

##Print Metropolis acceptance rates
print(ggt.sp.out$accept)

##Print DIC scores using samples starting at DIC.start in the ggt.sp function
print(ggt.sp.out$DIC)


##Get the effective range of spatial dependence.
eff.range <- -log(0.05)/ggt.sp.out$p.samples[,"phi"]

##Plot the chains with the mcmc function in coda.
mcmc.obj <- mcmc(cbind(ggt.sp.out$p.samples, eff.range))
plot(mcmc.obj)

##Summarize the post burn-in posterior distributions
##We use the "start" option in the mcmc function in coda
mcmc.obj.post.burnin <- mcmc(cbind(ggt.sp.out$p.samples, eff.range), start=501)
print(summary(mcmc.obj.post.burnin))

##ALTERNATIVE method of computing DIC -- using the sp.DIC function
DIC.stuff <- sp.DIC(ggt.sp.out, start=501)

##############################################
##recover random spatial effect
##Y = X^tB + w + e, where w is the random
##spatial effect and e is pure error
##############################################

##get the spatial random effect w
sp.effect.out <- sp.effect(ggt.sp.out, start=501, thin=1)
w <- rowMeans(sp.effect.out$pp.samples)

par(mfrow=c(1,2))
int.obs <- interp(BEF.subset$XUTM, BEF.subset$YUTM, BEF.subset$BE_BA_AREA)
image(int.obs, xlab="Longitude (meters)", ylab="Latitude (meters)", main="Observed density\nof American beech (Y)")
contour(int.obs, add=TRUE)

int.w <- interp(BEF.subset$XUTM, BEF.subset$YUTM, w)
image(int.w, xlab="Longitude (meters)", ylab="Latitude (meters)", main="Spatial random effect density (w)")
contour(int.w, add=TRUE)

###############################################
##prediction
###############################################

##resample
n.pred <- 100
BEF.pred <- BEF[sample(1:nrow(BEF),n.pred),]

pred.coords <- cbind(BEF.pred$XUTM, BEF.pred$YUTM)

pred <- sp.predict(ggt.sp.out, pred.joint=TRUE, pred.coords=pred.coords, pred.covars=BEF.pred, start=500, thin=2)

par(mfrow=c(1,2))
int.obs <- interp(BEF.subset$XUTM, BEF.subset$YUTM, BEF.subset$BE_BA_AREA)
image(int.obs, xlab="Longitude (meters)", ylab="Latitude (meters)", main="Observed density\nof American beech")
contour(int.obs, add=TRUE)

int.pred <- interp(BEF.pred$XUTM, BEF.pred$YUTM, rowMeans(pred$pp.samples))
image(int.pred, xlab="Longitude (meters)", ylab="Latitude (meters)", main="Predicted density\nof American beech")
contour(int.pred, add=TRUE)