Small Area Estimation with BRFSS Data
The Behavioral Risk Factor Surveillance System (BRFSS) is an annual telephone health survey conducted by the Centers for Disease Control and Prevention (CDC) that tracks health conditions and risk behaviors in the United States and its territories since 1984. The BRFSS sampling scheme is complex with high variability in the sampling weights. In this example, we estimate the prevalence of Type II diabetes in health reporting areas (HRAs) in King County, using BRFSS data. We will compare the weighted direct estimates, with simple spatial smoothed estimates (ignore weighting), and the smoothed and weighted estimates.
Load Package and Data
First, we load the package and the necessary data. INLA is not in a standard repository, so we check if it is available and install it if it is not.
library(SUMMER)
library(survey)
library(ggplot2)
library(patchwork)
if (!isTRUE(requireNamespace("INLA", quietly = TRUE))) {
install.packages("INLA", repos=c(getOption("repos"),
INLA="https://inla.r-inla-download.org/R/stable"), dep=TRUE)
}
data(BRFSS)
data(KingCounty)BRFSS contains the full BRFSS dataset with \(16,283\) observations. The diab2 variable is the binary indicator of Type II diabetes, strata is the strata indicator and rwt_llcp is the final design weight. For the purpose of this analysis, we first remove records with missing HRA code or diabetes status from this dataset.
KingCounty contains the map of the King County HRAs. In order to fit spatial smoothing model, we first need to compute the adjacency matrix for the HRAs, mat, and make sure both the column and row names correspond to the HRA names.
## [1] 48 48The Direct Estimates
Let \(y_i\) and \(m_i\) be the number of individuals flagged as having type II diabetes and the denominators in areas \(i = 1, ..., n\). Ignoring the survey design, the naive estimates for the prevalence of Type II diabetes can be easily calculated as \(\hat p_i = y_i / m_i\), with associated standard errors \(\sqrt{\hat p_i(1 - \hat p_i)/m_i}\). The design-based weighted estimates of \(p_i\) and the associated variances can be easily calculated using the survey package.
design <- svydesign(ids = ~1, weights = ~rwt_llcp, strata = ~strata, data = BRFSS)
direct <- svyby(~diab2, ~hracode, design, svymean)
head(direct)## hracode diab2 se
## Auburn-North Auburn-North 0.104 0.021
## Auburn-South Auburn-South 0.233 0.049
## Ballard Ballard 0.070 0.022
## Beacon/Gtown/S.Park Beacon/Gtown/S.Park 0.081 0.026
## Bear Creek/Carnation/Duvall Bear Creek/Carnation/Duvall 0.052 0.012
## Bellevue-Central Bellevue-Central 0.059 0.015The Smoothed Estimates
When the number of samples in each area is large, the design-based variance is usually small and the direct estimates will work well. However, when we have small sample from each area, we would like to perform some form of smoothing over the areas. For now, let us ignore the survey weights, we can consider the following Bayesian smoothing model:
\[\begin{eqnarray*}
y_i | p_i &\sim& \mbox{Binomial}(m_i,p_i)\\
\theta_i &=& \log \left( \frac{p_i}{1-p_i}\right) = \mu+ \epsilon_i+s_i,\\
\epsilon_i &\sim & N(0,\sigma_\epsilon^2)\\
s_i | s_j, j \in \text{ne}(i) &\sim& N\left(\bar{s_i}, \dfrac{\sigma_s^2}{n_i} \right).
\end{eqnarray*}\] where \(n_i\) is the number of neighbors for area \(i\), \(\bar{s_i} = \frac{1}{n_i}\sum_{j \in \text{ne}(i)} s_j\), and hyperpriors are put on \(\mu,\sigma_\epsilon^2,\sigma_s^2\). This simple smoothing model can be fitted using the smoothSurvey function by specifying NULL for the survey parameters
The smoothed estimates of \(p_i\) and \(\theta_i\) can be found in the smooth object returned by the function, and the direct estimates are stored in the HT object (without specifying survey weights, these are the simple binomial probabilities).
## Length Class Mode
## HT 8 data.frame list
## smooth 11 data.frame list
## smooth.overall 0 -none- NULL
## fit 51 inla list
## CI 1 -none- numeric
## Amat 2304 -none- numeric
## responseType 1 -none- character
## formula 3 formula callThe Weighted and Smoothed Estimates
To account for the survey designs in the smoothing, we instead model the logit of the design-based direct estimates \(\hat p_i \sim N(\theta_i, \hat V_i)\) directly, where both \(\hat p_i\) and the asymptotic variance on the logit scale, \(\hat V_i\), are assumed known. This model can be fit with
Again, the design-based direct estimates and the smoothed estimates accounting for design can be obtained by svysmoothed$smooth and svysmoothed$HT. We can now compare the three types of estimates and their associated variance and it can be seen that smoothing reduces the variance of estimates significantly.
est <- data.frame(naive = smoothed$HT$HT.est,
weighted = svysmoothed$HT$HT.est,
smooth = smoothed$smooth$mean,
weightedsmooth = svysmoothed$smooth$mean)
var <- data.frame(naive = smoothed$HT$HT.var,
weighted = svysmoothed$HT$HT.var,
smooth = smoothed$smooth$var,
weightedsmooth = svysmoothed$smooth$var)
l1 <- range(est)
l2 <- range(var)
g1 <- ggplot(est, aes(x = naive, y = smooth)) + geom_point() +
geom_abline(slope = 1, intercept = 0, color = "red") +
ggtitle("Naive Estimates") + xlab("Direct") +
ylab("Smoothed")+ xlim(l1) + ylim(l1)
g2 <- ggplot(var, aes(x = naive, y = weightedsmooth)) + geom_point() +
geom_abline(slope = 1, intercept = 0, color = "red") +
ggtitle("Naive Variance") + xlab("Direct") +
ylab("Smoothed") + xlim(l2) + ylim(l2)
g3 <- ggplot(est, aes(x = weighted, y = weightedsmooth)) + geom_point() +
geom_abline(slope = 1, intercept = 0, color = "red") +
ggtitle("Survey Weighted Estimates") + xlab("Direct") +
ylab("Smoothed") + xlim(l1) + ylim(l1)
g4 <- ggplot(var, aes(x = weighted, y = weightedsmooth)) + geom_point() +
geom_abline(slope = 1, intercept = 0, color = "red") +
ggtitle("Survey Weighed Variance") + xlab("Direct") +
ylab("Smoothed") + xlim(l2) + ylim(l2)
(g1 + g2) / (g3 + g4)
Area-level covariates
Covariates at the area level could be included in the smoothing model by specifying the X argument. The first column of the covariate matrix should be the names of regions that match the column and row names of the adjacency matrix.
Education <- svyby(~educau, ~hracode, design, svymean)
Education <- Education[, 1:5]
svysmoothed_with_covariate <- smoothSurvey(data = BRFSS, geo = KingCounty, Amat = mat, X = Education, responseType = "binary", responseVar="diab2", strataVar="strata", weightVar="rwt_llcp", regionVar="hracode", clusterVar = "~1", CI = 0.95)
svysmoothed_with_covariate$fit$summary.fixed## mean sd 0.025quant 0.5quant 0.97quant mode
## (Intercept) -2.67 0.071 -2.81 -2.67 -2.54 -2.67
## educauless.than.HS 1.93 3.206 -4.41 1.94 7.95 1.95
## educauHS.grad -1.36 1.350 -4.02 -1.37 1.19 -1.37
## educausome.college 1.59 1.144 -0.66 1.59 3.75 1.58
## educaucollege.grad -0.74 0.507 -1.73 -0.74 0.22 -0.75
## kld
## (Intercept) 2.9e-07
## educauless.than.HS 9.1e-07
## educauHS.grad 1.3e-06
## educausome.college 5.9e-07
## educaucollege.grad 5.5e-06Prior specification
The smoothSurvey function has some default hyperprior choices built in using the Penalising Complexity (PC) priors. To change the hyperparameters of the PC prior, we can simply use the pc.u and pc.alpha for priors on precision parameters, and pc.u.phi and pc.alpha.phi for mixing parameter \(\phi\) in the BYM2 model for spatial effects. The interpretation of PC prior parameters for the precision \(\tau\) is \[
\mbox{Prob}(\sigma > u) = \alpha
\] where \(\sigma = 1 / \sqrt{\tau}\) is the standard deviation. And for the mixing parameter \(\phi\) is \[
\mbox{Prob}(\phi < u) = \alpha
\]
For example, if we change pc.u and pc.alpha from the default value \(u = 1, \alpha = 0.01\) to \(u = 0.1, \alpha = 0.01\), we would assign more prior mass on smaller variance of the random effects.
We can visualize one or more metrics on the map is by the mapPlot function, for example,
toplot <- svysmoothed$smooth
toplot$mean.new <- svysmooth.1$smooth$mean
toplot$var.new <- svysmooth.1$smooth$var
variables <- c("mean", "mean.new",
"var", "var.new")
names <- c("Mean (default prior)", "Mean (new prior)",
"Variance (default prior)", "Variance (new prior)")
g1 <- mapPlot(data = toplot, geo = KingCounty, variables=variables[1:2], labels = names[1:2], by.data = "region", by.geo = "HRA2010v2_", size = 0.1)
g2 <- mapPlot(data = toplot, geo = KingCounty, variables=variables[3:4], labels = names[3:4], by.data = "region", by.geo = "HRA2010v2_", size = 0.1)
g1 / g2
A more complex way to change the random effect model entirely is to use the formula argument. To use this option, we will need to use the internal indicator for region, region.struct or region.unstruct as the index. This index variable name can also be observed from the fitted object by summary(svysmoothed$fit) or svysmoothed$formula. For example, we can reassign a different parameterization of the BYM model. For users familiar with INLA syntax, this allows more possibilities of expanding the modeling framework.
formula <- "f(region.struct, model = 'besag', graph = Amat,
constr = TRUE, scale.model = TRUE) +
f(region.unstruct, model = 'iid')"
svysmooth.2 <- smoothSurvey(data = BRFSS, geo = KingCounty, Amat = mat,
responseType = "binary", responseVar="diab2",
strataVar="strata", weightVar="rwt_llcp", regionVar="hracode",
clusterVar = "~1", CI = 0.95,
formula = formula)
