stan_car {geostan} | R Documentation |
Description
Use the CAR model as a prior on parameters, or fit data to a spatial Gaussian CAR model.
Usage
stan_car( formula, slx, re, data, car_parts, C, family = gaussian(), prior = NULL, ME = NULL, centerx = FALSE, prior_only = FALSE, censor_point, chains = 4, iter = 2000, refresh = 500, keep_all = FALSE, slim = FALSE, drop = NULL, pars = NULL, control = NULL, ...)
Arguments
formula | A model formula, following the R |
slx | Formula to specify any spatially-lagged covariates. As in, |
re | To include a varying intercept (or "random effects") term, alpha_re ~ N(0, alpha_tau)alpha_tau ~ Student_t(d.f., location, scale). With the CAR model, any |
data | A |
car_parts | A list of data for the CAR model, as returned by |
C | Optional spatial connectivity matrix which will be used to calculate residual spatial autocorrelation as well as any user specified |
family | The likelihood function for the outcome variable. Current options are |
prior | A named list of parameters for prior distributions (see
.
|
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 |
centerx | To center predictors on their mean values, use |
prior_only | Logical value; if |
censor_point | Integer value indicating the maximum censored value; this argument is for modeling censored (suppressed) outcome data, typically disease case counts or deaths. |
chains | Number of MCMC chains to use. |
iter | Number of samples per chain. |
refresh | Stan will print the progress of the sampler every |
keep_all | If |
slim | If |
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 CAR models:
If |
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 |
... | Other arguments passed to |
Details
CAR models are discussed in Cressie and Wikle (2011, p. 184-88), Cressie (2015, Ch. 6-7), and Haining and Li (2020, p. 249-51). It is often used for areal or lattice data.
Details for the Stan code for this implementation of the CAR model can be found in Donegan (2021).
The general scheme for the CAR model is as follows:
y \sim Gauss( \mu, ( I - \rho C)^{-1} M),
where I
is the identity matrix, \rho
is a spatial dependence parameter, C
is a spatial connectivity matrix, and M
is a diagonal matrix of variance terms. The diagonal of M
contains a scale parameter \tau
multiplied by a vector of weights (often set to be proportional to the inverse of the number of neighbors assigned to each site). The CAR model owes its name to the fact that this joint distribution corresponds to a set of conditional distributions that relate the expected value of each observation to a function of neighboring values, i.e., the Markov condition holds:
E(y_i | y_1, y_2, \dots, y_{i-1}, y_{i+1}, \dots, y_n) = \mu_i + \rho \sum_{j=1}^n c_{i,j} (y_j - \mu_j),
where entries of c_{i,j}
are non-zero only if j \in N(i)
and N(i)
indexes the sites that are neighbors of the i^{th}
site.
With the Gaussian probability distribution,
y_i | y_j: j \neq i \sim Gauss(\mu_i + \rho \sum_{j=1}^n c_{i,j} (y_j - \mu_j), \tau_i^2)
where \tau_i
is a scale parameter and \mu_i
may contain covariates or simply the intercept.
The covariance matrix of the CAR model contains two parameters: \rho
(car_rho
) which controls the kind (positive or negative) and degree of spatial autocorrelation, and the scale parameter \tau
(car_scale
). The range of permissible values for \rho
depends on the specification of \boldsymbol C
and \boldsymbol M
; for specification options, see prep_car_data and Cressie and Wikle (2011, pp. 184-188) or Donegan (2021).
Further details of the models and results depend on the family
argument, as well as on the particular CAR specification chosen (from prep_car_data).
Auto-Gaussian
When family = auto_gaussian()
(the default), the CAR model is applied directly to the data as follows:
y \sim Gauss( \mu, (I - \rho C)^{-1} M),
where \mu
is the mean vector (with intercept, covariates, etc.), C
is a spatial connectivity matrix, and M
is a known diagonal matrix containing the conditional variances \tau_i^2
. C
and M
are provided by prep_car_data.
The auto-Gaussian model contains an implicit spatial trend (i.e. autocorrelation) component \phi
which can be calculated as follows (Cressie 2015, p. 564):
\phi = \rho C (y - \mu).
This term can be extracted from a fitted auto-Gaussian model using the spatial method.
When applied to a fitted auto-Gaussian model, the residuals.geostan_fit method returns 'de-trended' residuals R
by default. That is,
R = y - \mu - \rho C (y - \mu).
To obtain "raw" residuals (y - \mu
), use residuals(fit, detrend = FALSE)
. Similarly, the fitted values obtained from the fitted.geostan_fit will include the spatial trend term by default.
Poisson
For family = poisson()
, the model is specified as:
y \sim Poisson(e^{O + \lambda})
\lambda \sim Gauss(\mu, (I - \rho C)^{-1} \boldsymbol M).
If the raw outcome consists of a rate \frac{y}{p}
with observed counts y
and denominator p
(often this will be the size of the population at risk), then the offset term O=log(p)
is the log of the denominator.
This is often written (equivalently) as:
y \sim Poisson(e^{O + \mu + \phi})
\phi \sim Gauss(0, (I - \rho C)^{-1} \boldsymbol M).
For Poisson models, the spatial method returns the parameter vector \phi
.
In the Poisson CAR model, \phi
contains a latent spatial trend as well as additional variation around it: \phi_i = \rho \sum_{i=1}^n c_{ij} \phi_j + \epsilon_i
, where \epsilon_i \sim Gauss(0, \tau_i^2)
. If you would like to extract the latent/implicit spatial trend from \phi
, you can do so by calculating (following Cressie 2015, p. 564):
\rho C \phi.
Binomial
For family = binomial()
, the model is specified as:
y \sim Binomial(N, \lambda)
logit(\lambda) \sim Gauss(\mu, (I - \rho C)^{-1} \boldsymbol M).
where outcome data y
are counts, N
is the number of trials, \lambda
is the 'success' rate, and \mu
contains the intercept and possibly covariates. Note that the model formula should be structured as: cbind(sucesses, failures) ~ x
, such that trials = successes + failures
.
This is often written (equivalently) as:
y \sim Binomial(N, \lambda)
logit(\lambda) = \mu + \phi
\phi \sim Gauss(0, (I - \rho C)^{-1} \boldsymbol M).
For fitted Binomial models, the spatial method will return the parameter vector phi
.
As is also the case for the Poisson model, \phi
contains a latent spatial trend as well as additional variation around it. If you would like to extract the latent/implicit spatial trend from \phi
, you can do so by calculating:
\rho C \phi.
Additional functionality
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.
Value
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.- stanfit
an object of class
stanfit
returned byrstan::stan
- data
a data frame containing the model data
- family
the user-provided or default
family
argument used to fit the model- formula
The model formula provided by the user (not including CAR component)
- slx
The
slx
formula- re
A list containing
re
, the varying intercepts (re
) formula if provided, andData
a data frame with columnsid
, the grouping variable, andidx
, the index values assigned to each group.- priors
Prior specifications.
- x_center
If covariates are centered internally (
centerx = TRUE
), thenx_center
is a numeric vector of the values on which covariates were centered.- spatial
A data frame with the name of the spatial component parameter (either "phi" or, for auto Gaussian models, "trend") and method ("CAR")
- ME
A list indicating if the object contains an ME model; if so, the user-provided ME list is also stored here.
- C
Spatial connectivity matrix (in sparse matrix format).
Author(s)
Connor Donegan, connor.donegan@gmail.com
Source
Besag, Julian (1974). Spatial interaction and the statistical analysis of lattice systems. Journal of the Royal Statistical Society B36.2: 192–225.
Cressie, Noel (2015 (1993)). Statistics for Spatial Data. Wiley Classics, Revised Edition.
Cressie, Noel and Wikle, Christopher (2011). Statistics for Spatio-Temporal Data. Wiley.
Donegan, Connor and Chun, Yongwan and Griffith, Daniel A. (2021). Modeling community health with areal data: Bayesian inference with survey standard errors and spatial structure. Int. J. Env. Res. and Public Health 18 (13): 6856. DOI: 10.3390/ijerph18136856 Data and code: https://github.com/ConnorDonegan/survey-HBM.
Donegan, Connor (2021). Building spatial conditional autoregressive (CAR) models in the Stan programming language. OSF Preprints. doi:10.31219/osf.io/3ey65.
Haining, Robert and Li, Guangquan (2020). Modelling Spatial and Spatial-Temporal Data: A Bayesian Approach. CRC Press.
Examples
# model mortality riskdata(georgia)C <- shape2mat(georgia, style = "B")cp <- prep_car_data(C)fit <- stan_car(deaths.male ~ offset(log(pop.at.risk.male)), car_parts = cp, data = georgia, family = poisson(), iter = 800, chains = 1 # for example speed only )rstan::stan_rhat(fit$stanfit)rstan::stan_mcse(fit$stanfit)print(fit)sp_diag(fit, georgia)## DCAR specification (inverse-distance based)library(sf)A <- shape2mat(georgia, "B")D <- sf::st_distance(sf::st_centroid(georgia))A <- D * Acp <- prep_car_data(A, "DCAR", k = 1)fit <- stan_car(deaths.male ~ offset(log(pop.at.risk.male)), data = georgia, car = cp, family = poisson(), iter = 800, chains = 1 # for example speed only )print(fit)
[Package geostan version 0.5.4 Index]