[This article was first published on r-spatial, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.* Summary * Installation + From github + From CRAN * Getting started * Adjacency matrix * Non-spatial regression + MCMC output + Methods * Spatial regression + A filtering approach + A bivariate model + Predicted values * Future work and support * Appendix * References
[view rawRmd]
SummaryThis post introduces thegeostan R package for spatialanalysis. The package is mainly oriented towards areal data, althoughsome models may also be used for other spatial data types. The packageimplements the spatial error/simultaneous spatial autoregressive (SAR)model, conditional autoregressive (CAR) model, and eigenvector spatialfilter (ESF) models for spatial regression. A version of ESF modellingalso appears in the ecology literature as principle coordinate analysisof neighbor matrices (PCNM) (Griffith and Peres-Neto 2006).
geostan also supports the application of the above regression methodsto hierarchical models for count data, as is common in analyses ofdisease incidence or mortality in small areas (‘disease mapping’).Additional features of the software include models forsampling/measurement error in covariates and methods for handlingcensored count data, such as when mortality or disease counts have beencensored for privacy. The models were built using theStan modeling language, so all inference iscompleted using Markov chain Monte Carlo (MCMC) sampling (StanDevelopment Team 2023; Gabry et al. 2024). The spatial autoregressivemodels use custom-built Stan functions that speed up MCMC samplingconsiderably (Donegan 2021).
This post will walk through an example analysis using international dataon life expectancy and per capita GDP. Package vignettes can be foundwith the onlinedocumentation,including an introduction to spatial weights matrices, exploratoryspatial data analysis, spatial measurement error models, rasterregression, and using geostan to build custom spatial models with Stan.A paper in the Journal of Open SourceSoftware reports these and otherfeatures and provides the recommended citation when using geostan(Donegan 2022).
Installationgeostan is currently on CRAN, although that may not always be thecase. You can also install directly from geostan’s githubrepository.
From githubYou can install from the package github repository:
if (!require('devtools')) install.packages('devtools')devtools::install\_github("connordonegan/geostan")
If you are using Windows and installing using install_github, you mayneed to install Rtoolsfirst. Rtools is not needed when installing from CRAN. You may alsocontact the author by e-mail for a pre-compiled version that you can usewithout Rtools.
If you are using Mac and installing with install_github then you mayneed to install Xcode Command Line Tools first.
From CRANUsing your R console, you can install from CRAN as follows:
install.packages("geostan")
Getting startedTo begin, load the geostan and sf packages into your R environment,as well as the world data:
library(geostan)library(sf)data(world, package = "spData")
The world data contains life expectancy and gross domestic product(GDP) per capita (presumably measured in current $US) for 161 countriesas of 2014, gathered from the World Bank. The rest of this post is goingto be structured around a bivariate analysis of these variables.
We are going to apply the Robinson map projection for the countries:
world <- st\_transform(world, crs = 'ESRI:54030')
At least a couple of the missing values can be filled in using WorldBank data, so we will do that:
```
``` And we will also remove Antarctica:
world <- subset(world, name\_long != "Antarctica")
Mapping the variables shows the familiar geography of high-, middle-,and low-income countries and a similar geography of longevity:
```
``` Choropleth maps of GDP per capita and life expectancy.
par(ogpar)
The map_pars function breaks the variables into quantiles and returnsbreaks, colors, labels for the maps; it can be found at the end of thepost. There will be no discussion of substantive (non-statistical)issues here, for which one can consult any number of texts on globalpower and inequality (e.g., Paul Farmer’s Infections andInequalities).
By conventional methods, the correlation coefficient for life expectancyand log GDP per capita is 0.81:
log\_x <- log10( world$gdpPercap )y <- world$lifeExpcor.test(log\_x, y)## ## Pearson's product-moment correlation## ## data: log\_x and y## t = 17.394, df = 160, p-value < 2.2e-16## alternative hypothesis: true correlation is not equal to 0## 95 percent confidence interval:## 0.7478143 0.8561778## sample estimates:## cor ## 0.8087528
The conventional assessment is based on the proposition that we have 161independent observations. The visible geography of the variables, andany level of social awareness, indicates that these are not independentobservations. Rather, there are various functional regions of countriesthat share basic political-economic conditions. A lot, but not all, ofthe variation can be described as variation across continents andregions. We will want to account for this dependence using a spatialmodel (for background see Chun and Griffith 2012; Donegan 2024). Thefirst step will be to construct a spatial weights matrix.
Adjacency matrixThis section will illustrate use of two geostan functions for creatingan revising a spatial weights matrix: shape2mat and edges. Theshape2mat function may be helpful for some users but one can always dothis using spdep or other methods, especially if shape2mat does notprovide the exact method you’re looking for.
We are going to start by removing the 15 countries that are missingvalues:
```
``
Now we can apply theshape2matfunction to obtain an adjacencymatrixthat encodes spatial adjacency relations for countries into a binaryN-by-N matrix. The function usesspdepto find adjacency relations andreturns results as a sparse matrix (using theMatrix` package):
A <- shape2mat(world, "B", method = "rook")## Contiguity condition: rook## Number of neighbors per unit, summary:## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0.000 2.000 3.000 3.605 5.000 13.000## ## Spatial weights, summary:## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 1 1 1 1 1 1
Visualizing the connections in the matrix is important for uncoveringunexpected results. geostan’s edges function converts the matrixinto a list of nodes and edges that we can plot. For this we need tosupply the function with the adjacency matrix, A, and the associatedspatial object, world:
```
```
par(ogpar)
This reveals quite a few unexpected results. French Guiana is stored inthe world data as part of France (a multi-part polygon); this iscorrect of course but it leads to Brazil and Suriname being listed asneighbors of France, which is not sensible. Besides removing thoseconnections, there are a number of island nations that we might want toconnect to nearby places.
To connect Mozambique to Madagascar, we just replace the zeroes withones in the slots that correspond to those countries. First we grabtheir index positions in the matrix:
moz\_idx <- grep("Mozambique", world$name\_long)mad\_idx <- grep("Madagascar", world$name\_long)
And then we assign the correct slots in the matrix a value of 1 (orTRUE), remembering that the adjacency matrix is symmetric:
A[moz\_idx, mad\_idx] <- A[mad\_idx, moz\_idx] <- TRUE
This can become tedious but it is important. Before moving on, we willmake a series of adjustments. This will be made a bit easier with thisconvenience function:
connect <- function(country\_a, country\_b, names\_vec = world$name\_long, matrix = A, add = TRUE) { stopifnot( country\_a %in% names\_vec ) stopifnot( country\_b %in% names\_vec ) a\_idx <- which(names\_vec == country\_a) b\_idx <- which( names\_vec == country\_b) matrix[a\_idx, b\_idx] <- matrix[b\_idx, a\_idx] <- add return( matrix )}
The following are at least reasonable changes to make; they also ensurethat every country has at least one neighbor:
A <- connect("Mozambique", "Madagascar")A <- connect("Australia", "New Zealand")A <- connect("Philippines", "Malaysia")A <- connect("Japan", "Republic of Korea")A <- connect("Fiji", "Vanuatu")A <- connect("Solomon Islands", "Vanuatu")A <- connect("Solomon Islands", "Papua New Guinea")A <- connect("Australia", "Papua New Guinea")A <- connect("Haiti", "Jamaica")A <- connect("Bahamas", "United States")A <- connect("Dominican Republic", "Puerto Rico")A <- connect("Trinidad and Tobago", "Venezuela")A <- connect("Sri Lanka", "India")A <- connect("Cyprus", "Turkey")A <- connect("Cyprus", "Lebanon")A <- connect("Norway", "Iceland")## remove connections between South American and FranceA <- connect("Suriname", "France", add = FALSE)A <- connect("Brazil", "France", add = FALSE)
We should look at the revised adjacency matrix:
graph <- st\_geometry( edges(A, shape = world) )ogpar <- par(mar = rep(0, 4))plot(world\_geom, lwd = .1)plot(graph, add = TRUE, type = 'p')plot(graph, add = TRUE, type = 'l')
par(ogpar)
Sometimes it can help to examine the edge list interactively using aproper geographic information system like QGIS. For those who arefamiliar with (non-R) GIS software, you can save the edges list as aGeoPackage and then open it up in your GIS to examine the connections‘by hand’ with a base map or other data:
E <- edges(A, shape = world)st\_write(E, "world.gpkg", layer = "edge list")
Non-spatial regressionFitting regression models with geostan is similar to using base R’sglm function: the user provides a model formula, data, and the modelfamily or distribution. We can fit a normal linear model using thestan_glm function:
fit\_lm <- stan\_glm(lifeExp ~ log(gdpPercap), data = world, quiet = TRUE)
And we can examine parameter estimates by printing to the console:
print(fit\_lm)## Spatial Model Results ## Formula: lifeExp ~ log(gdpPercap)## Spatial method (outcome): none ## Likelihood function: gaussian ## Link function: identity ## Residual Moran Coefficient: NA ## WAIC: 977.33 ## Observations: 162 ## Data models (ME): none## Inference for Stan model: foundation.## 4 chains, each with iter=2000; warmup=1000; thin=1; ## post-warmup draws per chain=1000, total post-warmup draws=4000.## ## mean se\_mean sd 2.5% 20% 50% 80% 97.5% n\_eff## intercept 20.660 0.089 2.934 15.106 18.148 20.705 23.117 26.431 1082## log(gdpPercap) 5.499 0.010 0.317 4.879 5.235 5.498 5.770 6.108 1081## sigma 4.902 0.007 0.268 4.408 4.676 4.894 5.119 5.468 1458## Rhat## intercept 1.002## log(gdpPercap) 1.002## sigma 1.001## ## Samples were drawn using NUTS(diag\_e) at Thu Aug 1 12:33:33 2024.## For each parameter, n\_eff is a crude measure of effective sample size,## and Rhat is the potential scale reduction factor on split chains (at ## convergence, Rhat=1).
The output printed to the console provides a summary of the posteriorprobability distributions of the model parameters. The distributions canalso be visualized using plot(fit_lm):
plot(fit\_lm)## `stat\_bin()` using `bins = 30`. Pick better value with `binwidth`.
The mean of the distribution is reported in the mean column. For thosewho are more familiar with concepts from sampling theory, the mean maybe understood as the estimate of the parameter. Each distribution’sstandard deviation is found in the sd column; this describes the widthof the posterior distribution. The sd is analogous to the standarderror of the estimate. The quantiles also summarize the width of theposterior distributions; the 2.5% and 97.5% values form a 95% credibleinterval for the parameter value.
MCMC outputThe effective sample size (ESS), n_eff, tells us how many independentMCMC samples the inference is based on after adjusting for serialautocorrelation in the MCMC samples. This is an important quantity topay attention to and generally one might like to see these numbers above400 (or around 100 samples per MCMC chain). The standard error of themean, se_mean, reports how much MCMC sampling error to expect in themean (se_mean is calculated using n_eff). The R-hat statistic,Rhat, should always be very close to 1, preferably less than 1.01. TheR-hat diagnostic tests that the MCMC chains are all depicting the samedistribution. If they diverge from one another, it either means that youneed to draw a higher number of MCMC samples (run the chains for longer)or that there is a problem fitting the model to your data.
By default, geostan models run four independent MCMC chains for 3,000iterations each, half of which is discarded as warm-up. The number ofiterations is controlled by the iter argument, the default beingiter = 3e3. For some models this may be too low and you will want toincrease this. Other times this might be more than is needed in whichcase you can reduce the computation time by using fewer iterations. Whatmatters most is not your number of iterations but your ESS and R-hatstatistics. When it comes to reporting results, it is generally best touse at least the default of four MCMC chains (chains = 4).
MethodsA number of familiar methods are available for working with geostanmodels including fitted, resid, and predict.
The fitted method returns a data.frame with summaries of the fittedvalues. The probability distribution for each fitted value is summarizedby its posterior mean, standard deviation, and quantiles:
fdf <- fitted(fit\_lm)head(fdf)## mean sd 2.5% 20% 50% 80% 97.5%## fitted[1] 70.22990 0.3969583 69.46571 69.88939 70.23141 70.55628 71.00631## fitted[2] 63.46359 0.5884022 62.31542 62.98235 63.45552 63.94977 64.64037## fitted[3] 79.33702 0.6206999 78.07451 78.81257 79.34308 79.86151 80.53088## fitted[4] 80.36367 0.6675295 79.01816 79.79764 80.37163 80.93076 81.65248## fitted[5] 76.02492 0.4883834 75.05750 75.61649 76.03594 76.44239 76.95946## fitted[6] 67.88820 0.4339864 67.04641 67.52555 67.88695 68.24651 68.74736
The resid method behaves similarly. Examining the Moran scatter plotusing the residuals shows a moderate degree of positive SA as well assome skewness:
rdf <- resid(fit\_lm)moran\_plot(rdf$mean, A)
Spatial regressionOptions for spatial regression models currently include conditionalautoregressive (CAR), simultaneous autoregressive (SAR/spatial error),and eigenvector spatial filtering (ESF). For count data, commonvariations on the intrinsic autoregressive (ICAR) model are alsoavailable.
All of the spatial models require at least a spatial weights matrix asinput. All additional requirements for data preparation are handled bygeostan’s prep_ functions: prep_car_data, prep_icar_data,prep_sar_data.
The make_EV function is used to create Moran’s eigenvectors for ESFregression; if you want to create your own eigenvectors (say, followingthe PCNM method) you can provide those directly to the ESF model (see?stan_esf).
For the CAR model, we always provide the binary adjacency matrix asinput to prep_car_data. See the prep_car_data documentation foroptions. Here we will fit an intercept-only CAR model to the lifeexpectancy data:
cars <- prep\_car\_data(A)fit\_car <- stan\_car(lifeExp ~ 1, data = world, car\_parts = cars, iter = 1e3, quiet = TRUE)print(fit\_car)## Spatial Model Results ## Formula: lifeExp ~ 1## Spatial method (outcome): CAR ## Likelihood function: auto\_gaussian ## Link function: identity ## Residual Moran Coefficient: -0.356467 ## WAIC: 994.65 ## Observations: 162 ## Data models (ME): none## Inference for Stan model: foundation.## 4 chains, each with iter=1000; warmup=500; thin=1; ## post-warmup draws per chain=500, total post-warmup draws=2000.## ## mean se\_mean sd 2.5% 20% 50% 80% 97.5% n\_eff Rhat## intercept 70.261 0.081 2.664 65.007 68.274 70.280 72.325 76.104 1080 1.002## car\_rho 0.981 0.000 0.010 0.958 0.973 0.983 0.990 0.996 1245 1.003## car\_scale 7.827 0.011 0.449 7.000 7.441 7.806 8.202 8.761 1564 1.001## ## Samples were drawn using NUTS(diag\_e) at Thu Aug 1 12:33:38 2024.## For each parameter, n\_eff is a crude measure of effective sample size,## and Rhat is the potential scale reduction factor on split chains (at ## convergence, Rhat=1).
Notice that using iter = 1000 was more than adequate for inference inthis case.
The CAR model has a spatial dependence parameter car_rho. Thisparameter does not have an interpretation similar to a correlationcoefficient, and it is often near 1; this is not a problem unless onemisinterprets it or desires a value similar to the correlationcoefficient. The spatial dependence parameter in the SAR model doesprovide that kind of interpretation.
A filtering approachReturning to the correlation coefficient estimated above, one way toadjust our estimate for spatial dependence is to filter out the spatialtrend from each of the two variables and then calculate the correlationcoefficient using the detrended values (Chun and Griffith 2012, 71).This spatial ‘filtering’ or ‘pre-whitening’ method is not particularlycommon in practice but its a good trick to know given the familiarity ofthe correlation coefficient. We will use it here to demonstrate somebasic features of the software.
The spatial trend term can be extracted from any spatial geostan modelusing the spatial method.
theta <- spatial(fit\_car)$meanpars <- map\_pars(theta)ogpar <- par(mar = rep(0, 4))plot(st\_geometry(world), col = pars$col, lwd = .2)legend("left", fill = pars$pal, title = 'Spatial trend (LE)', legend = pars$lbls, bty = 'n' )
par(ogpar)
We can obtain detrended values most simply by taking the residuals froman intercept-only spatial model:
```
``
Usingcor.test` with those provides an estimate of correlation adjustedfor spatial autocorrelation:
```
``` The adjusted estimate of .59 is considerably different from the naiveestimate of .80 and is outside the naive confidence intervals. (Theadjusted estimate is .62 if we use SAR models.)
A bivariate modelHere we will use the SAR model to illustrate its use. Fitting thespatial error or SAR model requires nearly the same steps as above.
Unlike prep_car_data, be sure to row-standardize the adjacency matrixbefore passing it to prep_sar_data.
W <- row\_standardize(A)sars <- prep\_sar\_data(W)
When fitting the model, we are going to add centerx = TRUE to centerthe covariate. (Internally this will callcenter(x, center = TRUE, scale = FALSE).) This will not changecoefficient estimates but it does often improve MCMC samplingefficiency, sometimes considerably so. It does change interpretation ofthe intercept: the intercept will be an estimate of the average lifeexpectancy (or the expected life expectancy when all covariates are attheir average values).
fit\_sar <- stan\_sar(lifeExp ~ log(gdpPercap), data = world, sar\_parts = sars, centerx = TRUE, iter = 1e3, quiet = TRUE)
Lets plot the results this time:
plot(fit\_sar)## `stat\_bin()` using `bins = 30`. Pick better value with `binwidth`.
The spatial dependence parameter is around 0.7, which indicatesmoderately strong SA. The mean life expectancy is about 71 (probablysomewhere between about 69 and 74). And the coefficient for log GDP isaround 4 (or somewhere between 3 and 5). The residual variation has astandard deviation of around 3.6 years.
If we scale both variables before fitting the bivariate spatialregression model (so that their variances both equal 1) then we getapproximately the same estimate as the adjusted correlation coefficient(above). The credible interval is slightly wider because uncertainty inrho is (appropriately) incorporated here:
world <- transform(world, sx = scale(log(gdpPercap), scale = T, center = T), sy = scale(lifeExp, scale = T, center = T) )fit\_scaled <- stan\_sar(sy ~ sx, data = world, sar\_parts = sars, iter = 1e3, quiet = TRUE)print(fit\_scaled)## Spatial Model Results ## Formula: sy ~ sx## Spatial method (outcome): SAR ## Likelihood function: auto\_gaussian ## Link function: identity ## Residual Moran Coefficient: -0.0555795 ## WAIC: 224.16 ## Observations: 162 ## Data models (ME): none## Inference for Stan model: foundation.## 4 chains, each with iter=1000; warmup=500; thin=1; ## post-warmup draws per chain=500, total post-warmup draws=2000.## ## mean se\_mean sd 2.5% 20% 50% 80% 97.5% n\_eff Rhat## intercept 0.017 0.003 0.121 -0.213 -0.081 0.015 0.110 0.259 1374 1.001## sx 0.584 0.002 0.062 0.462 0.533 0.586 0.636 0.706 1550 1.001## sar\_rho 0.703 0.002 0.057 0.577 0.656 0.708 0.753 0.801 1224 1.000## sar\_scale 0.438 0.001 0.027 0.391 0.415 0.436 0.461 0.496 1541 1.000## ## Samples were drawn using NUTS(diag\_e) at Thu Aug 1 12:33:49 2024.## For each parameter, n\_eff is a crude measure of effective sample size,## and Rhat is the potential scale reduction factor on split chains (at ## convergence, Rhat=1).
Predicted valuesWe can visualize the model results by plotting the expected lifeexpectancy across the full range of GDP per capita. We use the predictfunction for this. As input, it requires our fitted model and adata.frame with covariate values.
We will start by creating a data.frame with GDP per capita values thatspan from the minimum to maximum values in the world data:
gdp <- range(world$gdpPercap)min\_gdp <- gdp[1]max\_gdp <- gdp[2]pdf <- data.frame(gdpPercap = seq(min\_gdp, max\_gdp, length.out = 200))
The column names in this data.frame have to match the variable namesthat were present in the data that we first provided to the model. Inthis case, the name of the columns should match those from the worlddata. Likewise, we provide the new GDP data on its original(un-transformed) scale, just as we did when we fit the model usingstan_sar (the log transformation will be applied by predict becauseit is specified in the model formula). Because we centered the covariateusing the centerx = TRUE argument, we will also allow the predictfunction to handle the centering automatically using information that isstored with the fitted model (stan_sar$x_center).
Now we pass this new data to predict:
preds <- predict(fit\_sar, newdata = pdf)
The output includes our pdf data plus some new columns. The newcolumns provide a summary of the predicted values. As usual, the meanis the estimate and the estimate is accompanied by other values that canbe used to taken as credible intervals for the predicted value. Theoutput reflects uncertainty in the model parameter estimates.
head(preds)## gdpPercap mean sd 2.5% 20% 50% 80% 97.5%## 1 597.1352 60.10281 1.525231 57.31512 58.81438 60.02981 61.40785 63.23217## 2 1201.4715 62.90912 1.324012 60.47457 61.78458 62.84984 64.02423 65.60880## 3 1805.8079 64.54460 1.222436 62.29812 63.50132 64.50656 65.56321 67.07337## 4 2410.1442 65.70331 1.159626 63.56161 64.71873 65.67025 66.65184 68.13579## 5 3014.4805 66.60137 1.117121 64.51219 65.64677 66.57454 67.52918 68.94410## 6 3618.8169 67.33477 1.086904 65.28404 66.40366 67.29468 68.21982 69.56564
These ‘predicted’ values represent the expectation of the outcomevariable at the given level of the covariates. So we would expect actualobservations to form a cloud of points around the ‘predicted’ values. Tocalculate these predicted values, the predict function only includescovariates and the intercept, it does not include any spatialautocorrelation components. Its purpose is to examine implications ofthe coefficient estimates on recognizable scales of variation, not topredict values for particular places. (The log-linear model can also beinterpreted in terms of percent changes in the covariate, such as ’a 10%increase in GDP per capita, e.g., from 10,000 to 11,000, is associatedwith around 4 * log(11/10) = 0.38 additional years of life expectancyon average.)
```
2.5%), max(preds$97.5%))plot(preds$gdpPercap, preds$mean, t = 'l', ylim = yrange, axes = F, xlab = "GDP per capita ($1,000s)", ylab = "Life expectancy")axis(1)axis(2)# add credible intervalslines(preds$gdpPercap, preds$2.5%, lty = 3)lines(preds$gdpPercap, preds$97.5%, lty = 3)``` Per this dataset, about 50% of the world population lives in countrieswith GDP per capita below $12,300.
Future work and supportYou can submit any questions, requests, or issues on the package issuespage or thediscussionspage. geostanis still actively being developed so users are encouraged to check thepackage news page forupdates.
If you are interesting contributing to the package you are encouraged tosend an e-mail to the author or use the discussions page. You can submita pull request with any bug fixes. Contributions that would make thepackage more useful to fields other than geostan’s current focus(human geography and public health), such as ecology, would beespecially welcome.
Appendix ```
``` ReferencesChun, Yongwan, and Daniel A Griffith. 2012. “Spatial Statistics andGeostatistics: Theory and Applications for Geographic InformationScience and Technology.”
Donegan, Connor. 2021. “Building Spatial Conditional Autoregressive(CAR) Models in the Stan Programming Language.” https://osf.io/3ey65/.
———. 2022. “Geostan: An R Package for Bayesian Spatial Analysis.”Journal of Open Source Software 7 (79): 4716.https://doi.org/10.21105/joss.04716.
———. 2024. “Plausible Reasoning and Spatial-Statistical Theory: ACritique of Recent Writings on ‘Spatial Confounding’.” GeographicalAnalysis Early view. https://doi.org/10.1111/gean.12408.
Gabry, Jonah, Ben Goodrich, Martin Lysy, and Andrew Johnson. 2024.Rstantools: Tools for Developing R Packages Interfacing with ’Stan’.https://CRAN.R-project.org/package=rstantools.
Griffith, Daniel A, and Pedro R Peres-Neto. 2006. “Spatial Modeling inEcology: The Flexibility of Eigenfunction Spatial Analyses.” Ecology87 (10): 2603–13.
Stan Development Team. 2023. Stan User’s Guide. https://mc-stan.org.
To leave a comment for the author, please follow the link and comment on their blog: r-spatial.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.Continue reading: Spatial analysis with geostan