set.seed(1234)

library(akima)
library(geoR)
library(geoRglm)
library(fields)

FullData <- read.table("http://www.biostat.umn.edu/~sudiptob/DukeSummerCourse/MeshoDataCo1.txt", header=T)
##FullData <- read.table("BerberisMeshoDataBt1.txt", header=T)

X.dummy <- model.matrix(~factor(HabitatClass)+factor(LULCChange)+factor(cat1970), data=FullData)
X <- cbind(X.dummy,FullData$CanopyClosure,FullData$HeavilyManagedPts,FullData$LogEdgeDistance)

##table(FullData$Copresence)
##  0   1 
## 352 251 

##table(FullData$Coabund)
##  0   1   2   3   4   5 
## 352  12  83 100  53   3 

##> table(FullData$HabitatClass) (categorical)
##  1   2   3   4 
## 130 218  13 242 

##> table(FullData$CanopyClosure) (ordinal)
##  0   1   2   3 
## 173 119 105 206 

##> table(FullData$LULCChange) (categorical)
##  1   2   3   4   5 
## 89  60 234 178  42 

##> table(FullData$cat1970) (categorical)
##  1   2   3   4   5   6 
## 222   1  87  93  63 137 

##table(FullData$HeavilyManagedPts)
##  0  1 
## 564  39 

Copresence.glm <- glm(Copresence~factor(HabitatClass)+factor(LULCChange)+factor(cat1970)+FullData$CanopyClosure+HeavilyManagedPts+LogEdgeDistance, family=binomial(link=logit),data=FullData)

Copresence.glm.summary <- summary(Copresence.glm)

sites.original <- cbind(FullData$Xcoord, FullData$Ycoord)
max.dist <- max(dist(sites.original))
x.limits <- c(min(sites.original[,1]),max(sites.original[,1]))
y.limits <- c(min(sites.original[,2]),max(sites.original[,2]))

sites <- sites.original/max.dist

x.grid <- seq(from=x.limits[1],to=x.limits[2],length=50)
y.grid <- seq(from=y.limits[1],to=y.limits[2],length=50)
sites.predict <-  as.matrix(expand.grid(x.grid,y.grid))

obj <- cbind(sites,FullData$Copresence,FullData$HabitatClass,FullData$LULCChange,FullData$cat1970,FullData$CanopyClosure,FullData$HeavilyManagedPts,FullData$LogEdgeDistance)

geodat <- as.geodata(obj,coords.col=1:2,data.col=3,covar.col=4:9)

model.spec <- model.glm.control(trend.d = ~ factor(FullData$HabitatClass)+factor(FullData$LULCChange)+factor(FullData$cat1970)+FullData$CanopyClosure+FullData$HeavilyManagedPts+FullData$LogEdgeDistance, cov.model="exponential")

prior.spec <-
  prior.glm.control(beta.prior="flat",sigmasq.prior="sc.inv.chisq",sigmasq=1.0,df.sigmasq=2,phi.prior="exponential",phi=0.2,phi.discrete=seq(from=0.1,to=0.5,by=0.01))

mcmc.spec <- mcmc.control(S.scale=0.25,phi.scale=0.0015,burn.in=1000,thin=1,n.iter=5000)

output.spec <- output.glm.control(sim.posterior=TRUE,quantile=c(0.50,0.025,0.975))

bayes1 <- binom.krige.bayes(geodat, units.m="default", locations="no", model=model.spec, prior=prior.spec, mcmc.input=mcmc.spec, output=output.spec)

out.full <- bayes1$posterior

start <- 1
end <- mcmc.spec$n.iter

beta.posterior <- out.full$beta$sample[,start:end]

beta.qnt <- matrix(rep(0,times=3*nrow(beta.posterior)), ncol=3, byrow=T)
myqnt <- c(0.50,0.025,0.975)
for (i in 1:nrow(beta.posterior)) {
beta.qnt[i,] <- quantile(beta.posterior[i,], myqnt)
}

phi.posterior <- out.full$phi$sample[start:end]
phi.qnt <- quantile(out.full$phi$sample, c(0.50,0.025,0.975))

sigmasq.posterior <- out.full$sigmasq$sample[start:end]
sigmasq.qnt <- quantile(sigmasq.posterior, c(0.50,0.025,0.975))

S.posterior <- out.full$simulations[,start:end]

S.median <- rep(0,times=nrow(S.posterior))
S.mean <- rep(0,times=nrow(S.posterior))
S.sd <- rep(0,times=nrow(S.posterior))
S.quantiles <- matrix(nrow=nrow(S.posterior),ncol=length(myqnt))

for (i in 1:nrow(S.posterior)) {
  S.median[i] <- quantile(S.posterior[i,], 0.50)
  S.mean[i] <- mean(S.posterior[i,])
  S.sd[i] <- sqrt(var(S.posterior[i,]))
  S.quantiles[i,] <- quantile(S.posterior[i,], myqnt) 
}

write(t(cbind(sites.original,S.median)), "SpatialProbMedians.txt", ncolumn=3)
write(t(cbind(sites.original,S.mean)), "SpatialProbMeans.txt", ncolumn=3)
write(t(cbind(sites.original,S.mean,S.sd,S.quantiles)), "SpatialProbSummaries.txt", ncolumn=7)

##pretty, but no good for B&W
jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", 
"#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))

##postscript("Pics/TrendSurfaceProbScale.ps",horizontal=F)
S.interp <- interp.new(sites.original[,1], sites.original[,2], S.mean)
image(S.interp, xlab="X-coord", ylab="Y-coord", col=jet.colors(200))
contour(S.interp, add=T)
##dev.off()

S.gp.posterior <- log(S.posterior/(1-S.posterior))

S.median <- rep(0,times=nrow(S.gp.posterior))
S.mean <- rep(0,times=nrow(S.gp.posterior))

for (i in 1:nrow(S.gp.posterior)) {
  S.median[i] <- quantile(S.gp.posterior[i,], 0.50)
  S.mean[i] <- mean(S.gp.posterior[i,])
}

##postscript("Pics/TrendSurfaceGPScale.ps",horizontal=F)
S.interp <- interp.new(sites.original[,1], sites.original[,2], S.mean)
image(S.interp, xlab="X-coord", ylab="Y-coord", col=jet.colors(200))
contour(S.interp, add=T)
##dev.off()

S.resid.posterior <- S.gp.posterior - X%*%beta.posterior
write(t(S.resid.posterior), "SpatialEffects.txt", ncolumn=ncol(S.resid.posterior))

##S.predictive <- bayes1$predictive$simulations
##S.predictive.median <- S.predictive$median
##S.gp.predictive <- log(S.predictive.median/(1-S.predictive.median))

S.median <- rep(0,times=nrow(S.resid.posterior))
S.mean <- rep(0,times=nrow(S.resid.posterior))

for (i in 1:nrow(S.resid.posterior)) {
  S.median[i] <- quantile(S.resid.posterior[i,], 0.50)
  S.mean[i] <- mean(S.resid.posterior[i,])
}



##postscript("Pics/ResidualSurface.ps",horizontal=F)
S.interp <- interp.old(sites.original[,1], sites.original[,2], S.mean, extrap=T)
image(S.interp, xlab="X-coord", ylab="Y-coord", xlim=x.limits, ylim=y.limits, col=jet.colors(200))
contour(S.interp, add=T)
points(sites.original[,1], sites.original[,2], pch=3)
##dev.off()