Spatial measurement error models

0.1 Modeling errors of observation

The measurement error (ME) models implemented in geostan are hierarchical Bayesian models (HBMs) that incorporate two sources of information:

1. A sampling distribution for the survey estimates.
2. Spatial autocorrelation and generic background knowledge on social variables.

The model treats the true value \(\boldsymbol x\) as an unknown parameter or a latent variable. The goal is to obtain a probability distribution for \(\boldsymbol x\). That information can then be incorporated into any of geostan’s regression or disease mapping models.

The sampling distribution for the model states that the true value \(x_i\) is probably within two standard errors of the estimate, or \[z_i = x_i + \epsilon_i, \epsilon_i \sim Normal(0, s_i^2),\] where \(z_i\) is the estimate, \(s_i\) is its standard error, and \(x_i\) is the true value.

Our background knowledge on \(x\) is the second component of the model. This is simply the knowledge that social variables tend to be spatially patterned: relatively extreme values are not implausible but they tend to be clustered together. This information is encoded into a prior probability distribution for the unknown values \(\boldsymbol x\). Our prior distribution is the spatial conditional autoregressive (CAR) model (see Donegan 2022 Table 1, WCAR specification).

If \(\boldsymbol x\) is a rate variable—e.g., the poverty rate—it can be important to use a logit-transformation before applying the CAR prior (Logan et al. 2020).

The CAR model has three scalar parameters (mean \(\mu\), spatial dependence \(\rho\), and scale \(\tau\)) that require prior probability distributions; geostan uses the following by default:

\[\begin{equation} \begin{split} \mu &\sim Gauss(0, 100) \\ \tau &\sim Student_t(10, 0, 40) \\ \rho &\sim Uniform(\text{lower_bound}, \text{upper_bound}) \end{split} \end{equation}\]

The default prior for \(\rho\) is uniform across its entire support (determined by the extreme eigenvalues of the spatial weights matrix).

0.2 Getting started

From the R console, load the geostan package.

library(geostan)
data(georgia)

The line data(georgia) loads the georgia data set from the geostan package into your working environment. You can learn more about the data by entering ?georgia to the R console.

0.3 Preparing the data

Users can add a spatial ME model to any geostan model. A list of data for the ME models can be prepared using the prep_me_data function. The function requires at least two inputs from the user:

1. se A data frame with survey standard errors.
2. car_parts A list of data created by the prep_car_data function.

It also has optional inputs:

1. prior A list of prior distributions for the CAR model parameters.
2. logit A vector of logical values (TRUE or FALSE) specifying if the variate should be logit transformed. Only use this for rates (proportions between 0 and 1).

The default prior distributions are designed to be fairly vague relative to ACS variables (including percentages from 0-100), assuming that magnitudes such as income have been log-transformed already.

The logit-transformation is often required for skewed rates, such as the poverty rate. The user does not do this transformation on their own (i.e., before passing the data to the model), simply use the logit argument in prep_me_data.¹

As an example, we will be using estimates of the percent covered by health insurance in Georgia counties. We are going to convert this to a rate so that we can use the logit transform:

data(georgia)
georgia$insurance <- georgia$insurance / 100
georgia$insurance.se <- georgia$insurance.se / 100

Now we need to gather the standard errors georgia$insurance.se into a data frame. The column name for the standard errors must match the name of the variable it refers to (“insurance”):

SE <- data.frame(insurance = georgia$insurance.se)

If we intended to model additional survey-based covariates, we would simply add their standard errors to this data frame as another column.

We create a binary spatial connectivity matrix, and then pass it to prep_car_data to prepare it for the CAR model:

C <- shape2mat(georgia, "B", quiet = TRUE)
cars <- prep_car_data(C)

## Range of permissible rho values: -1.661, 1

Now we use prep_me_data to combine these items. We are going to use the logit-transform, since we are working with rates:

ME_list <- prep_me_data(se = SE,
                        car_parts = cars,
                        logit = TRUE
                        )

Now we are ready to begin modeling the data.

0.4 Spatial ME model

The following code chunk sets up a spatial model for male county mortality rates (ages 55-64) using the insurance variable as a covariate (without the spatial ME model):

fit <- stan_car(deaths.male ~ offset(log(pop.at.risk.male)) + insurance,
                centerx = TRUE,
                data = georgia, 
                family = poisson(),
                car_parts = cars)

To add the spatial ME model to the above specification, we simply provide our ME_list to the ME argument:

fit <- stan_car(deaths.male ~ offset(log(pop.at.risk.male)) + insurance,
                centerx = TRUE,
                ME = ME_list,
                family = poisson(),
                data = georgia, 
                car_parts = cars,
                iter = 650, # for demo speed 
                quiet = TRUE, 
                )

