spWombling: An R-package to perform Bayesian curvilinear curvature Wombling on spatial data

Reference to the paper titled, "Bayesian Modeling with Spatial Curvature Processes". Journal of the American Statistical Association (2023) (https://doi.org/10.1080/01621459.2023.2177166).

Illustration of curvilinear Bayesian Wombling on spatial data.

Code for performing curvilinear Bayesian wombling on data

Contents

1. Load the data and separate (a) co-ordinates (b) response (c) covariates
2. Fit a spatial Bayesian hierarchical model to the data
3. Perform gradient, curvature and posterior surface analysis
4. Locate or Annotate curves of interest
5. Perform rectilinear wombling at line segments

Demonstration with synthetic data.

Load the data and separate (a) co-ordinates (b) response (c) covariates

``` if(!require(devtools)) install.packages("devtools") devtools::install_github('arh926/spWombling') require(spWombling)

loading additional packages

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

N <- 100 tau <- 1

synthetic location

coords <- matrix(runif(2*N),nr=N,nc=2)

synthetic pattern

y <- rnorm(N,10(sin(3picoords[,1])+cos(3pi*coords[,2])),tau)

y <- rnorm(N,10(sin(3picoords[,1])cos(3picoords[,2])),tau)

y=y # response X=matrix(1,nr=N) # intercept ```

Fit a spatial Bayesian hierarchical model to the data

``` niter <- 1e4 nburn <- niter/2 report <- 1e2 mcsp <- hlmBayessp(coords = coords, y = y, X = X, niter = niter, nburn = nburn, report = report, cov.type = "matern2", steps_init = 1, verbose = T)

inference from model

modelsummary <- hlmsummary(chain = mcsp, nburn = nburn, niter = niter, thin = 1) coef <- modelsummary$summary.pars; round(coef,4)

z <- model_summary$summary.latent ```

Perform gradient, curvature and posterior surface analysis

``` grid.points <- as.matrix(expand.grid(seq(0,1,by=0.05)[-c(1,21)], seq(0,1,by=0.05)[-c(1,21)]),nc=2) samples <- (nburn+1):niter gradientest <- spatialgradient(coords=coords, grid.points = grid.points, samples = samples, chain = mc_sp, cov.type = "matern2", niter = niter, nburn = nburn, nbatch = 100, return.mcmc = T)

Gradient and Curvature Assessment

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(gradientest$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(gradientest$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

detest <- eigenest1 <- eigenest2 <- laplaceest <- divest <- matrix(NA,nr=niter,nc=nrow(grid.points)) for(i in 1:nburn){ for(j in 1:nrow(grid.points)){ hess.mat <- matrix(c(gradientest$grad.s11.mcmc[i,j], gradientest$grad.s12.mcmc[i,j], gradientest$grad.s12.mcmc[i,j], gradientest$grad.s22.mcmc[i,j]),nr=2,nc=2,byrow = T) detest[i,j] <- det(hess.mat) eigen.hess.mat <- eigen(hess.mat)$values eigenest1[i,j] <- eigen.hess.mat[1]; eigenest2[i,j] <- eigen.hess.mat[2] laplaceest[i,j] <- sum(diag(hess.mat)) divest[i,j] <- gradientest$grad.s1.mcmc[i,j]+gradientest$grad.s2.mcmc[i,j] } } slist <- split(1:nburn, ceiling(seq_along(1:nburn)/(length(1:nburn)/100)))

Determinant Surface

detestval <- cbind.data.frame(apply(do.call(rbind,lapply(slist, function(x) apply(detest[x,],2,median))),2,median), HPDinterval(as.mcmc(do.call(rbind,lapply(slist, function(x) apply(detest[x,],2,median)))))) detestval$sig <- apply(detestval,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 Value Surface

eigenestval1 <- cbind.data.frame(apply(do.call(rbind,lapply(slist, function(x) apply(eigenest1[x,],2,median))),2,median), HPDinterval(as.mcmc(do.call(rbind,lapply(slist, function(x) apply(eigenest1[x,],2,median)))))) eigenestval1$sig <- apply(eigenestval1,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 Operator Surface

laplaceestval <- cbind.data.frame(apply(do.call(rbind,lapply(slist, function(x) apply(laplaceest[x,],2,median))),2,median), HPDinterval(as.mcmc(do.call(rbind,lapply(slist, function(x) apply(laplaceest[x,],2,median)))))) laplaceestval$sig <- apply(laplaceestval,1,function(x){ if(x[2]>0 & x[3]>0) return (1) if(x[2]<0 & x[3]<0) return (-1) else return(0) })

Divergence Operator

divestval <- cbind.data.frame(apply(do.call(rbind,lapply(slist, function(x) apply(divest[x,],2,median))),2,median), HPDinterval(as.mcmc(do.call(rbind,lapply(slist, function(x) apply(divest[x,],2,median)))))) divestval$sig <- apply(divestval,1,function(x){ if(x[2]>0 & x[3]>0) return (1) if(x[2]<0 & x[3]<0) return (-1) else return(0) }) ```

True and Estimated Gradients

True and Estimated Curvature

Locate or Annotate curves of interest

```

locate curves:: level sets/contours

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

rastr.obj <- sp_plot(11,"Spectral",cbind(coords,y), contour.plot = T,raster.surf = T, legend = F) x <- raster::rasterToContour(rastr.obj,nlevels=20)

Test for wombling

x.levels <- as.numeric(as.character(x$level)) level.choice <- c(x.levels[2],x.levels[length(x.levels)-1]) subset.points.1 <- subset(x,level==level.choice[1]) subset.points.2 <- subset(x,level==level.choice[2]) x <- raster::rasterToContour(rastr.obj,nlevels=10) subset.points.3 <- subset(x,level==15)

curve-1

lines(subset.points.1@lines[[1]]@Lines[[1]]@coords,lwd=2.5)

curve-2

lines(subset.points.2@lines[[1]]@Lines[[3]]@coords,lwd=2.5)

curve 3

lines(subset.points.3@lines[[1]]@Lines[[3]]@coords,lwd=2.5)

locate curves:: annotate curves

pts.hull <- locator(20) subset.points.4 <- bezierCurve(pts.hull$x,pts.hull$y,100)

curve 4

lines(subset.points.4, type="l",lwd=2.5) points(pts.hull, cex=0.5, col="white") points(pts.hull, pch="+", cex=0.5)

legendimage <- as.raster(matrix(colorRampPalette(RColorBrewer::brewer.pal(11,"Spectral"))(100), ncol=1)) plot(c(0,3),c(0,1),type = 'n', axes = F,xlab = '', ylab = '', main = '') text(x=2, y = seq(0.01,0.99,l=6), labels = sprintf("%.2f",round(seq(min(y),max(y),l=6),2))) rasterImage(legendimage, 0, 0, 1,1) ```

Perform rectilinear wombling at line segments

curve <- subset.points.1@lines[[1]]@Lines[[1]]@coords womb.measure <- bayes_cwomb(coords = coords, chain = mc_sp, cov.type = "matern2", niter = niter, nburn = nburn, nbatch = 100, curve = curve, type = "rectilinear") str(womb.measure) womb.measure$womb.measure.inf

Curvilinear Gradient

Curvilinear Curvature

Authors

| Name | Email | | |:------ |:----------- | :----------- | | Aritra Halder (maintainer)| aritra.halder@drexel.edu | Assistant Professor, Department of Biostatistics, Drexel University|
| Sudipto Banerjee | sudipto@ucla.edu | Professor and Chair, Department of Biostatistics, UCLA | | Dipak K. Dey | dipak.dey@uconn.edu | Board of Trustees Distinguished Professor, Department of Statistics, UCONN | <!--- --->