# Chapter 16 Spatial Regression

Even though it may be tempting to focus on interpreting the map pattern of an areal support response variable of interest, the pattern may largely derive from covariates (and their functional forms), as well as the respective spatial footprints of the variables in play. Spatial autoregressive models in two dimensions began without covariates and with clear links to time series (Whittle 1954). Extensions included tests for spatial autocorrelation in linear model residuals, and models applying the autoregressive component to the response or the residuals, where the latter matched the tests for residuals (Cliff and Ord 1972, 1973). These “lattice” models of areal data typically express the dependence between observations using a graph of neighbours in the form of a contiguity matrix.

Of course, handling a spatial correlation structure in a generalised least squares model or a (generalized) linear or nonlinear mixed effects model such as those provided in the nlme and many other packages does not have to use a graph of neighbours (Pinheiro and Bates 2000). These models are also spatial regression models, using functions of the distance between observations, and fitted variograms to model the spatial autocorrelation present; such models have been held to yield a clearer picture of the underlying processes (Wall 2004), building on geostatistics. For example, the glmmTMB package successfully uses this approach to spatial regression (Brooks et al. 2017). Here we will only consider spatial regression using spatial weights matrices.

## 16.1 Markov random field and multilevel models with spatial weights

There is a large literature in disease mapping using conditional autoregressive (CAR) and intrinsic CAR (ICAR) models in spatially structured random effects. These extend to multilevel models, in which the spatially structured random effects may apply at different levels of the model (Roger S. Bivand et al. 2017). In order to try out some of the variants, we need to remove the no-neighbour observations from the tract level, and from the model output zone aggregated level, in two steps as reducing the tract level induces a no-neighbour outcome at the model output zone level. Many of the model estimating functions take family= arguments, and fit generalized linear mixed effects models with per-observation spatial random effects structured using a Markov random field representation of relationships between neighbours. In the multilevel case, the random effects may be modelled at the group level, which is the case presented in the following examples.

We follow Gómez-Rubio (2019) in summarizing Pinheiro and Bates (2000) and McCulloch and Searle (2001) to describe the mixed-effects model representation of spatial regression models. In a Gaussian linear mixed model setting, a random effect $$u$$ is added to the model, with response $$Y$$, fixed covariates $$X$$, their coefficients $$\beta$$ and error term $$\varepsilon_i \sim N(0, \sigma^2), i=1,\dots, n$$:

$Y = X \beta + Z u + \varepsilon$ $$Z$$ is a fixed design matrix for the random effects. If there are $$n$$ random effects, it will be an $$n \times n$$ identity matrix, if instead the observations are aggregated into $$m$$ groups, so with $$m < n$$ random effects, it will be an $$n \times m$$ matrix showing which group each observation belongs to. The random effects are modelled as a multivariate Normal distribution $$u \sim N(0, \sigma^2_u \Sigma)$$, and $$\Sigma$$ is the square variance-covariance matrix of the random effects.

A division has grown up, possibly unhelpfully, between scientific fields using CAR models (Besag 1974), and simultaneous autoregressive models (SAR) (Ord 1975; Hepple 1976). Although CAR and SAR models are closely related, these fields have found it difficult to share experience of applying similar models, often despite referring to key work summarising the models (Ripley 1981, 1988; Cressie 1993). Ripley gives the SAR variance as (1981, 89), here shown as the inverse $$\Sigma^{-1}$$ (also known as the precision matrix):

