Spatial simultaneous autoregressive terms: sar()

sar(M, type = ) in a formula makes the residuals of a gaussian() or student() response one spatially autoregressive draw over a spatial weights matrix M, with brms's spelling and its two types. It is available once frmtmb.spatial is loaded: library(frmtmb.spatial). Without the package, frmtmb::frm() refuses the term and names this package.

Usage

sar(M, type = "lag")

Arguments

M: The spatial weights matrix; see "Arguments of the term".
type: "lag" or "error".

Value

sar() returns the term as an object of class frmtmb_ac_term. Its use is inside a formula, where frmtmb::frm() reads it.

Details

With \(A = I - \rho W\), the two types are, as in brms 2.23.0:

type = "lag" (the default): SAR of the response, \(y = \rho W y + \eta + e\), so \(A y - \eta = e\). The mean of the response is \(A^{-1} \eta\), which fitted() reports, as brms's posterior_epred() does. The parameter is lagsar.
type = "error": SAR of the residuals, \(y = \eta + u\), \(u = \rho W u + e\), so \(A (y - \eta) = e\). The mean is \(\eta\). The parameter is errorsar.

e holds independent normal residuals with standard deviation sigma, or for student() one multivariate-t draw with scale sigma and nu. The log-likelihood has the log-determinant \(\log|I - \rho W| = \sum_i \log(1 - \rho \lambda_i)\), with \(\lambda_i\) the eigenvalues of W, computed once.

Arguments of the term

M: The spatial weights matrix, one row and column per row of the data, in data2 = list(M = M): a matrix, a sparse Matrix, or a listw or nb object of the spdep package, which is read as brms reads it (spdep::listw2mat(), and spdep::nb2mat() with its row-standardized default) without needing spdep. The rows the model drops (missing values) are dropped from M, as brms drops them. M is not checked for row standardization or symmetry, as in brms. Its eigenvalues must be real, with at least one negative and one positive.
type: "lag" or "error". A unique abbreviation is read as brms reads it (match.arg()).

The range of rho

\(I - \rho W\) is invertible and its determinant positive on the interval \((1 / \lambda_{min}, 1 / \lambda_{max})\), the default interval of spatialreg::lagsarlm(). frmtmb estimates \(\rho\) on that interval through a logistic map, starting at 0. brms 2.23.0 bounds \(\rho\) by \((\lambda_{min}, \lambda_{max})\) instead. For a row-standardized W (\(\lambda_{max} = 1\)) that is a narrower interval with the same upper end, which cuts off the negative values below \(\lambda_{min}\). For a W that is not row-standardized, brms's interval reaches past \(1 / \lambda_{max}\), where the determinant changes sign; there brms's log-determinant is the log of a negative number and Stan rejects the proposal, so its support is frmtmb's interval there.

The consequence on a row-standardized W: when the data favor rho below the smallest eigenvalue, frmtmb and brms fit different models. On a row-standardized 8 x 8 queen lattice (eigenvalues -0.5006 to 1) with data drawn at rho = -0.8, frmtmb's estimate is -1.103, with a log-likelihood 3.74 above the best inside brms's bounds, and brms's posterior piles at its bound, -0.5006 (median -0.443; a fifth of the draws within 0.02 of it).

A W whose extreme eigenvalue is repeated in all but one direction, as on a complete graph, gives a likelihood with no maximum at that end when the model has an intercept; frmtmb warns.

What else the term does

Only gaussian() and student() take it, with a constant sigma and nu, as in brms ("SAR models are not implemented when predicting 'sigma' or 'nu'"). It is written in the formula of mu. The likelihood is one joint density over all the rows, so weights(), cens(), trunc(), se(), mi(), rescor = TRUE, mixtures and quadrature are refused, and a response carries one residual correlation term. Random effects work.

summary() and frmtmb::confint_varcorr() report lagsar or errorsar with a delta-method interval. fitted() and residuals() use the mean above; pearson residuals divide by the marginal standard deviation, \(\sigma \sqrt{[(A'A)^{-1}]_{ii}}\). simulate() and predict() draw the whole response through the model. frmtmb::autocor_matrix() returns \((A'A)^{-1}\). fitted(), predict() and simulate() take newdata with one row per row of M, read as the locations in their order, as brms takes it: a counterfactual such as x + 1 everywhere gives the lag model's spillover, \(A^{-1} \eta_{new}\). When the fit dropped rows for missing values, newdata covers them too, and A is then over every row of M, at the fitted rho; a missing predictor gives NA, which the lag form spreads to every row, as in brms. Other newdata is refused; brms refuses it too unless a matrix over the new rows comes in its newdata2. frmtmb::frm_linpred() gives the linear predictor \(\eta\) on any new data.

brms's class priors lagsar and errorsar work: set_prior("normal(0, 0.5)", class = "lagsar") is a density on rho, carried to its internal parameter with the log-Jacobian of the map. frmtmb.sample's default is brms's, flat on rho over its interval. frmtmb.sample::log_lik() gives brms's pointwise log-density: the leave-one-out conditional density of each row given the other rows, normal (or Student-t with \(\nu + n - 1\) degrees of freedom) with precision \(A'A / \sigma^2\), which is sparse, so its cost per draw is the nonzeros of W and not the dense n x n inverse brms takes. The row sum is not the log-likelihood.

Examples

set.seed(2)
grid <- expand.grid(r = 1:6, c = 1:6)
n <- nrow(grid)
A <- outer(seq_len(n), seq_len(n), function(i, j) {
  as.numeric(abs(grid$r[i] - grid$r[j]) + abs(grid$c[i] - grid$c[j]) == 1)
})
W <- A / rowSums(A)
d <- data.frame(x = rnorm(n))
d$y <- as.vector(solve(diag(n) - 0.5 * W, 1 + 0.5 * d$x + rnorm(n)))
fit <- frm(bf(y ~ x + sar(W, type = "lag")), family = gaussian(),
           data = d, data2 = list(W = W))
summary(fit)
#>  Family: gaussian 
#>  Links: mu = identity; sigma = log
#> 
#> Formula: y ~ x + sar(W, type = "lag") 
#>    Data: d (Number of observations: 36) 
#>  Method: ML   logLik: -56.5871   AIC: 121.174   BIC: 127.508 
#> 
#> Correlation Structures:
#>        Estimate Est.Error l-95% CI u-95% CI
#> lagsar      0.7      0.13     0.44     0.96
#> 
#> Regression Coefficients:
#>           Estimate Est.Error l-95% CI u-95% CI z value Pr(>|z|)
#> Intercept     0.47      0.32    -0.15     1.10    1.49     0.14
#> x             0.26      0.16    -0.06     0.59    1.61     0.11
#> 
#> Further Distributional Parameters:
#>       Estimate Est.Error l-95% CI u-95% CI
#> sigma     1.06      0.13     0.83     1.36

Spatial simultaneous autoregressive terms: `sar()`