library(spBayes)

set.seed(1234)

data(BEF)

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

##############################################
##General ggt 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))

##Get the distance matrix.
D <- as.matrix(dist(cbind(BEF.subset$XUTM, BEF.subset$YUTM)))

##Define the variance function (e.g., exponential covariance)
var.function <- function(){
  sigmasq*exp(-phi*D)+tausq*diag(n.subset)
}

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

##Specify the 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 alternatively if you want a flat prior
##do not include the beta.prior tag in the beta.control list.
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)

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

##Print the acceptance ratio
print(ggt.out$accept)

##Print the DIC score (marginalized focus) with a burn-in of 500 given by DIC.start=501 in ggt.out.
print(ggt.out$DIC)


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

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

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