Fit a spatial regression model using eigenvector spatial filtering (ESF).

```
stan_esf(
formula,
slx,
re,
data,
C,
EV = make_EV(C, nsa = nsa, threshold = threshold),
nsa = FALSE,
threshold = 0.25,
family = gaussian(),
prior = NULL,
ME = NULL,
centerx = FALSE,
censor_point,
prior_only = FALSE,
chains = 4,
iter = 2000,
refresh = 500,
keep_all = FALSE,
slim = FALSE,
drop = NULL,
pars = NULL,
control = NULL,
...
)
```

Chun, Y., D. A. Griffith, M. Lee and P. Sinha (2016). Eigenvector selection with stepwise regression techniques to construct eigenvector spatial filters. *Journal of Geographical Systems*, 18(1), 67-85. doi:10.1007/s10109-015-0225-3
.

Dray, S., P. Legendre & P. R. Peres-Neto (2006). Spatial modelling: a comprehensive framework for principal coordinate analysis of neighbour matrices (PCNM). *Ecological Modeling*, 196(3-4), 483-493.

Donegan, C., Y. Chun and A. E. Hughes (2020). Bayesian estimation of spatial filters with Moran’s Eigenvectors and hierarchical shrinkage priors. *Spatial Statistics*. doi:10.1016/j.spasta.2020.100450
(open access: doi:10.31219/osf.io/fah3z
).

Donegan, Connor (2021). Building spatial conditional autoregressive (CAR) models in the Stan programming language. *OSF Preprints*. doi:10.31219/osf.io/3ey65
.

Griffith, Daniel A., and P. R. Peres-Neto (2006). Spatial modeling in ecology: the flexibility of eigenfunction spatial analyses. *Ecology* 87(10), 2603-2613.

Griffith, D., and Y. Chun (2014). Spatial autocorrelation and spatial filtering, Handbook of Regional Science. Fischer, MM and Nijkamp, P. eds.

Griffith, D., Chun, Y. and Li, B. (2019). *Spatial Regression Analysis Using Eigenvector Spatial Filtering*. Elsevier.

Piironen, J and A. Vehtari (2017). Sparsity information and regularization in the horseshoe and other shrinkage priors. In *Electronic Journal of Statistics*, 11(2):5018-5051.

- formula
A model formula, following the R

`formula`

syntax. Binomial models are specified by setting the left hand side of the equation to a data frame of successes and failures, as in`cbind(successes, failures) ~ x`

.- slx
Formula to specify any spatially-lagged covariates. As in,

`~ x1 + x2`

(the intercept term will be removed internally). When setting priors for`beta`

, remember to include priors for any SLX terms.- re
To include a varying intercept (or "random effects") term,

`alpha_re`

, specify the grouping variable here using formula syntax, as in`~ ID`

. Then,`alpha_re`

is a vector of parameters added to the linear predictor of the model, and:`alpha_re ~ N(0, alpha_tau) alpha_tau ~ Student_t(d.f., location, scale).`

- data
A

`data.frame`

or an object coercible to a data frame by`as.data.frame`

containing the model data.- C
Spatial connectivity matrix which will be used to calculate eigenvectors, if

`EV`

is not provided by the user. Typically, the binary connectivity matrix is best for calculating eigenvectors (i.e., using`C = shape2mat(shape, style = "B")`

). This matrix will also be used to calculate residual spatial autocorrelation and any user specified`slx`

terms; it will be row-standardized before calculating`slx`

terms. See`shape2mat`

.- EV
A matrix of eigenvectors from any (transformed) connectivity matrix, presumably spatial (see

`make_EV`

). If`EV`

is provided, still also provide a spatial weights matrix`C`

for other purposes;`threshold`

and`nsa`

are ignored for user provided`EV`

.- nsa
Include eigenvectors representing negative spatial autocorrelation? Defaults to

`nsa = FALSE`

. This is ignored if`EV`

is provided.- threshold
Eigenvectors with standardized Moran coefficient values below this

`threshold`

value will be excluded from the candidate set of eigenvectors,`EV`

. This defaults to`threshold = 0.25`

, and is ignored if`EV`

is provided.- family
The likelihood function for the outcome variable. Current options are

`family = gaussian()`

,`student_t()`

and`poisson(link = "log")`

, and`binomial(link = "logit")`

.- prior
A named list of parameters for prior distributions (see

`priors`

):- intercept
The intercept is assigned a Gaussian prior distribution (see

`normal`

.
- beta
Regression coefficients are assigned Gaussian prior distributions. Variables must follow their order of appearance in the model

`formula`

. Note that if you also use`slx`

terms (spatially lagged covariates), and you use custom priors for`beta`

, then you have to provide priors for the slx terms. Since slx terms are*prepended*to the design matrix, the prior for the slx term will be listed first.- sigma
For

`family = gaussian()`

and`family = student_t()`

models, the scale parameter,`sigma`

, is assigned a (half-) Student's t prior distribution. The half-Student's t prior for`sigma`

is constrained to be positive.- nu
`nu`

is the degrees of freedom parameter in the Student's t likelihood (only used when`family = student_t()`

).`nu`

