In arh926/spWombling: Bayesian Wombling for Spatial Processes

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

require(spWombling)
require(sp)
require(coda)
require(Matrix)

data(boston)

boston.bdry <- chull(boston.c$LON,boston.c$LAT)
boston.c[boston.bdry,c("LON","LAT")]
boston.shp <- raster::spPolygons(as.matrix(boston.c[boston.bdry,c("LON","LAT")],nc=2), crs=sp::CRS("+proj=longlat +datum=WGS84 +no_defs"))

mat <- matrix(c(1,2,3), nr=1,nc=3, byrow=T)
layout(mat,
       widths = c(5,5,2.5),
       heights = c(3,3))
hist(boston.c$CMEDV, xlab="Median House Prices (in USD 1000)", main="", col="lightblue",breaks=50)
sp_plot(11,"Spectral",data_frame = cbind(boston.c$LON,
                                         boston.c$LAT,
                                         boston.c$CMEDV),
        xlab=latex2exp::TeX("Longitude$\\degree$"),
        ylab = latex2exp::TeX("Latitude$\\degree$"),
        shape=boston.shp,
        contour.plot = T,points.plot = T)

````r coords <- cbind(boston.c$LON, boston.c$LAT)*100

cnames <- c("(Intercept)","CRIM","ZN","INDUS","CHAS","NOX.2","RM.2","AGE","lDIS","lRAD","TAX","PTRATIO","B","lLSTAT")

N <- nrow(coords) y <- boston.c$CMEDV X <- matrix(1,nrow=N) niter <- 2e3 nburn <- niter/2 report <- 1e2 mc_sp <- hlmBayes_sp(coords = coords, y = y, X = X, niter = niter, nburn = nburn, report = report, cov.type = "matern2", steps_init = 1, verbose = T)

