Getting Started with saeHB.Spatial.Beta

Introduction

The saeHB.Spatial.Beta package provides several functions to estimate small area proportions using Hierarchical Bayesian (HB) methods under spatial and non-spatial models for area-level random effects. This package is specifically designed to accommodate survey design effects (DEFF) for sampling variances.

In this vignette, we will demonstrate a complete analytical workflow:

1. Preparation and Data Loading, including calculating the variance and RSE of the direct estimates to serve as a baseline.
2. Fitting the HB Beta Deff Non-Spatial Model.
3. Spatial Autocorrelation Diagnostic on the random effects using Moran’s I.
4. Fitting the HB Beta Deff Spatial Models (SAR and Leroux CAR).
5. Evaluation and Comparison of estimation precision using Relative Standard Error (RSE).

Model Specifications

Before diving into the workflow, let’s briefly review the underlying models. Let \(\hat{\theta}_i\) denote the direct estimator of the proportion for area \(i\) (\(i = 1, \dots, m\)), and \(\theta_i\) denote the true proportion parameter. Under the Hierarchical Bayesian framework, the Beta spatial model accommodating survey design effects is specified as follows:

1. Sampling Model: \[\hat{\theta}_i \mid \theta_i \sim \text{Beta}(a_i, b_i)\] To accommodate the effective sample size from complex surveys, the shape parameters are strictly defined using the area’s sample size (\(n_i\)) and survey design effect (\(d_i\)). Thus, the parameters are formulated as: \[a_i = \theta_i \left(\frac{n_i}{d_i} - 1\right) \quad \text{and} \quad b_i = (1 - \theta_i) \left(\frac{n_i}{d_i} - 1\right)\] With this specification, the expected value holds \(E(\hat{\theta}_i \mid \theta_i) = \theta_i\). The sampling variance of the direct estimator can be expressed as \(\psi_i = \left[ \frac{\theta_i (1-\theta_i)}{n_i} \right] d_i\), indicating that the sampling variance directly depends on the proportion parameter \(\theta_i\) and the design effect \(d_i\).

2. Linking Model: The true proportion \(\theta_i\) is modeled using a logit link function: \[\text{logit}(\theta_i) = \mathbf{x}_i^T \boldsymbol{\beta} + v_i\] Here, \(\mathbf{x}_i^T\) represents the vector of auxiliary variables, \(\boldsymbol{\beta}\) the regression coefficients, and \(v_i\) the area-specific random effect.

3. Spatial Random Effects: To account for spatial dependency, the random effects vector \(\mathbf{v} = (v_1, \dots, v_m)^T\) can be modeled using two different spatial structures:

1. SAR Model: The random effects follow a Simultaneous Autoregressive process, mathematically defined as:

   \[\mathbf{v} = \rho \mathbf{W} \mathbf{v} + \mathbf{u}, \quad \mathbf{u} \sim N_m(\mathbf{0}, \sigma_u^2 \mathbf{I})\]

   Assuming the matrix \((\mathbf{I} - \rho \mathbf{W})\) is non-singular, it can be rewritten explicitly as \(\mathbf{v} = (\mathbf{I} - \rho \mathbf{W})^{-1} \mathbf{u}\), which implies \(\mathbf{v} \sim N_m(\mathbf{0}, \mathbf{G})\) with the covariance dispersion matrix defined as:

   \[\mathbf{G} = \sigma_u^2 [(\mathbf{I} - \rho \mathbf{W})^T (\mathbf{I} - \rho \mathbf{W})]^{-1}\]

   Where \(\rho\) is the spatial autocorrelation parameter, \(\mathbf{W}\) is a row-standardized spatial weights matrix, and \(\mathbf{u}\) is the independent random error vector.

2. Leroux CAR Model: The random effects follow a Conditional Autoregressive structure, \(\mathbf{v} \sim N_m(\mathbf{0}, \mathbf{Q}(\rho)^{-1})\), with the precision matrix defined as:

   \[\mathbf{Q}(\rho) = \frac{1}{\sigma_v^2} [(1 - \rho)\mathbf{I} + \rho(\mathbf{D} - \mathbf{W})]\]

   Where \(\mathbf{W}\) is the binary spatial adjacency matrix, \(\mathbf{D}\) is a diagonal matrix containing the number of neighbors for each area, and \(\rho \in [0,1)\) controls the spatial dependence.

Step 1: Preparation and Data Loading

First, load the package along with the provided synthetic dataset (databeta). We will also load ggplot2 for visualization. After loading the data, we calculate the variance and the Relative Standard Error (RSE) of the direct estimates to serve as our baseline for comparison.

library(saeHB.Spatial.Beta)
library(ggplot2)

# Load data
data("databeta")

# Calculate Variance of Direct Estimator for proportion data considering DEFF
# var(y) = [y * (1 - y) / n_i] * deff
var_direct <- (databeta$y * (1 - databeta$y) / databeta$n_i) * databeta$deff

# Calculate Relative Standard Error (RSE) of Direct Estimation
databeta$rse_direct <- (sqrt(var_direct) / databeta$y) * 100

Step 2: Fitting the HB Beta Deff Non-Spatial Model

We begin by fitting a baseline non-spatial model that accommodates the survey design effect. Here, we use the default Markov Chain Monte Carlo (MCMC) iterations.

Note: By default, the function runs MCMC with a predefined number of iterations and burn-in periods. For practical applications, you may need to adjust parameters such as iter.mcmc, burn.in, thin, and chains to ensure proper mixing and convergence.

