--- title: "Spatial and Spatio-Temporal Models" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Spatial and Spatio-Temporal Models} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5 ) ``` ## Overview When small areas are geographic units (such as counties, regencies, or districts), spatial correlation frequently occurs: neighboring areas tend to have more similar outcomes than distant ones. Furthermore, when surveys are repeated over multiple years, temporal correlation between time periods is also present. **fastsae** provides specialized, high-performance implementations for both: 1. **Spatial Fay-Herriot (`eblup_sfh`)**: Simultaneous Autoregressive SAR(1) random effects. 2. **Spatio-Temporal Fay-Herriot (`eblup_stfh`)**: Combined spatial SAR(1) and temporal AR(1) random effects. --- ## 1. Spatial Fay-Herriot Model (`eblup_sfh`) ### Model Formulation The Spatial Fay-Herriot model (Pratesi & Salvati, 2008) models area random effects using a Simultaneous Autoregressive (SAR(1)) process: $$y = X\beta + u + e$$ $$u = \rho_1 W u + \epsilon_1 \implies u = (I - \rho_1 W)^{-1} \epsilon_1$$ where: - $W$ is a known row-standardized spatial proximity matrix ($D \times D$). - $\rho_1 \in (-1, 1)$ is the spatial autoregressive parameter. - $\epsilon_1 \sim N(0, \sigma_{u1}^2 I_D)$ is independent innovation noise. - $e \sim N(0, V_e)$ is sampling error with diagonal matrix $V_e = \text{diag}(D_1, \dots, D_D)$. ### Fitting the Model ```{r sfh_fit} library(fastsae) data("mys") data("mys_proxmat") # Fit Spatial Fay-Herriot model with REML fit_sfh <- eblup_sfh( formula = y ~ x1 + x2 + x3, vardir = ~vardir, data = mys, W = mys_proxmat, method = "REML", print_result = FALSE ) summary(fit_sfh) ``` ### Unsampled Domains and Spatial Kriging A major advantage of `eblup_sfh` in **fastsae** is its automatic handling of unsampled domains (domains with `y = NA`). Instead of throwing an error or requiring manual subsetting, `eblup_sfh` automatically performs **full-spatial kriging**: ```{r unsampled} # Count sampled vs unsampled domains table(is.na(mys$y)) # Inspect estimates for unsampled domains head(fit_sfh$df_eblup[is.na(mys$y), ]) ``` For unsampled areas, prediction borrows strength from both the regression synthetic component $x_d^\top \hat{\beta}$ and spatial proximity to neighboring sampled areas through $(I - \hat{\rho}_1 W)^{-1}$. ### Multi-Threaded Bootstrap MSE In addition to analytical Prasad-Rao style MSE, `eblup_sfh` supports multi-threaded OpenMP bootstrap MSE estimation: - **Parametric Bootstrap (`mse_method = "pbmse"`)** - **Non-Parametric Bootstrap (`mse_method = "npbmse"`)** ```{r sfh_pbmse} fit_sfh_pb <- eblup_sfh( formula = y ~ x1 + x2 + x3, vardir = ~vardir, data = mys, W = mys_proxmat, mse_method = "pbmse", B = 50, n_threads = 2, seed = 123, print_result = FALSE ) head(fit_sfh_pb$df_eblup[, c("domain", "y", "eblup", "mse", "mse_pb", "mse_pbbc")]) ``` --- ## 2. Spatio-Temporal Fay-Herriot Model (`eblup_stfh`) ### Model Formulation The Spatio-Temporal model (Marhuenda, Molina, & Morales, 2013) combines spatial correlation and temporal autoregression for $D$ areas over $T$ time periods. The panel does not need to be complete: cells without a direct estimate are predicted (see below). $$y_{dt} = x_{dt}^\top \beta + u_{1d} + u_{2dt} + e_{dt}$$ where: - $u_1 = (u_{11}, \dots, u_{1D})^\top$ captures static spatial effects: $u_1 = \rho_1 W u_1 + \epsilon_1$. - $u_{2d} = (u_{2d1}, \dots, u_{2dT})^\top$ captures dynamic temporal effects for area $d$: $u_{2dt} = \rho_2 u_{2d,t-1} + \epsilon_{2dt}$. - $e_{dt} \sim \text{ind. } N(0, D_{dt})$ are sampling errors. The model is defined on the full $D \times T$ grid: $u_1$ on all domains of $W$, $u_2$ on all periods. ### Fitting the Spatio-Temporal Model ```{r stfh_fit} data("mys_panel") # 42 domains x 11 years; y is NA for the two domains that are never sampled table(sampled = !is.na(mys_panel$y)) # Fit Spatio-Temporal model with the analytical MSE, on the full panel and full W fit_stfh <- eblup_stfh( formula = y ~ x1 + x2 + x3, data = mys_panel, vardir = ~vardir, domain = ~area, time = ~year, W = mys_proxmat, model = "ST", compute_mse = TRUE, print_result = FALSE ) summary(fit_stfh) ``` The estimated variance and correlation components include: - `sigma21`: Spatial random effect variance ($\sigma_{u1}^2$). - `rho1`: Spatial autocorrelation ($\rho_1$). - `sigma22`: Temporal random effect variance ($\sigma_{u2}^2$). - `rho2`: Temporal AR(1) autoregression coefficient ($\rho_2$). The analytical MSE is $g_1 + g_2 + 2 g_3 - g_4$ (Prasad-Rao type, with the exact derivative of the EBLUP in $g_3$ and the curvature term $g_4$ for the nonlinear $\rho_1$, $\rho_2$); the components are returned as `fit_stfh$g1`, ..., `fit_stfh$g4`. A parametric bootstrap MSE is available with `mse_type = "bootstrap"` (and `B`, `seed`). ### Unsampled Domains and Missing Periods Rows whose response is `NA` are kept and predicted, as long as the auxiliary variables are known there. The variance components are estimated from the sampled cells only (the marginal likelihood with the missing cells integrated out), and each missing cell gets the best linear unbiased predictor given the sampled data (Saei & Chambers, 2005; Benedetti et al., 2024): - a **missing period** of a sampled domain borrows from the other periods of that domain (AR(1) correlation) and from the neighbouring domains (SAR); - a **domain that is never sampled** is predicted from its neighbours only (spatial kriging), so `W` must keep it. ```{r stfh_unsampled} est <- fit_stfh$df_eblup # Cells predicted without a direct estimate head(est[!est$sampled, c("domain", "time", "eblup", "mse", "rse")]) # Their MSE is larger: no direct estimate, only borrowed strength tapply(est$mse, ifelse(est$sampled, "sampled", "unsampled"), median) ``` For these cells the analytical MSE is $g_1 + g_2 + g_3$ ($g_3$ from the exact derivative of the predictor in the variance components); the bootstrap generates the random effects on the full grid and covers every cell too. --- ## 3. Comparing Models with `autoplot` You can directly compare EBLUP estimates and uncertainty across different models using `autoplot()`: ```{r compare_models} # Compare Fay-Herriot vs Spatial Fay-Herriot fit_fh <- eblup_fh(y ~ x1 + x2 + x3, vardir = ~vardir, data = na.omit(mys), print_result = FALSE) W_clean <- mys_proxmat[!is.na(mys$y), !is.na(mys$y)] fit_sfh_clean <- eblup_sfh(y ~ x1 + x2 + x3, vardir = ~vardir, data = na.omit(mys), W = W_clean, print_result = FALSE) autoplot(list("Standard FH" = fit_fh, "Spatial FH" = fit_sfh_clean), type = "comparison") ``` --- ## 4. Bayesian Spatio-Temporal Modeling (`hb_area`) In addition to frequentist C++ EBLUP implementations (`eblup_sfh` and `eblup_stfh`), **fastsae** provides `hb_area()` for Bayesian spatio-temporal modeling powered by **INLA**. This supports: - Bounded proportions via Beta regression (matching `tipsae`) - Discrete counts via Poisson and Negative Binomial regression - Binary aggregates via Binomial regression - Flexible space-time structures including domain-specific random walks, separable dynamic fields, and Knorr-Held Types I–IV. For full details and a case study on Italian poverty data, see the vignette [Bayesian Area-Level Small Area Estimation with INLA](hb-area-inla.html). --- ## References - Benedetti, R., Piersimoni, F., Pratesi, M., Salvati, N., & Suesse, T. (2024). Handling out-of-sample areas to estimate the unemployment rate at local labour market areas in Italy. *International Statistical Review*. https://doi.org/10.1111/insr.12596 - Marhuenda, Y., Molina, I., & Morales, D. (2013). Small area estimation with spatio-temporal Fay-Herriot models. *Computational Statistics & Data Analysis*, 58, 308–325. - Pratesi, M., & Salvati, N. (2008). Small area estimation for spatially correlated data: A Fay-Herriot with the spatial linear spline model. *Journal of Applied Statistics*, 35(7), 781–794. - Saei, A., & Chambers, R. (2005). Out of sample estimation for small areas using area level data. Southampton Statistical Sciences Research Institute, Methodology Working Paper M05/11.