print(fit)

## Spatial Model Results 
## Formula: deaths.male ~ offset(log(pop.at.risk.male)) + insurance
## Spatial method (outcome):  CAR 
## Likelihood function:  poisson 
## Link function:  log 
## Residual Moran Coefficient:  -0.0006792308 
## WAIC:  1321.42 
## Observations:  159 
## Data models (ME): insurance
##  Data model (ME prior): CAR (auto Gaussian)
## Inference for Stan model: foundation.
## 4 chains, each with iter=650; warmup=325; thin=1; 
## post-warmup draws per chain=325, total post-warmup draws=1300.
## 
##             mean se_mean    sd   2.5%    20%    50%    80%  97.5% n_eff  Rhat
## intercept -4.165   0.002 0.051 -4.274 -4.201 -4.164 -4.129 -4.069  1062 1.003
## insurance -2.718   0.021 0.738 -4.308 -3.336 -2.681 -2.091 -1.313  1245 1.001
## car_rho    0.860   0.003 0.087  0.663  0.790  0.876  0.937  0.982  1110 1.007
## car_scale  0.449   0.001 0.032  0.391  0.421  0.449  0.474  0.518  1051 1.003
## 
## Samples were drawn using NUTS(diag_e) at Wed Sep 18 07:07:05 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 CAR models in geostan can be highly efficient in terms of MCMC sampling, but the required number of iterations will depend on many characteristics of the model and data. Often the default iter = 2000 is more than sufficient (with cores = 4). To speed up sampling with multi-core processing, use cores = 4 (to sample 4 chains in parallel).

In the following section, methods for examining MCMC samples of the modeled covariate values \(\boldsymbol x\) will be illustrated. Note that the process of storing MCMC samples for \(\boldsymbol x\) can become computationally burdensome if you have multiple covariates and moderately large N. If you do not need the samples of \(\boldsymbol x\) you can use the drop argument, as in:

fit2 <- stan_car(deaths.male ~ offset(log(pop.at.risk.male)) + insurance,
                centerx = TRUE,
                ME = ME_list,
        drop = 'x_true', 
                family = poisson(),
                data = georgia, 
                car_parts = cars
                )

Using slim = TRUE may be faster (particularly for larger data sets), but will prevent the collection of MCMC samples for quantities of interest such as fitted values, pointwise log-likelihood values (required for WAIC), as well as modeled covariates \(\boldsymbol x\).

0.5 Visual diagnostics

The me_diag provides some useful diagnostics for understanding how the raw estimates compare with the modeled values.

First, it provides a point-interval scatter plot that compares the raw survey estimates (on the horizontal axis) to the modeled values (on the vertical axis). The (posterior) probability distributions for the modeled values are represented by their 95% credible intervals. The intervals provide a sense of the quality of the ACS data—if the credible intervals on the modeled value are large, this tells us that the data are not particularly reliable.

The me_diag function also provides more direct visual depictions of the difference between the raw survey estimates \(z_i\) an the modeled values \(x_i\): \[\delta_i = z_i - x_i.\] We want to look for any strong spatial pattern in these \(\delta_i\) values, because that would be an indication of a bias. However, the magnitude of the \(\delta_i\) value is important to consider—there may be a pattern, but if the amount of shrinkage is very small, that pattern may not matter. (The model returns \(M\) samples from the posterior distribution of each \(x_i\); or, indexing by MCMC sample \(x_i^m\) (\(m\) is an index, not an exponent). The reported \(\delta_i\) values is the MCMC mean \(\delta_i = \frac{1}{M} \sum_m z_i - x_i^m\).)

Two figures are provided to evaluate patterns in \(\delta_i\): first, a Moran scatter plot and, second, a map.

me_diag(fit, 'insurance', georgia)

In this case, the results do not look too concerning insofar as there are no conspicuous patterns. However, a number of the credible intervals on the modeled values are large, which indicates low data reliability. The fact that some of the \(\delta_i\) values are substantial also points to low data quality for those estimates.

It can be important to understand that the ME models are not independent from the rest of the model. When the modeled covariate \(\boldsymbol x\) enters a regression relationship, that regression relationship can impact the probability distribution for \(\boldsymbol x\) itself. (This is a feature of all Bayesian ME models). To examine the ME model in isolation, you can set the prior_only argument to TRUE (see ?stan_glm or any of the spatial models).

mu.x <- as.matrix(fit, pars = "mu_x_true")
dim(mu.x)

## [1] 1300    1

head(mu.x)