is assigned a gamma prior distribution. The default prior is`prior = list(nu = gamma(alpha = 3, beta = 0.2))`

.- tau
The scale parameter for random effects, or varying intercepts, terms. This scale parameter,

`tau`

, is assigned a half-Student's t prior. To set this, use, e.g.,`prior = list(tau = student_t(df = 20, location = 0, scale = 20))`

.- beta_ev
The eigenvector coefficients are assigned the horseshoe prior (Piironen and Vehtari, 2017), parameterized by

`global_scale`

(to control overall prior sparsity), plus the degrees of freedom and scale of a Student's t model for any large coefficients (see`priors`

). To allow the spatial filter to account for a greater amount of spatial autocorrelation (i.e., if you find the residuals contain spatial autocorrelation), increase the global scale parameter (to a maximum of`global_scale = 1`

).

- ME
To model observational uncertainty (i.e. measurement or sampling error) in any or all of the covariates, provide a list of data as constructed by the

`prep_me_data`

function.- centerx
To center predictors on their mean values, use

`centerx = TRUE`

. If the ME argument is used, the modeled covariate (i.e., latent variable), rather than the raw observations, will be centered. When using the ME argument, this is the recommended method for centering the covariates.- censor_point
Integer value indicating the maximum censored value; this argument is for modeling censored (suppressed) outcome data, typically disease case counts or deaths. For example, the US Centers for Disease Control and Prevention censors (does not report) death counts that are nine or fewer, so if you're using CDC WONDER mortality data you could provide

`censor_point = 9`

.- prior_only
Draw samples from the prior distributions of parameters only.

- chains
Number of MCMC chains to estimate. Default

`chains = 4`

.- iter
Number of samples per chain. Default

`iter = 2000`

.- refresh
Stan will print the progress of the sampler every

`refresh`

number of samples. Defaults to`500`

; set`refresh=0`

to silence this.- keep_all
If

`keep_all = TRUE`

then samples for all parameters in the Stan model will be kept; this is necessary if you want to do model comparison with Bayes factors and the`bridgesampling`

package.- slim
If

`slim = TRUE`

, then the Stan model will not collect the most memory-intensive parameters (including n-length vectors of fitted values, log-likelihoods, and ME-modeled covariate values). This will disable many convenience functions that are otherwise available for fitted`geostan`

models, such as the extraction of residuals, fitted values, and spatial trends, WAIC, and spatial diagnostics, and ME diagnostics; many quantities of interest, such as fitted values and spatial trends, can still be calculated manually using given parameter estimates. The "slim" option is useful for data-intensive routines, such as regression with raster data, Monte Carlo studies, and measurement error models. For more control over which parameters are kept or dropped, use the`drop`

argument instead of`slim`

.- drop
Provide a vector of character strings to specify the names of any parameters that you do not want MCMC samples for. Dropping parameters in this way can improve sampling speed and reduce memory usage. The following parameter vectors can potentially be dropped from ESF models:

- fitted
The N-length vector of fitted values

- log_lik
The N-length vector of pointwise log-likelihoods, which is used to calculate WAIC.

- alpha_re
Vector of 'random effects'/varying intercepts.

- x_true
N-length vector of 'latent'/modeled covariate values created for measurement error (ME) models.

- esf
The N-length eigenvector spatial filter.

- beta_ev
The vector of coefficients for the eigenvectors.

Using

`drop = c('fitted', 'log_lik', 'alpha_re', 'x_true', 'esf', 'beta_ev')`

is equivalent to`slim = TRUE`

. Note that if`slim = TRUE`

, then`drop`

will be ignored---so only use one or the other.- pars
Optional; specify any additional parameters you'd like stored from the Stan model.

- control
A named list of parameters to control the sampler's behavior. See stan for details.

- ...
Other arguments passed to sampling.

An object of class class `geostan_fit`

(a list) containing:

- summary
Summaries of the main parameters of interest; a data frame

- diagnostic
Widely Applicable Information Criteria (WAIC) with a measure of effective number of parameters (

`eff_pars`

) and mean log pointwise predictive density (`lpd`

), and mean residual spatial autocorrelation as measured by the Moran coefficient.- data
a data frame containing the model data

- EV
A matrix of eigenvectors created with

`w`

and`geostan::make_EV`

- C
The spatial weights matrix used to construct EV

- family
the user-provided or default

`family`

