###################################################
### chunk number 1: libs
###################################################
options(width = 65) 
library(spBayes)
library(MBA)
library(geoR)
library(fields)
library(cluster)
library(sp)
library(maptools)
library(rgdal)
library(classInt)
library(gstat)
library(lattice)
library(xtable)


###################################################
### chunk number 2: coords
###################################################
MN.shp <- readShapePoly("MN-data/minnesota_aea.shp")
MN.spp <- read.table("MN-data/MN_spp.txt", header=TRUE)

n.spp <- MN.spp[,"n.spp"]
coords <- MN.spp[,c("x","y")]

plot(coords, pch=19, cex=0.5, xlab="Easting (m)", ylab="Northing (m)")


###################################################
### chunk number 3: surf
###################################################
res <- 100
surf <- mba.surf(cbind(coords, n.spp), no.X=res, no.Y=res, h=7, extend=TRUE, sp=TRUE)$xyz.est
surf  <- surf[!is.na(overlay(surf, MN.shp)),]
surf  <- as.image.SpatialGridDataFrame(surf)

image.plot(surf, xaxs = "r", yaxs = "r", xlab="Easting (m)", ylab="Northing (m)")
plot(MN.shp, add=TRUE)


###################################################
### chunk number 4: knots
###################################################
x.range <- range(coords[,1])
y.range <- range(coords[,2])

knots <- expand.grid(x=seq(x.range[1], x.range[2], length.out=15),
                     y=seq(y.range[1], y.range[2], length.out=15))

coordinates(knots) <- c("x","y")

knots <- knots[!is.na(overlay(knots, MN.shp))]@coords

plot(MN.shp, axes=TRUE, xlab="Easting (m)", ylab="Northing (m)")
points(knots)


###################################################
### chunk number 5: spGLM
###################################################
coords <- as.matrix(coords/1000)
knots <- as.matrix(knots/1000)

max.dist <- max(iDist(coords))
max.dist

beta.starting <- coefficients(glm(n.spp~1, family="poisson"))
beta.tuning <- t(chol(vcov(glm(n.spp~1, family="poisson"))))
            
n.samples <- 1000

spp.spGLM <- spGLM(n.spp~1, family="poisson", coords=coords, knots=knots,
             starting=list("beta"=beta.starting, "phi"=3/100,"sigma.sq"=0.2, "w"=0),
             tuning=list("beta"=beta.tuning, "phi"=0.1, "sigma.sq"=0.005, "w"=0.0001),
             priors=list("beta.Flat", "phi.Unif"=c(3/400, 3/1), "sigma.sq.IG"=c(2, 0.2)),
             cov.model="exponential",
             n.samples=n.samples, verbose=TRUE, n.report=100)

spp.spGLM$p.samples[,"phi"] <- 3/spp.spGLM$p.samples[,"phi"]

n.samples <- nrow(spp.spGLM$p.samples)

burn.in <- 0.75*n.samples
summary(mcmc(spp.spGLM$p.samples[burn.in:n.samples,]))$quantiles[,c(1,3,5)]

beta.hat <- mean(spp.spGLM$p.samples[burn.in:n.samples,1])
w.hat <- rowMeans(spp.spGLM$sp.effects[,burn.in:n.samples])

n.spp.fitted <-exp(beta.hat+w.hat)



###################################################
### chunk number 6: fitted
###################################################
coords <- MN.spp[,c("x","y")]

par(mfrow=c(2,1))
surf <- mba.surf(cbind(coords,n.spp), no.X=res, no.Y=res, extend=TRUE, sp=TRUE)$xyz.est
surf  <- surf[!is.na(overlay(surf, MN.shp)),]
surf  <- as.image.SpatialGridDataFrame(surf)
image.plot(surf, xaxs = "r", yaxs = "r", xlab="Easting (m)", ylab="Northing (m)", main="Interpolated species count")

surf <- mba.surf(cbind(coords,n.spp.fitted), no.X=res, no.Y=res, extend=TRUE, sp=TRUE)$xyz.est
surf  <- surf[!is.na(overlay(surf, MN.shp)),]
surf  <- as.image.SpatialGridDataFrame(surf)
image.plot(surf, xaxs = "r", yaxs = "r", xlab="Easting (m)", ylab="Northing (m)", main="Fitted species counts")