##           parameters
## iterations mu_x_true[1]
##       [1,]     1.852153
##       [2,]     1.787614
##       [3,]     1.781939
##       [4,]     1.787128
##       [5,]     1.731832
##       [6,]     1.885786

mean(mu.x)

## [1] 1.790968

print(fit$stanfit, pars = c("mu_x_true", "car_rho_x_true", "sigma_x_true"))

## Inference for Stan model: foundation.
## 4 chains, each with iter=650; warmup=325; thin=1; 
## post-warmup draws per chain=325, total post-warmup draws=1300.
## 
##                   mean se_mean   sd 2.5%  25%  50%  75% 97.5% n_eff Rhat
## mu_x_true[1]      1.79       0 0.09 1.62 1.75 1.79 1.83  1.95   429    1
## car_rho_x_true[1] 0.91       0 0.07 0.75 0.88 0.92 0.96  0.99  1101    1
## sigma_x_true[1]   0.48       0 0.04 0.42 0.46 0.48 0.51  0.56  1176    1
## 
## Samples were drawn using NUTS(diag_e) at Wed Sep 18 07:07:05 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).

x <- as.matrix(fit, pars = "x_true")
dim(x)

## [1] 1300  159

x.mu <- apply(x, 2, mean)
head(x.mu)

## x_insurance[1] x_insurance[2] x_insurance[3] x_insurance[4] x_insurance[5] 
##      0.8365083      0.8002242      0.8518988      0.8522241      0.9178264 
## x_insurance[6] 
##      0.8702460

res <- resid(fit)$mean
plot(x.mu, res,
     xlab = 'Insurance rate',
     ylab = 'Residual',
     type = 'n',
     bty = 'n')
abline(h = 0)
points(x.mu, res,
       col = 4,
       pch = 20)

0.6 Non-spatial ME models

If the ME list doesn’t have a slot with car_parts, geostan will automatically use a non-spatial Student’s t model instead of the CAR model: \[ p( x | \mathcal{M}) = Student(\boldsymbol x | \nu, \mu, \sigma), \] with degrees of freedom \(\nu\), mean \(\mu\), and scale \(\sigma\).

ME_nsp <- prep_me_data(
  se = data.frame(insurance = georgia$insurance.se),
  logit = TRUE
)
fit_nsp <- stan_glm(log(rate.male) ~ insurance, data = georgia, ME = ME_nsp)

0.7 Multiple covariates

To model multiple covariates, simply add them to the data frame of standard errors:

georgia$college <- georgia$college / 100
georgia$college.se <- georgia$college.se / 100

se = data.frame(insurance = georgia$insurance.se,
                college = georgia$college.se)

ME <- prep_me_data(se = se, logit = c(TRUE, TRUE))

fit <- stan_glm(deaths.male ~ offset(log(pop.at.risk.male)) + insurance + college,
                ME = ME,
        re = ~ GEOID,
                family = poisson(),
                data = georgia, 
                iter = 700
                )       

Results can be examined one at a time usign me_diag(fit, 'insurance', georgia) and me_diag(fit, 'college', georgia).

0.8 Log transforms

Income and other magnitudes are often log transformed. In that case, the survey standard errors also need to be transformed. The se_log function provide approximate standard errors for this purpose.

When using se_log, pass in the variable itself (untransformed) and its standard errors. Here is a full workflow for a spatial mortality model with covariate measurement error:

data(georgia)

# income in $1,000s
georgia$income <- georgia$income / 1e3
georgia$income.se <- georgia$income.se /1e3

# create log income
georgia$log_income <- log( georgia$income )

# create SEs for log income
log_income_se <- se_log( georgia$income, georgia$income.se )

# prepare spatial CAR data
C <- shape2mat(georgia, "B")
cars <- prep_car_data(C)

# prepare ME data 
se <- data.frame( log_income = log_income_se )
ME <- prep_me_data( se = se, car_parts = cars )

# fit model 
fit <- stan_car(deaths.male ~ offset(log(pop.at.risk.male)) + log_income,
                ME = ME,
        car_parts = cars,
        family = poisson(),
        data = georgia,
        cores = 4
                )

# check ME model
me_diag(fit, 'log_income', georgia)

# check disease model
sp_diag(fit, georgia)

# coefficient estimates
plot(fit)

# plot the income-mortality gradient
values <- seq( min(georgia$log_income), max(georgia$log_income), length.out = 100 )
new_data <- data.frame(pop.at.risk.male = 1,
                       log_income = values)

preds <- predict(fit, new_data, type = 'response')
preds$Income <- exp( preds$log_income )
preds$Mortality <- preds$mean * 100e3

plot(preds$Income, preds$Mortality,
    type = 'l',
    xlab = "Income",
    ylab = "Deaths per 100,000")