This vignette provides a tutorial for fitting spatial regression models to raster data using **geostan**. The term “raster” is used here to refer to any regularly spaced set of observations such that the data can be represented spatially by a rectangular grid. Remotely sensed imagery is a common form of raster data.

**geostan** can be used for spatial regression with fairly large raster data layers, although the functionality of these models will often be limited to the estimation of regression coefficients and spatial autocorrelation parameters. Limited experience thus far finds that **geostan**’s spatial autoregressive models can be fit to raster layers with two hundred thousand observations using a laptop computer and fewer than ten minutes of sampling time.

Start by loading some necessary R packages.

```
library(geostan)
library(sf)
set.seed(1127)
```

We will create a small raster data layer for the purpose of illustration.

```
40
row <- 30
col <-c(N <- row * col)
```

`## [1] 1200`

```
st_sfc(st_polygon(list(rbind(c(0,0), c(col,0), c(col,row), c(0,0)))))
sfc = st_make_grid(sfc, cellsize = 1, square = TRUE)
grid <- st_as_sf(grid)
grid <- shape2mat(grid, style = "W", queen = FALSE)
W <-$z <- sim_sar(w = W, rho = 0.9)
grid$y <- -0.5 * grid$z + sim_sar(w = W, rho = .9, sigma = .3) grid
```

`plot(grid[,'z'])`

The following R code will fit a spatial autoregressive model to these data:

` stan_sar(y ~ z, data = grid, C = W) fit <-`

The `stan_sar`

function will take the spatial weights matrix `W`

and pass it through a function called `prep_sar_data`

which will calculate the eigenvalues of the spatial weights matrix using `base::eigen`

, as required for computational reasons. This step is prohibitive for large data sets (e.g., \(N = 100,000\)).

The following code would normally be used to fit a conditional autoregressive (CAR) model:

```
shape2mat(grid, style = "B", queen = FALSE)
C <- prep_car_data(C, "WCAR")
car_list <- stan_car(y ~ z, data = grid, car_parts = car_list) fit <-
```

Here, the `prep_car_data`

function calculates the eigenvalues of the spatial weights matrix using `base::eigen`

, which is not feasible for large N.

The `prep_sar_data2`

and `prep_car_data2`

functions are designed for raster layers. As input, they require the dimensions of the grid (number of rows and number of columns). The eigenvalues are produced very quickly using Equation 5 from Griffith (2000). The methods have certain restrictions. First, this is only applicable to raster layers—regularly spaced, rectangular grids of observations. Second, to define which observations are adjacent to one another, the “rook” criteria is used (spatially, only observations that share an edge are defined as neighbors to one another). Third, the spatial adjacency matrix will be row-standardized. This is standard (and required) for SAR models, and it corresponds to the “WCAR” specification of the CAR model (see Donegan 2022).

The following code will fit a SAR model to our `grid`

data, and is suitable for much larger raster layers:

` prep_sar_data2(row = row, col = col) sar_list <-`

`## Range of permissible rho values: -1 1`

```
stan_sar(y ~ z,
fit <-data = grid,
centerx = TRUE,
sar_parts = sar_list,
iter = 500,
chains = 4,
slim = TRUE #,
# cores = 4, # for multi-core processing
)
```

```
##
## *Setting prior parameters for intercept
```

`## Distribution: normal`

```
## location scale
## 1 0.035 5
```

```
##
## *Setting prior parameters for beta
## Distribution: normal
```

```
## location scale
## 1 0 5
```

```
##
## *Setting prior for SAR scale parameter (sar_scale)
```

`## Distribution: student_t`

```
## df location scale
## 1 10 0 3
```

```
##
## *Setting prior for SAR spatial autocorrelation parameter (sar_rho)
```

`## Distribution: uniform`