mod_ns_deff <- betadeff_nonspatial(
  formula = y ~ x1 + x2,
  deff = "deff",
  n_i = "n_i",
  data = databeta
)

Extract the RSE for the non-spatial estimates:

rse_ns_deff <- (mod_ns_deff$est$Est.Error / mod_ns_deff$est$Estimate) * 100

Step 3: Spatial Autocorrelation Diagnostic (Moran’s I)

To determine whether a spatial model is warranted, we evaluate the spatial autocorrelation of the random effects (\(v\)) obtained from the non-spatial model. We use a row-standardized spatial weight matrix for Moran’s I testing.

data("weight_mat")
W_listw <- spdep::mat2listw(weight_mat, style = "W")

# Extract the mean of the random effects (v)
v_ns_deff <- as.numeric(mod_ns_deff$randeff$Estimate)

Since our dataset has a relatively small number of areas (\(m = 36\)), we use the Monte Carlo permutation approach with 999 permutations, which computes the p-value empirically by randomly permuting the observed values. The analytical randomisation approach can be used as an alternative (mc = FALSE).

set.seed(123) 
moran_result <- moran_test(x = v_ns_deff, listw = W_listw, mc = TRUE, nsim = 999)
print(moran_result)
#> 
#>  Monte-Carlo simulation of Moran I
#> 
#> data:  v_ns_deff
#> statistic = 0.43499, observed rank = 1000, p-value = 0.001
#> alternative hypothesis: greater

A significant p-value confirms that the non-spatial model left unexplained spatial structure, heavily justifying the use of spatial models.

Step 4: Fitting the HB Beta Deff Spatial Models

We will now fit two spatial models: the Simultaneous Autoregressive (SAR) model and the Leroux Conditional Autoregressive (CAR) model. The SAR model requires a row-standardized weight matrix, while the Leroux CAR model requires a binary adjacency matrix.

# 1. Fit Spatial SAR Model
# Load the spatial weight matrix for the SAR model
data("weight_mat")
mod_sar_deff <- betadeff_sar(
  formula = y ~ x1 + x2,
  deff = "deff",
  n_i = "n_i",
  proxmat = weight_mat,
  data = databeta
)

# 2. Fit Spatial Leroux CAR Model
# Load the binary adjacency matrix for the Leroux CAR model
data("adjacency_mat")
mod_leroux_deff <- betadeff_lerouxcar(
  formula = y ~ x1 + x2,
  deff = "deff",
  n_i = "n_i",
  proxmat = adjacency_mat,
  data = databeta
)

Extract the RSE for both spatial models:

rse_sar_deff <- (mod_sar_deff$est$Est.Error / mod_sar_deff$est$Estimate) * 100
rse_leroux_deff <- (mod_leroux_deff$est$Est.Error / mod_leroux_deff$est$Estimate) * 100

Step 5: Evaluation and Comparison

We compare the performance of the Direct Estimator, the HB Beta Deff Non-Spatial Model, and the HB Beta Deff Spatial Models (SAR and Leroux CAR). A lower RSE indicates a more reliable and precise estimate.

First, let’s look at the summary statistics of the RSEs provided by the SAE models compared to the direct estimates.

We can also visualize this comparison across all areas using a line and point plot to observe the fluctuation of RSE values.

# Combine RSEs into a single data frame for plotting
df_rse <- data.frame(
  Area = seq_along(databeta$y),
  Direct = databeta$rse_direct,
  Non_Spatial = rse_ns_deff,
  Spatial_SAR = rse_sar_deff,
  Spatial_Leroux = rse_leroux_deff
)

# Order by Direct RSE for better visualization
df_rse <- df_rse[order(df_rse$Direct), ]
df_rse$Area_Index <- seq_len(nrow(df_rse))

# Plotting the RSE Comparison
ggplot(df_rse, aes(x = Area_Index)) +
  # Direct Estimation
  geom_line(aes(y = Direct, color = "Direct"), linewidth = 0.8, alpha = 0.6) +
  geom_point(aes(y = Direct, color = "Direct"), size = 2, alpha = 0.6) +
  
  # Non-Spatial
  geom_line(aes(y = Non_Spatial, color = "HB Beta Deff Non-Spatial"), linewidth = 0.8, alpha = 0.8) +
  geom_point(aes(y = Non_Spatial, color = "HB Beta Deff Non-Spatial"), size = 2, alpha = 0.8) +
  
  # Spatial SAR
  geom_line(aes(y = Spatial_SAR, color = "HB Beta Deff Spatial SAR"), linewidth = 1) +
  geom_point(aes(y = Spatial_SAR, color = "HB Beta Deff Spatial SAR"), size = 2) +
  
  # Spatial Leroux CAR
  geom_line(aes(y = Spatial_Leroux, color = "HB Beta Deff Spatial Leroux CAR"), linewidth = 1) +
  geom_point(aes(y = Spatial_Leroux, color = "HB Beta Deff Spatial Leroux CAR"), size = 2) +
  
  scale_color_manual(
    name = "Estimator",
    values = c("Direct" = "#E69F00", 
               "HB Beta Deff Non-Spatial" = "#56B4E9", 
               "HB Beta Deff Spatial SAR" = "#009E73",
               "HB Beta Deff Spatial Leroux CAR" = "#D55E00")
  ) +
  labs(
    title = "Comparison of Relative Standard Error (RSE)",
    subtitle = "Lower RSE indicates higher precision",
    x = "Area (Ordered by Direct RSE)",
    y = "RSE (%)"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")