argument used to fit the model- formula
The model formula provided by the user (not including ESF component)

- slx
The

`slx`

formula- re
A list containing

`re`

, the random effects (varying intercepts) formula if provided, and`data`

a data frame with columns`id`

, the grouping variable, and`idx`

, the index values assigned to each group.- priors
Prior specifications.

- x_center
If covariates are centered internally (

`centerx = TRUE`

), then`x_center`

is a numeric vector of the values on which covariates were centered.- ME
The

`ME`

data list, if one was provided by the user for measurement error models.- spatial
A data frame with the name of the spatial component parameter ("esf") and method ("ESF")

- stanfit
an object of class

`stanfit`

returned by`rstan::stan`

Eigenvector spatial filtering (ESF) is a method for spatial regression analysis. ESF is extensively covered in Griffith et al. (2019). This function implements the methodology introduced in Donegan et al. (2020), which uses Piironen and Vehtari's (2017) regularized horseshoe prior.

ESF decomposes spatial autocorrelation into a linear combination of various patterns, typically at different scales (such as local, regional, and global trends). By adding a spatial filter to a regression model, these spatial autocorrelation patterns are shifted from the residuals to the spatial filter. ESF models take the spectral decomposition of a transformed spatial connectivity matrix, \(C\). The resulting eigenvectors, \(E\), are mutually orthogonal and uncorrelated map patterns. The spatial filter equals \(E \beta_{E}\) where \(\beta_{E}\) is a vector of coefficients.

ESF decomposes the data into a global mean, \(\alpha\), global patterns contributed by covariates \(X \beta\), spatial trends \(E \beta_{E}\), and residual variation. Thus, for `family=gaussian()`

,

$$ y \sim Gauss(\alpha + X * \beta + E \beta_{E}, \sigma). $$

An ESF component can be incorporated into the linear predictor of any generalized linear model. For example, a spatial Poisson model for rare disease incidence may be specified as follows:

$$y \sim Poisson(e^{O + \mu})$$ $$\mu = \alpha + E \beta_{E} + A $$ $$ A \sim Guass(0, \tau) $$ $$ \tau \sim student(20, 0, 2) $$ $$ \beta_{E} \sim horseshoe(.) $$

The form of this model is similar to the BYM model (see stan_icar), in the sense that it contains a spatially structured trend term (\(E \beta_{E}\)) and an unstructured ('random effects') term (\(A\)).

The spatial.geostan_fit method will return \(E \beta_{E}\).

The model can also be extended to the space-time domain; see shape2mat to specify a space-time connectivity matrix.

The coefficients \(\beta_{E}\) are assigned the regularized horseshoe prior (Piironen and Vehtari, 2017), resulting in a relatively sparse model specification. In addition, numerous eigenvectors are automatically dropped because they represent trace amounts of spatial autocorrelation (this is controlled by the `threshold`

argument). By default, `stan_esf`

will drop all eigenvectors representing negative spatial autocorrelation patterns. You can change this behavior using the `nsa`

argument.

The CAR models can also incorporate spatially-lagged covariates, measurement/sampling error in covariates (particularly when using small area survey estimates as covariates), and censored outcomes (such as arise when a disease surveillance system suppresses data for privacy reasons). For details on these options, please see the Details section in the documentation for stan_glm.

```
# \donttest{
data(sentencing)
# spatial weights matrix with binary coding scheme
C <- shape2mat(sentencing, style = "B")
# log-expected number of sentences
## expected counts are based on county racial composition and mean sentencing rates
log_e <- log(sentencing$expected_sents)
# fit spatial Poisson model with ESF + unstructured 'random effects'
fit.esf <- stan_esf(sents ~ offset(log_e),
re = ~ name,
family = poisson(),
data = sentencing,
C = C,
chains = 2, iter = 800) # for speed only
# spatial diagnostics
sp_diag(fit.esf, sentencing)
plot(fit.esf)
# plot marginal posterior distributions of beta_ev (eigenvector coefficients)
plot(fit.esf, pars = "beta_ev")
# plot the marginal posterior distributions of the spatial filter
plot(fit.esf, pars = "esf")
# calculate log-standardized incidence ratios
library(ggplot2)
library(sf)
f <- fitted(fit.esf, rates = FALSE)$mean
SSR <- f / sentencing$expected_sents
log.SSR <- log( SSR, base = 2 )
# map the log-SSRs
st_as_sf(sentencing) %>%
ggplot() +
geom_sf(aes(fill = log.SSR)) +
scale_fill_gradient2(
midpoint = 0,
name = NULL,
breaks = seq(-3, 3, by = 0.5)
) +
labs(title = "Log-Standardized Sentencing Ratios",
subtitle = "log( Fitted/Expected ), base 2"
) +
theme_void()
# }
```