```
## lower upper
## 1 -1 1
##
## SAMPLING FOR MODEL 'foundation' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.000508 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 5.08 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 500 [ 0%] (Warmup)
## Chain 1: Iteration: 251 / 500 [ 50%] (Sampling)
## Chain 1: Iteration: 500 / 500 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 3.553 seconds (Warm-up)
## Chain 1: 0.77 seconds (Sampling)
## Chain 1: 4.323 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'foundation' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0.000524 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 5.24 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 500 [ 0%] (Warmup)
## Chain 2: Iteration: 251 / 500 [ 50%] (Sampling)
## Chain 2: Iteration: 500 / 500 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 1.489 seconds (Warm-up)
## Chain 2: 0.806 seconds (Sampling)
## Chain 2: 2.295 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'foundation' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0.000443 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 4.43 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 500 [ 0%] (Warmup)
## Chain 3: Iteration: 251 / 500 [ 50%] (Sampling)
## Chain 3: Iteration: 500 / 500 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 1.436 seconds (Warm-up)
## Chain 3: 0.773 seconds (Sampling)
## Chain 3: 2.209 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'foundation' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0.000763 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 7.63 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 500 [ 0%] (Warmup)
## Chain 4: Iteration: 251 / 500 [ 50%] (Sampling)
## Chain 4: Iteration: 500 / 500 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 1.446 seconds (Warm-up)
## Chain 4: 0.691 seconds (Sampling)
## Chain 4: 2.137 seconds (Total)
## Chain 4:
```

`print(fit) `

```
## Spatial Model Results
## Formula: y ~ z
## Spatial method (outcome): SAR
## Likelihood function: auto_gaussian
## Link function: identity
## Residual Moran Coefficient: NA
## Observations: 1200
## Data models (ME): none
## Inference for Stan model: foundation.
## 4 chains, each with iter=500; warmup=250; thin=1;
## post-warmup draws per chain=250, total post-warmup draws=1000.
##
## mean se_mean sd 2.5% 20% 50% 80% 97.5% n_eff Rhat
## intercept 0.035 0.002 0.073 -0.104 -0.026 0.035 0.099 0.168 1009 1.000
## z -0.510 0.000 0.009 -0.527 -0.517 -0.510 -0.501 -0.493 1350 0.998
## sar_rho 0.874 0.000 0.015 0.844 0.862 0.874 0.886 0.902 924 1.004
## sar_scale 0.312 0.000 0.006 0.299 0.306 0.311 0.317 0.324 1026 1.001
##
## Samples were drawn using NUTS(diag_e) at Tue Apr 16 08:24:30 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 user first creates the data list using `prep_sar_data2`

and then passes it to `stan_sar`

using the `sar_parts`

argument. Also, `slim = TRUE`

is invoked to prevent the model from collecting N-length parameter vectors and quantities of interest (such as fitted values and log-likelihoods).

For large data sets and complex models, `slim = TRUE`

can bring about computational improvements at the cost of losing some functionality (including the loss of convenience functions like `sp_diag`

, `me_diag`

, `spatial`

, `resid`

, and `fitted`

). Many quantities of interest, such as fitted values and spatial trend terms, can still be calculated manually using the data and parameter estimates (intercept, coefficients, and spatial autocorrelation parameters).

The favorable MCMC diagnostics for this model (sufficiently large effective sample sizes `n_eff`

, and `Rhat`

values very near to 1), based on just 250 post-warmup iterations per chain with four MCMC chains, provides some indication as to how computationally efficient these spatial autoregressive models can be.

Also, note that Stan usually samples more efficiently when variables have been mean-centered. Using the `centerx = TRUE`

argument in `stan_sar`

(or any other model-fitting function in **geostan**) can be very helpful in this respect. Also note that the SAR models in **geostan** are (generally) no less computationally-efficient than the CAR models, and may even be slightly more efficient.

Donegan, Connor. 2022. “Building Spatial Conditional Autoregressive Models in the Stan Programming Language.” *OSF Preprints*. https://doi.org/10.31219/osf.io/3ey65.

Griffith, Daniel A. 2000. “Eigenfunction Properties and Approximations of Selected Incidence Matrices Employed in Spatial Analyses.” *Linear Algebra and Its Applications* 321 (1-3): 95–112.