#################################################################################
# Bayesian model selection and statistical modeling 
# Chapman & Hall/CRC Taylor and Francis Group 
#
# Chapter5: Section5.8.3
#
# Bayesian spatial modeling
#
# Author : Tomohiro Ando
#
################################################################################

library(spBayes)

# Bartlett experimental forest inventory data

data(BEF.dat)
Data <- BEF.dat[BEF.dat$ALLBIO02_KGH>0,]

bio     <- Data$ALLBIO02_KGH*0.001;
log.bio <- log(bio)
coords  <- as.matrix(Data[,c("XUTM","YUTM")])


#MCMC estimation

bef.sp <- spLM(
#log.bio~ELEV+SLOPE+SUM_02_TC1+SUM_02_TC2+SUM_02_TC3, #Model 1
log.bio~ELEV+SLOPE,  #Model 2
data=BEF.dat, coords=coords,
starting=list("phi"=0.01,"sigma.sq"=0.05, "tau.sq"=0.03),
sp.tuning=list("phi"=0.01, "sigma.sq"=0.05, "tau.sq"=0.05),
priors=list("phi.Unif"=c(0.0001,0.05), "sigma.sq.IG"=c(10^-5,10^-5),
"tau.sq.IG"=c(10^-5,10^-5)),
cov.model="exponential",
n.samples=6000, 
sub.samples=c(1000, 5000, 10),
verbose=TRUE, n.report=1)


# PL score calculation

dic.m2 <- spDIC(bef.sp)

Likelihood <- dic.m2$DIC.marg

PL2 <- -1/2*Likelihood-ncol(bef.sp$p.samples)/2