$\Sigma^{-1} = [(I - \rho W)'(I - \rho W)]$

where $$\rho$$ is a spatial autocorrelation parameter and $$W$$ is a nonsingular spatial weights matrix that represents spatial dependence. The CAR variance is:

$\Sigma^{-1} = (I - \rho W)$ where $$W$$ is a symmetric and strictly positive definite spatial weights matrix. In the case of the intrinsic CAR model, avoiding the estimation of a spatial autocorrelation parameter, we have:

$\Sigma^{-1} = M = \mathrm{diag}(n_i) - W$ where $$W$$ is a symmetric and strictly positive definite spatial weights matrix as before and $$n_i$$ are the row sums of $$W$$. The Besag-York-Mollié model includes intrinsic CAR spatially structured random effects and an unstructured random effects. The Leroux model combines matrix components for unstructured and spatially structured random effects, where the spatially structured random effects are taken as following an intrinsic CAR specification:

$\Sigma^{-1} = [(1 - \rho) I_n + \rho M]$ References to the definitions of these models may be found in V. Gómez-Rubio (2020), and estimation issues affecting the Besag-York-Mollié and Leroux models are reviewed by Gerber and Furrer (2015).

More recent books expounding the theoretical bases for modelling with areal data simply point out the similarities between SAR and CAR models in relevant chapters (Gaetan and Guyon 2010; Lieshout 2019); the interested reader is invited to consult these sources for background information.

### 16.1.1 Boston house value data set

Here we shall use the Boston housing data set, which has been restructured and furnished with census tract boundaries (R. Bivand 2017). The original data set used 506 census tracts and a hedonic model to try to estimate willingness to pay for clean air. The response was constructed from counts of ordinal answers to a 1970 census question about house value. The response is left and right censored in the census source and has been treated as Gaussian. The key covariate was created from a calibrated meteorological model showing the annual nitrogen oxides (NOX) level for a smaller number of model output zones. The numbers of houses responding also varies by tract and model output zone. There are several other covariates, some measured at the tract level, some by town only, where towns broadly correspond to the air pollution model output zones.

We can start by reading in the 506 tract data set from spData (R. Bivand, Nowosad, and Lovelace 2021), and creating a contiguity neighbour object and from that again a row standardized spatial weights object.

library(sf)
library(spData)
boston_506 <- st_read(system.file("shapes/boston_tracts.shp", package = "spData")[1])
# Reading layer boston_tracts' from data source
#   /home/edzer/R/x86_64-pc-linux-gnu-library/4.0/spData/shapes/boston_tracts.shp'
#   using driver ESRI Shapefile'
# Simple feature collection with 506 features and 36 fields
# Geometry type: POLYGON
# Dimension:     XY
# Bounding box:  xmin: -71.5 ymin: 42 xmax: -70.6 ymax: 42.7
# Geodetic CRS:  NAD27
nb_q <- spdep::poly2nb(boston_506)
lw_q <- spdep::nb2listw(nb_q, style = "W")

If we examine the median house values, we find that those for censored values have been assigned as missing, and that 17 tracts are affected.

table(boston_506$censored) # # left no right # 2 489 15 summary(boston_506$median)
#    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's
#    5600   16800   21000   21749   24700   50000      17

Next, we can subset to the remaining 489 tracts with non-censored house values, and the neighbour object to match. The neighbour object now has one observation with no neighbours.

boston_506$CHAS <- as.factor(boston_506$CHAS)
boston_489 <- boston_506[!is.na(boston_506$median),] nb_q_489 <- spdep::poly2nb(boston_489) lw_q_489 <- spdep::nb2listw(nb_q_489, style = "W", zero.policy = TRUE) The NOX_ID variable specifies the upper level aggregation, letting us aggregate the tracts to air pollution model output zones. We can create aggregate neighbour and row standardized spatial weights objects, and aggregate the NOX variable taking means, and the CHAS Charles River dummy variable for observations on the river. Here we follow the principles outlined in section 5.3.1 for spatially extensive and intensive variables; neither NOX nor CHAS can be summed as they are not count variables. agg_96 <- list(as.character(boston_506$NOX_ID))
boston_96 <- aggregate(boston_506[, "NOX_ID"], by = agg_96, unique)
nb_q_96 <- spdep::poly2nb(boston_96)
lw_q_96 <- spdep::nb2listw(nb_q_96)
boston_96$NOX <- aggregate(boston_506$NOX, agg_96, mean)$x boston_96$CHAS <- aggregate(as.integer(boston_506$CHAS)-1, agg_96, max)$x

The response is aggregated using the weightedMedian() function in matrixStats, and midpoint values for the house value classes. Counts of houses by value class were punched to check the published census values, which can be replicated using weightedMedian() at the tract level. Here we find two output zones with calculated weighted medians over the upper census question limit of USD 50,000, and remove them subsequently as they also are affected by not knowing the appropriate value to insert for the top class by value. This is a case of spatially extensive aggregation, for which the summation of counts is appropriate:

nms <- names(boston_506)
ccounts <- 23:31
for (nm in nms[c(22, ccounts, 36)]) {
boston_96[[nm]] <- aggregate(boston_506[[nm]], agg_96, sum)$x } br2 <- c(3.50, 6.25, 8.75, 12.50, 17.50, 22.50, 30.00, 42.50, 60.00)*1000 counts <- as.data.frame(boston_96)[, nms[ccounts]] f <- function(x) matrixStats::weightedMedian(x = br2, w = x, interpolate = TRUE) boston_96$median <- apply(counts, 1, f)
is.na(boston_96$median) <- boston_96$median > 50000
summary(boston_96$median) # Min. 1st Qu. Median Mean 3rd Qu. Max. NA's # 9009 20417 23523 25263 30073 49496 2 Before subsetting, we aggregate the remaining covariates by weighted mean using the tract population counts punched from the census (R. Bivand 2017); these are spatially intensive variables, not count data. boston_94 <- boston_96[!is.na(boston_96$median),]
nb_q_94 <- spdep::subset.nb(nb_q_96, !is.na(boston_96$median)) lw_q_94 <- spdep::nb2listw(nb_q_94, style="W") We now have two data sets at each level, at the lower, census tract level, and at the upper, air pollution model output zone level, one including the censored observations, the other excluding them. boston_94a <- aggregate(boston_489[,"NOX_ID"], list(boston_489$NOX_ID), unique)
nb_q_94a <- spdep::poly2nb(boston_94a)
NOX_ID_no_neighs <- boston_94a$NOX_ID[which(spdep::card(nb_q_94a) == 0)] boston_487 <- boston_489[is.na(match(boston_489$NOX_ID, NOX_ID_no_neighs)),]
boston_93 <- aggregate(boston_487[, "NOX_ID"], list(ids = boston_487$NOX_ID), unique) row.names(boston_93) <- as.character(boston_93$NOX_ID)

### 16.2.2 IID and CAR random effects with hglm

The same model may be estimated using the hglm package (Alam, Ronnegard, and Shen 2019), which also permits the modelling of discrete responses, this time using an extra one-sided formula to express the random effects term:

suppressPackageStartupMessages(library(hglm))
suppressWarnings(HGLM_iid <- hglm(fixed = form, random = ~1 | NOX_ID, data = boston_487,
family = gaussian()))
boston_93$HGLM_re <- unname(HGLM_iid$ranef)

The same package has been extended to spatially structured SAR and CAR random effects, for which a sparse spatial weights matrix is required (Alam, Rönnegård, and Shen 2015); we choose binary spatial weights:

W <- as(spdep::nb2listw(nb_q_93, style = "B"), "CsparseMatrix")

We fit a CAR model at the upper level, using the rand.family= argument, where the values of the indexing variable NOX_ID match the row names of $$W$$:

suppressWarnings(HGLM_car <- hglm(fixed = form, random = ~ 1 | NOX_ID, data = boston_487,
family = gaussian(), rand.family = CAR(D=W)))
boston_93$HGLM_ss <- HGLM_car$ranef[,1]

### 16.2.3 SAR random effects with HSAR

The HSAR package (Dong, Harris, and Mimis 2020) is restricted to the Gaussian response case (G. Dong and Harris 2015; G. Dong et al. 2015), and requires the specification of a sparse design matrix mapping the upper level entities onto lower level entities:

library(HSAR)
library(MatrixModels, warn.conflicts = FALSE)
Z <- as(model.Matrix(~ -1 + as.factor(NOX_ID), data = boston_487, sparse = TRUE),
"dgCMatrix")

hsar() fits an upper level SAR random effect using MCMC; if W= is a lower level weights matrix rather than NULL, it will also fit a lower level spatial econometrics-style spatial lag model, adding the lower level spatially lagged response to the model:

set.seed(123)
suppressWarnings(HSAR_ss <- hsar(form, data=boston_487, W = NULL, M = W, Delta = Z,
burnin = 2000, Nsim = 12000, thinning = 2))
boston_93$HSAR_ss <- HSAR_ss$Mus[1,]

### 16.2.4 IID and ICAR random effects with R2BayesX

The R2BayesX package (Umlauf et al. 2021) provides flexible support for structured additive regression models, including spatial multilevel models. The models include an IID unstructured random effect at the upper level using the "re" specification in the sx() model term (Umlauf et al. 2015); we choose the "MCMC" method:

suppressPackageStartupMessages(library(R2BayesX))
BX_iid <- bayesx(update(form, . ~ . + sx(NOX_ID, bs = "re")), family = "gaussian",
data = boston_487, method = "MCMC", iterations = 12000, burnin = 2000, step = 2, seed = 123)
boston_93$BX_re <- BX_iid$effects["sx(NOX_ID):re"][[1]]$Mean and the "mrf" (Markov Random Field) spatially structured intrinsic CAR random effect specification based on a graph derived from converting a suitable "nb" object for the upper level. The "region.id" attribute of the "nb" object needs to contain values corresponding to the indexing variable in the sx() effects term, to facilitate the internal construction of design matrix $$Z$$: RBX_gra <- nb2gra(nb_q_93) all.equal(row.names(RBX_gra), attr(nb_q_93, "region.id")) # [1] TRUE As we saw above in the intrinsic CAR model definition, the counts of neighbours are entered on the diagonal, but the current implementation uses a dense, not sparse, matrix: all.equal(unname(diag(RBX_gra)), spdep::card(nb_q_93)) # [1] TRUE The sx() model term continues to include the indexing variable, and now passes through the intrinsic CAR precision matrix: BX_mrf <- bayesx(update(form, . ~ . + sx(NOX_ID, bs = "mrf", map = RBX_gra)), family = "gaussian", data = boston_487, method = "MCMC", iterations = 12000, burnin = 2000, step = 2, seed = 123) boston_93$BX_ss <- BX_mrf$effects["sx(NOX_ID):mrf"][[1]]$Mean

### 16.2.5 IID, ICAR and Leroux random effects with INLA

R. Bivand, Gómez-Rubio, and Rue (2015) and V. Gómez-Rubio (2020) present the use of the INLA package (Rue et al. 2021) and the inla() model fitting function with spatial regression models:

suppressPackageStartupMessages(library(INLA))

Although differing in details, the approach by updating the fixed model formula with an unstructured random effects term is very similar to that seen above:

INLA_iid <- inla(update(form, . ~ . + f(NOX_ID, model = "iid")), family = "gaussian",
data = boston_487)
boston_93$INLA_re <- INLA_iid$summary.random$NOX_ID$mean

As with most implementations, care is needed to match the indexing variable with the spatial weights; in this case using indices $$1, \dots, 93$$ rather than the NOX_ID variable directly:

ID2 <- as.integer(as.factor(boston_487$NOX_ID)) The same sparse binary spatial weights matrix is used, and the intrinsic CAR representation is constructed internally: INLA_ss <- inla(update(form, . ~ . + f(ID2, model = "besag", graph = W)), family = "gaussian", data = boston_487) boston_93$INLA_ss <- INLA_ss$summary.random$ID2$mean The sparse Leroux representation as given by V. Gómez-Rubio (2020) can be constructed in the following way: M <- Diagonal(nrow(W), rowSums(W)) - W Cmatrix <- Diagonal(nrow(M), 1) - M This model can be estimated using the "generic1" model with the specified precision matrix: INLA_lr <- inla(update(form, . ~ . + f(ID2, model = "generic1", Cmatrix = Cmatrix)), family = "gaussian", data = boston_487) boston_93$INLA_lr <- INLA_lr$summary.random$ID2$mean ### 16.2.6 ICAR random effects with mgcv::gam() In a very similar way, the gam() function in the mgcv package (Wood 2022) can take an "mrf" term using a suitable "nb" object for the upper level. In this case the "nb" object needs to have the contents of the "region.id" attribute copied as the names of the neighbour list components, and the indexing variable needs to be a factor (Wood 2017): library(mgcv) names(nb_q_93) <- attr(nb_q_93, "region.id") boston_487$NOX_ID <- as.factor(boston_487$NOX_ID) The specification of the spatially structured term again differs in details from those above, but achieves the same purpose. The "REML" method of bayesx() gives the same results as gam() using "REML" in this case: GAM_MRF <- gam(update(form, . ~ . + s(NOX_ID, bs = "mrf", xt = list(nb = nb_q_93))), data = boston_487, method = "REML") The upper level random effects may be extracted by predicting terms; as we can see, the values in all lower-level tracts belonging to the same upper-level air pollution model output zones are identical: ssre <- predict(GAM_MRF, type = "terms", se = FALSE)[, "s(NOX_ID)"] all(sapply(tapply(ssre, list(boston_487$NOX_ID), c), function(x) length(unique(x)) == 1))
# [1] TRUE

so we can return the first value for each upper-level unit:

boston_93$GAM_ss <- aggregate(ssre, list(boston_487$NOX_ID), head, n=1)\$x

### 16.2.7 Upper level random effects: summary

In the cases of hglm(), bayesx(), inla() and gam(), we could also model discrete responses without further major difficulty, and bayesx(), inla() and gam() also facilitate the generalization of functional form fitting for included covariates.

Unfortunately, the coefficient estimates for the air pollution variable for these multilevel models are not helpful. All are negative as expected, but the inclusion of the model output zone level effects, IID or spatially structured, makes it is hard to disentangle the influence of the scale of observation from that of covariates observed at that scale rather than at the tract level.

Figure 16.1 shows that the air pollution model output zone level IID random effects are very similar across the four model fitting functions reported. In all the maps, the central downtown zones have stronger negative random effect values, but strong positive values are also found in close proximity; suburban areas take values closer to zero.

library(tmap, warn.conflicts=FALSE)
tm_shape(boston_93) + tm_fill(c("MLM_re", "HGLM_re", "INLA_re", "BX_re"), midpoint = 0,
title = "IID")  + tm_facets(free.scales = FALSE) + tm_borders(lwd = 0.3, alpha = 0.4) +
tm_layout(panel.labels = c("lmer", "hglm", "inla", "bayesx"))

Figure 16.2 shows that the spatially structured random effects are also very similar to each other, with the "SAR" spatial smooth being perhaps a little smoother than the "CAR" smooths when considering the range of values taken by the random effect term.

tm_shape(boston_93) + tm_fill(c("HGLM_ss", "HSAR_ss", "INLA_lr", "INLA_ss", "BX_ss",
"GAM_ss"), midpoint = 0, title = "SSRE")  + tm_facets(free.scales = FALSE) +
tm_borders(lwd = 0.3, alpha = 0.4) + tm_layout(panel.labels = c("hglm CAR", "hsar SAR",
"inla Leroux", "inla ICAR", "bayesx ICAR", "gam ICAR"))

Although there is still a great need for more thorough comparative studies of model fitting functions for spatial regression including multilevel capabilities, there has been much progress over recent years. Vranckx, Neyens, and Faes (2019) offer a recent comparative survey of disease mapping spatial regression, typically set in a Poisson regression framework offset by an expected count. In Roger S Bivand and Gómez-Rubio (2021), methods for estimating spatial survival models using spatial weights matrices are compared with spatial probit models.

## 16.3 Exercises

1. Construct a multilevel dataset using the Athens housing data in HSAR, as in the vignette: https://cran.r-project.org/web/packages/HSAR/vignettes/PropertiesAthens.html. At which point do the municipality department attribute values get copied out to all the point observations within each municipality department?
2. Create neighbour objects at both levels. Test greensp for spatial autocorrelation at the upper level, and then at the lower level. What has been the chief consequence of copying out the area of green spaces in square meters for the municipality departments to the point support property level?
3. Using the formula object from the vignette, assess whether adding the copied out upper level variables seems sensible. Use mgcv::gam() to fit a linear mixed effects model (IID of num_dep identifying the municipality departments) using just the lower level variables and the lower and upper level variables. Do your conclusions differ?
4. Complete the analysis by replacing the IID random effects with an "mrf" Markov random field and the contiguity neighbour object created above. Do you think that it is reasonable to for example draw any conclusions based on the municipality department level variables such as greensp`?

### References

Alam, Moudud, Lars Ronnegard, and Xia Shen. 2019. Hglm: Hierarchical Generalized Linear Models. https://CRAN.R-project.org/package=hglm.

Alam, Moudud, Lars Rönnegård, and Xia Shen. 2015. “Fitting Conditional and Simultaneous Autoregressive Spatial Models in Hglm.” The R Journal 7 (2): 5–18. https://doi.org/10.32614/RJ-2015-017.

Bates, Douglas, Martin Maechler, Ben Bolker, and Steven Walker. 2021. Lme4: Linear Mixed-Effects Models Using Eigen and S4. https://github.com/lme4/lme4/.

Besag, Julian. 1974. “Spatial Interaction and the Statistical Analysis of Lattice Systems.” Journal of the Royal Statistical Society. Series B (Methodological) 36: pp. 192–236.

Bivand, Roger. 2017. “Revisiting the Boston Data Set — Changing the Units of Observation Affects Estimated Willingness to Pay for Clean Air.” REGION 4 (1): 109–27. https://doi.org/10.18335/region.v4i1.107.

Bivand, Roger, Virgilio Gómez-Rubio, and Håvard Rue. 2015. “Spatial Data Analysis with R-Inla with Some Extensions.” Journal of Statistical Software, Articles 63 (20): 1–31. https://doi.org/10.18637/jss.v063.i20.

Bivand, Roger, Jakub Nowosad, and Robin Lovelace. 2021. SpData: Datasets for Spatial Analysis. https://nowosad.github.io/spData/.

Bivand, Roger S, and Virgilio Gómez-Rubio. 2021. “Spatial Survival Modelling of Business Re-Opening After Katrina: Survival Modelling Compared to Spatial Probit Modelling of Re-Opening Within 3, 6 or 12 Months.” Statistical Modelling 21 (1-2): 137–60. https://doi.org/10.1177/1471082X20967158.

Bivand, Roger S., Zhe Sha, Liv Osland, and Ingrid Sandvig Thorsen. 2017. “A Comparison of Estimation Methods for Multilevel Models of Spatially Structured Data.” Spatial Statistics. https://doi.org/10.1016/j.spasta.2017.01.002.

Brooks, Mollie E., Kasper Kristensen, Koen J. van Benthem, Arni Magnusson, Casper W. Berg, Anders Nielsen, Hans J. Skaug, Martin Maechler, and Benjamin M. Bolker. 2017. “glmmTMB Balances Speed and Flexibility Among Packages for Zero-Inflated Generalized Linear Mixed Modeling.” The R Journal 9 (2): 378–400. https://journal.r-project.org/archive/2017/RJ-2017-066/index.html.

Cliff, A. D., and J. K. Ord. 1973. Spatial Autocorrelation. London: Pion.

Cliff, A., and J.K. Ord. 1972. “Testing for Spatial Autocorrelation Among Regression Residuals.” Geographical Analysis 4: 267–84.

Cressie, N. A. C. 1993. Statistics for Spatial Data. New York:Wiley.

Dong, G., and R. Harris. 2015. “Spatial Autorgressive Models for Geographically Hierarchical Data Structures.” Geographical Analysis 47 (2): 173–91.

Dong, G., R. Harris, K. Jones, and J. Yu. 2015. “Multilevel Modeling with Spatial Interaction Effects with Application to an Emerging Land Market in Beijing, China.” PLOS One 10 (6). https://doi.org/doi:10.1371/journal.pone.0130761.

Dong, Guanpeng, Richard Harris, and Angelos Mimis. 2020. HSAR: Hierarchical Spatial Autoregressive Model. https://CRAN.R-project.org/package=HSAR.

Gaetan, Carlo, and Xavier Guyon. 2010. Spatial Statistics and Modeling. New York: Springer.

Gerber, Florian, and Reinhard Furrer. 2015. “Pitfalls in the Implementation of Bayesian Hierarchical Modeling of Areal Count Data: An Illustration Using Bym and Leroux Models.” Journal of Statistical Software, Code Snippets 63 (1): 1–32. https://doi.org/10.18637/jss.v063.c01.

Gómez-Rubio, V. 2019. “Spatial Data Analysis with Inla. Coding Club Uc3m Tutorial Series. Universidad Carlos Iii de Madrid.” https://codingclubuc3m.rbind.io/talk/2019-11-05/.

Gómez-Rubio, V. 2020. Bayesian Inference with Inla. Boca Raton, FL: CRC Press.

Hepple, Leslie W. 1976. “A Maximum Likelihood Model for Econometric Estimation with Spatial Series.” In Theory and Practice in Regional Science, edited by I. Masser, 90–104. London Papers in Regional Science. London: Pion.

Lieshout, M. N. M. van. 2019. Theory of Spatial Statistics. Boca Raton, FL: Chapman; Hall/CRC.

McCulloch, Charles E., and Shayle R. Searle. 2001. Generalized, Linear, and Mixed Models. New York: Wiley.

Ord, J. K. 1975. “Estimation Methods for Models of Spatial Interaction.” Journal of the American Statistical Association 70 (349): 120–26.

Pinheiro, Jose C., and Douglas M. Bates. 2000. Mixed-Effects Models in S and S-Plus. New York: Springer.

Ripley, B. D. 1981. Spatial Statistics. New York: Wiley.

Ripley, B. 1988. Statistical Inference for Spatial Processes. Cambridge: Cambridge University Press.

Rue, Havard, Finn Lindgren, Daniel Simpson, Sara Martino, Elias Teixeira Krainski, Haakon Bakka, Andrea Riebler, and Geir-Arne Fuglstad. 2021. INLA: Full Bayesian Analysis of Latent Gaussian Models Using Integrated Nested Laplace Approximations.

Umlauf, Nikolaus, Daniel Adler, Thomas Kneib, Stefan Lang, and Achim Zeileis. 2015. “Structured Additive Regression Models: An R Interface to BayesX.” Journal of Statistical Software 63 (21): 1–46. http://www.jstatsoft.org/v63/i21/.

Umlauf, Nikolaus, Thomas Kneib, Stefan Lang, and Achim Zeileis. 2021. R2BayesX: Estimate Structured Additive Regression Models with Bayesx. https://CRAN.R-project.org/package=R2BayesX.

Vranckx, M., T. Neyens, and C. Faes. 2019. “Comparison of Different Software Implementations for Spatial Disease Mapping.” Spatial and Spatio-Temporal Epidemiology 31: 100302. https://doi.org/10.1016/j.sste.2019.100302.

Wall, M. M. 2004. “A Close Look at the Spatial Structure Implied by the CAR and SAR Models.” Journal of Statistical Planning and Inference 121: 311–24.

Whittle, P. 1954. “On Stationary Processes in the Plane.” Biometrika 41 (3-4): 434–49. https://doi.org/10.1093/biomet/41.3-4.434.

Wood, Simon. 2022. Mgcv: Mixed Gam Computation Vehicle with Automatic Smoothness Estimation. https://CRAN.R-project.org/package=mgcv.

Wood, S.N. 2017. Generalized Additive Models: An Introduction with R. 2nd ed. Chapman; Hall/CRC.