```r
model_summary <- hlm_summary(chain = mc_sp,nburn = nburn,niter = niter,thin = 1)
coef <- model_summary$summary.pars; round(coef,4)

z <- model_summary$summary.latent
z$sig <- apply(z,1,function(x){
  if(x[2]>0 & x[3]>0) return (1)
  if(x[2]<0 & x[3]<0) return (-1)
  else return(0)
})

1-sum(apply(cbind(y,z[,"lower.hpd"]+coef["post_beta","lower.hpd"], z[,"upper.hpd"]+coef["post_beta","upper.hpd"]),
            1,
            function(x){
              if(x[1]>=x[2] & x[1]<=x[3]) return(0)
              else return(1)
            }))/N
y.hat <- z[,"median"]+coef["post_beta","median"]
y.upper.hpd <- z[order(y.hat),"upper.hpd"]+coef["post_beta","upper.hpd"]
y.lower.hpd <- z[order(y.hat),"lower.hpd"]+coef["post_beta","lower.hpd"]


# MSE
round(sqrt(mean((y-y.hat)^2)),2)
# Check with LM
gp0 <- lm(CMEDV ~ 1, data = boston.c)
summary(gp0)
coef;coefficients(gp0)

mat <- matrix(c(1,2,3), nr=1,nc=3, byrow=T)
layout(mat,
       widths = c(5,5,2.5),
       heights = c(3,3))
sp_plot(11,"Spectral",cbind(coords/100,y.hat),contour.plot = T,shape=boston.shp, points.plot = T,sig = z$sig, legend = F) 
sp_plot(11,"Spectral",cbind(coords/100,y), contour.plot = T,shape=boston.shp, points.plot = T)

# RMSE
par(mfcol=c(1,1))
sqrt(mean((y-y.hat)^2))
plot(y,y.hat)
plot(y-y.hat)

grid.points <- as.matrix(expand.grid(seq(min(boston.c$LON),max(boston.c$LON),by=0.03),
                                     seq(min(boston.c$LAT),max(boston.c$LAT),by=0.03)),nc=2)
tmp <- over(SpatialPoints(grid.points,proj4string = CRS(proj4string(boston.shp))),boston.shp)
grid.points <- grid.points[!is.na(tmp),]
sp_plot(11,"Spectral",cbind(coords/100,y),shape=boston.shp,contour.plot=T,legend=F)
points(grid.points, cex=0.2)
dim(grid.points)
grid.points <- grid.points*100
samples <- (nburn+1):niter
gradient_est <- spatial_gradient(coords=coords,
                                 grid.points = grid.points,
                                 samples = samples,
                                 chain = mc_sp,
                                 cov.type = "matern2",
                                 niter = 2e3,
                                 nburn = 1e3,
                                 nbatch = 100,
                                 return.mcmc = T)

grad.s1.hpd <- data.frame(gradient_est$grad1.est)
grad.s1.hpd$signif <- apply(grad.s1.hpd,1,function(x){
  if(x[2]>0 & x[3]>0) return (1)
  if(x[2]<0 & x[3]<0) return (-1)
  else return(0)
})  

grad.s2.hpd <- data.frame(gradient_est$grad2.est)
grad.s2.hpd$signif <- apply(grad.s2.hpd,1,function(x){
  if(x[2]>0 & x[3]>0) return (1)
  if(x[2]<0 & x[3]<0) return (-1)
  else return(0)
})  

grad.s11.hpd <- data.frame(gradient_est$grad11.est)
grad.s11.hpd$signif <- apply(grad.s11.hpd,1,function(x){
  if(x[2]>0 & x[3]>0) return (1)
  if(x[2]<0 & x[3]<0) return (-1)
  else return(0)
})  
grad.s12.hpd <- data.frame(gradient_est$grad12.est)
grad.s12.hpd$signif <- apply(grad.s12.hpd,1,function(x){
  if(x[2]>0 & x[3]>0) return (1)
  if(x[2]<0 & x[3]<0) return (-1)
  else return(0)
})  
grad.s22.hpd <- data.frame(gradient_est$grad22.est)
grad.s22.hpd$signif <- apply(grad.s22.hpd,1,function(x){
  if(x[2]>0 & x[3]>0) return (1)
  if(x[2]<0 & x[3]<0) return (-1)
  else return(0)
})  

# Posterior Surface Analysis
det_est <- eigen_est_1 <- eigen_est_2 <- laplace_est <- div_est <- matrix(NA,nr=niter,nc=nrow(grid.points))
for(i in 1:1000){
  for(j in 1:nrow(grid.points)){
    hess.mat <- matrix(c(gradient_est$grad.s11.mcmc[i,j],
                         gradient_est$grad.s12.mcmc[i,j],
                         gradient_est$grad.s12.mcmc[i,j],
                         gradient_est$grad.s22.mcmc[i,j]),nr=2,nc=2,byrow = T)
    det_est[i,j] <- det(hess.mat)
    eigen.hess.mat <- eigen(hess.mat)$values
    eigen_est_1[i,j] <- eigen.hess.mat[1]; eigen_est_2[i,j] <- eigen.hess.mat[2]
    laplace_est[i,j] <- sum(diag(hess.mat))
    div_est[i,j] <- gradient_est$grad.s1.mcmc[i,j]+gradient_est$grad.s2.mcmc[i,j]
  }
}
slist <- split(1:1000, ceiling(seq_along(1:1000)/(length(1:1000)/100)))
det_est_val <- cbind.data.frame(apply(do.call(rbind,lapply(slist, function(x) apply(det_est[x,],2,median))),2,median),
                                HPDinterval(as.mcmc(do.call(rbind,lapply(slist, function(x) apply(det_est[x,],2,median))))))
det_est_val$sig <- apply(det_est_val,1,function(x){
  if(x[2]>0 & x[3]>0) return (1)
  if(x[2]<0 & x[3]<0) return (-1)
  else return(0)
})
eigen_est_val1 <- cbind.data.frame(apply(do.call(rbind,lapply(slist, function(x) apply(eigen_est_1[x,],2,median))),2,median),
                                   HPDinterval(as.mcmc(do.call(rbind,lapply(slist, function(x) apply(eigen_est_1[x,],2,median))))))
eigen_est_val1$sig <- apply(eigen_est_val1,1,function(x){
  if(x[2]>0 & x[3]>0) return (1)
  if(x[2]<0 & x[3]<0) return (-1)
  else return(0)
})
eigen_est_val2 <- cbind.data.frame(apply(do.call(rbind,lapply(slist, function(x) apply(eigen_est_2[x,],2,median))),2,median),
                                   HPDinterval(as.mcmc(do.call(rbind,lapply(slist, function(x) apply(eigen_est_2[x,],2,median))))))
eigen_est_val2$sig <- apply(eigen_est_val2,1,function(x){
  if(x[2]>0 & x[3]>0) return (1)
  if(x[2]<0 & x[3]<0) return (-1)
  else return(0)
})
laplace_est_val <- cbind.data.frame(apply(do.call(rbind,lapply(slist, function(x) apply(laplace_est[x,],2,median))),2,median),
                                    HPDinterval(as.mcmc(do.call(rbind,lapply(slist, function(x) apply(laplace_est[x,],2,median))))))
laplace_est_val$sig <- apply(laplace_est_val,1,function(x){
  if(x[2]>0 & x[3]>0) return (1)
  if(x[2]<0 & x[3]<0) return (-1)
  else return(0)
})
div_est_val <- cbind.data.frame(apply(do.call(rbind,lapply(slist, function(x) apply(div_est[x,],2,median))),2,median),
                                HPDinterval(as.mcmc(do.call(rbind,lapply(slist, function(x) apply(div_est[x,],2,median))))))
div_est_val$sig <- apply(div_est_val,1,function(x){
  if(x[2]>0 & x[3]>0) return (1)
  if(x[2]<0 & x[3]<0) return (-1)
  else return(0)
})

mat <- matrix(c(1:6), nr=1,nc=6, byrow=T)
layout(mat,
       widths = c(rep(c(5,3),3)),
       heights = c(rep(3,6)))
rster.obj <- sp_plot(11,"Spectral",cbind(coords/100,z[,1]),
                     shape=boston.shp,
                     contour.plot = T,
                     xlab=latex2exp::TeX("Longitude$\\degree$"),
                     ylab = latex2exp::TeX("Latitude$\\degree$"),
                     points.plot=T,
                     raster.surf = T,
                     sig = z$sig)

sp_plot(11,"Spectral",cbind(grid.points/100,grad.s1.hpd[,1]),
        shape = boston.shp,contour.plot = T,
        xlab=latex2exp::TeX("Longitude$\\degree$"),
        ylab = latex2exp::TeX("Latitude$\\degree$"),
        points.plot = T, sig=grad.s1.hpd$signif)


sp_plot(11,"Spectral",cbind(grid.points/100,grad.s2.hpd[,1]), shape=boston.shp,contour.plot = T,
        xlab=latex2exp::TeX("Longitude$\\degree$"),
        ylab = latex2exp::TeX("Latitude$\\degree$"),
        points.plot = T, sig=grad.s2.hpd$signif)

sp_plot(11,"Spectral",cbind(grid.points/100,grad.s11.hpd[,1]),contour.plot = T, shape=boston.shp,
        xlab=latex2exp::TeX("Longitude$\\degree$"),
        ylab = latex2exp::TeX("Latitude$\\degree$"),
        points.plot = T, sig=grad.s11.hpd$signif)

sp_plot(11,"Spectral",cbind(grid.points/100,grad.s12.hpd[,1]),contour.plot = T, shape=boston.shp,
        xlab=latex2exp::TeX("Longitude$\\degree$"),
        ylab = latex2exp::TeX("Latitude$\\degree$"),
        points.plot = T, sig=grad.s12.hpd$signif)


sp_plot(11,"Spectral",cbind(grid.points/100,grad.s22.hpd[,1]),contour.plot = T, shape=boston.shp,
        xlab=latex2exp::TeX("Longitude$\\degree$"),
        ylab = latex2exp::TeX("Latitude$\\degree$"),
        points.plot = T, sig=grad.s22.hpd$signif)

mat <- matrix(c(1:4), nr=1,nc=4, byrow=T)
layout(mat,
       widths = c(rep(c(5,3),2)),
       heights = c(rep(3,4)))

sp_plot(11,"Spectral",cbind(grid.points/100,eigen_est_val1[,1]),
        contour.plot = T, shape=boston.shp,
        xlab=latex2exp::TeX("(a) $\\widehat{\\lambda_1}$"),
        points.plot=T,sig=eigen_est_val1$sig)


sp_plot(11,"Spectral",cbind(grid.points/100,eigen_est_val2[,1]),
        contour.plot = T, shape=boston.shp,
        xlab=latex2exp::TeX("(b) $\\widehat{\\lambda_2}$"),
        points.plot=T,sig=eigen_est_val2$sig)

sp_plot(11,"Spectral",cbind(grid.points/100,det_est_val[,1]),
        contour.plot = T, shape=boston.shp,
        xlab=latex2exp::TeX("(c) $\\widehat{|\\nabla^2Y|}$"),
        points.plot=T,sig=det_est_val$sig)

sp_plot(11,"Spectral",cbind(grid.points/100,div_est_val[,1]),
        contour.plot = T, shape=boston.shp,
        xlab=latex2exp::TeX("Longitude$\\degree$"),
        ylab = latex2exp::TeX("Latitude$\\degree$"),
        points.plot=T,sig=div_est_val$sig,
        grid=F)

sp_plot(11,"Spectral",cbind(grid.points/100,laplace_est_val[,1]),
        contour.plot = T, shape=boston.shp,
        xlab=latex2exp::TeX("Longitude$\\degree$"),
        ylab = latex2exp::TeX("Latitude$\\degree$"),
        points.plot=T,sig=laplace_est_val$sig,
        grid=F)

rster.obj <- sp_plot(11,"Spectral",cbind(coords/100,z[,1]),
                     shape=boston.shp,
                     contour.plot = T,
                     xlab=latex2exp::TeX("Longitude$\\degree$"),
                     ylab = latex2exp::TeX("Latitude$\\degree$"),
                     points.plot=T,
                     raster.surf = T,
                     sig = z$sig,
                     legend = F)
# Selecting a Curve
x <- raster::rasterToContour(rster.obj,nlevels=20)
# Test for wombling
x.levels <- as.numeric(as.character(x$level))
level.choice <- c(x.levels[2],x.levels[18]) #c(-10,10)
#c(x.levels[2],x.levels[length(x.levels)-1])
subset.points.1 <- subset(x,level==level.choice[1])

sp_plot(11,"Spectral",cbind(coords/100,z[,1]),
        shape=boston.shp,
        contour.plot = T,
        xlab=latex2exp::TeX("Longitude$\\degree$"),
        ylab = latex2exp::TeX("Latitude$\\degree$"),
        legend = F)
lines(subset.points.1@lines[[1]]@Lines[[1]]@coords,lwd=1.5, lty="dashed")
curve <- subset.points.1@lines[[1]]@Lines[[1]]@coords

womb.measure <- bayes_cwomb(coords = coords,
                            chain = mc_sp,
                            cov.type = "matern2",
                            niter = 1000,
                            nburn = 900,
                            nbatch = 10,
                            curve = curve,
                            type = "rectilinear")
str(womb.measure)
womb.measure$womb.measure.inf