--- title: "Sub-Area Level Estimation with Two-Fold Models" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Sub-Area Level Estimation with Two-Fold Models} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5 ) # INLA is not on CRAN: the hb_twofold() chunks are shown but not run without it has_inla <- requireNamespace("INLA", quietly = TRUE) ``` ## Overview Surveys are often designed for large areas (provinces, regions) but the estimates are needed for the sub-areas inside them (districts, municipalities). The sub-areas of one area share the area's characteristics, so their estimates should borrow strength not only from the regression on auxiliary variables, but also from the other sub-areas of the same area. The two-fold sub-area level model of Torabi and Rao (2014) does this with two nested random effects. **fastsae** fits it in two ways: | Function | Approach | Uncertainty | Area aggregates | | :--- | :--- | :--- | :--- | | `eblup_twofold()` | EBLUP, REML or ML | Analytical (Prasad-Rao type) or parametric bootstrap MSE | `df_area`, weighted by `weight` | | `hb_twofold()` | Hierarchical Bayes with INLA | Posterior SD and credible intervals | `df_area`, weighted by `weight` | `hb_twofold()` also fits binomial and Poisson responses and a spatial (BYM2, Besag) area effect. --- ## Model For sub-area $j$ ($j = 1, \dots, N_i$) of area $i$ ($i = 1, \dots, m$), the direct estimator $y_{ij}$ with known sampling variance $\psi_{ij}$ follows $$y_{ij} = \theta_{ij} + e_{ij}, \qquad \theta_{ij} = x_{ij}^\top \beta + v_i + u_{ij},$$ with $v_i \sim N(0, \sigma_v^2)$ the area effect, $u_{ij} \sim N(0, \sigma_u^2)$ the sub-area effect and $e_{ij} \sim N(0, \psi_{ij})$ the sampling error, all independent. Compared with the Fay-Herriot model, which treats every sub-area as unrelated, the area effect $v_i$ is estimated from all the sampled sub-areas of area $i$. This matters most for sub-areas that are not sampled: they get $x_{ij}^\top \hat\beta + \hat v_i$ instead of the synthetic $x_{ij}^\top \hat\beta$ alone. Only a sub-area of an area without any sample is predicted synthetically. The output follows the notation above: `random_effect_var` holds `sigma2_v` (area) and `sigma2_u` (sub-area), and `df_eblup` (or `df_hb`) holds `random_effect_area` and `random_effect_subarea`. --- ## Example with simulated data We simulate 30 areas with 3 to 6 sub-areas each, so that the true sub-area means $\theta_{ij}$ are known. Fifteen sub-areas have no sample (their direct estimate is `NA`), and all the sub-areas of area `A30` are unsampled. ```{r simulate} library(fastsae) set.seed(2024) m <- 30 dat <- do.call(rbind, lapply(seq_len(m), function(i) { ni <- sample(3:6, 1) data.frame( area = sprintf("A%02d", i), subarea = sprintf("A%02d-%d", i, seq_len(ni)), x1 = rnorm(ni, mean = 2), pop = sample(200:2000, ni) ) })) v <- rnorm(m, sd = 1) names(v) <- unique(dat$area) dat$theta <- 10 + 2 * dat$x1 + v[dat$area] + rnorm(nrow(dat), sd = sqrt(0.5)) dat$vardir <- runif(nrow(dat), 0.5, 2) dat$y <- dat$theta + rnorm(nrow(dat), sd = sqrt(dat$vardir)) # Unsampled sub-areas: 15 at random, plus every sub-area of A30 unsampled <- union(sample(which(dat$area != "A30"), 15), which(dat$area == "A30")) dat$y[unsampled] <- NA c(subareas = nrow(dat), sampled = sum(!is.na(dat$y))) head(dat) ``` Each row is one sub-area. Unsampled sub-areas stay in the data with `y = NA` and their covariates, so that they are predicted. ### Fitting the EBLUP ```{r fit_eblup} fit_tf <- eblup_twofold( y ~ x1, vardir = "vardir", domain = "area", subarea = "subarea", data = dat, weight = "pop", print_result = FALSE ) fit_tf$random_effect_var fit_tf$estcoef head(fit_tf$df_eblup) ``` The default `mse_type = "analytical"` gives the Prasad-Rao type MSE of the REML fit (for ML it adds the bias correction). Sub-area `A01-3` is not sampled: it still gets the area effect of `A01`, while its sub-area effect is 0. `weight` (here the population of each sub-area) aggregates the sub-area EBLUPs to area means, $\tilde\theta_i = \sum_j w_{ij}\tilde\theta_{ij} / \sum_j w_{ij}$, in `df_area`. Their MSE accounts for the errors the sub-areas of an area share (the estimated area effect and $\hat\beta$), so it is larger than $\sum_j w_{ij}^2\,\mathrm{mse}_{ij}$ would suggest. ```{r area_eblup} head(fit_tf$df_area) ``` ### How much does the area effect help? Because the data are simulated, we can compare each estimator with the true $\theta_{ij}$. The Fay-Herriot model, fitted with `eblup_fh()` on the same rows, ignores the areas. ```{r compare_truth} fit_fh <- eblup_fh(y ~ x1, vardir = "vardir", data = dat, print_result = FALSE) sampled <- !is.na(dat$y) in_sampled_area <- !sampled & dat$area != "A30" emp_mse <- function(est, rows) mean((est[rows] - dat$theta[rows])^2) data.frame( rows = c("sampled sub-areas", "unsampled sub-areas of sampled areas"), direct = c(emp_mse(dat$y, sampled), NA), fay_herriot = c( emp_mse(fit_fh$df_eblup$eblup, sampled), emp_mse(fit_fh$df_eblup$eblup, in_sampled_area) ), two_fold = c( emp_mse(fit_tf$df_eblup$eblup, sampled), emp_mse(fit_tf$df_eblup$eblup, in_sampled_area) ) ) ``` Both models improve on the direct estimates of the sampled sub-areas, and in this data set they are about equally accurate there: a sampled sub-area has its own direct estimate, so the area effect adds little. The difference is in the unsampled sub-areas of sampled areas. The Fay-Herriot model can only give them the synthetic estimate, while the two-fold model adds the estimated area effect and roughly halves the error. This is one simulated data set, not a simulation study: the gain depends on how large $\sigma_v^2$ is relative to $\sigma_u^2$. ### Bootstrap MSE The parametric bootstrap refits the model on `B` data sets generated from the fitted model and predicts with the same estimator as the point estimate. `seed` makes it reproducible without changing the random number stream of the session. ```{r bootstrap} fit_boot <- eblup_twofold( y ~ x1, vardir = "vardir", domain = "area", subarea = "subarea", data = dat, mse_type = "bootstrap", B = 100, seed = 1, print_result = FALSE ) summary(fit_boot$df_eblup$mse / fit_tf$df_eblup$mse) ``` The two MSE estimators agree on average; with `B = 100` single bootstrap values are still noisy. Use `mse_type = "none"` when only the point estimates are needed. ### Diagnostics and plots The usual tools of the package work on the result: ```{r diagnostics} autoplot(fit_tf, type = "comparison") diagnose(fit_tf) ``` --- ## Hierarchical Bayes with `hb_twofold()` `hb_twofold()` fits the same model with INLA. It needs the **INLA** package, which is not on CRAN; see . The chunks below are only run when INLA is installed. The argument `weight` (here the population of each sub-area) aggregates the sub-area posteriors to area estimates, $\hat\theta_i = \sum_j w_{ij} \theta_{ij} / \sum_j w_{ij}$, with their posterior uncertainty in `df_area`. For the Gaussian model these are linear combinations of the sub-area predictors, whose posterior INLA computes exactly; the binomial and Poisson aggregates come from posterior draws (`inla_engine(draws = )`). ```{r fit_hb, eval = has_inla} fit_hb <- hb_twofold( y ~ x1, vardir = "vardir", domain = "area", subarea = "subarea", weight = "pop", data = dat, print_result = FALSE ) fit_hb$random_effect_var head(fit_hb$df_hb) head(fit_hb$df_area) ``` With the default priors the posterior means are very close to the EBLUP (mean difference, limits of agreement and Lin's concordance correlation below): ```{r compare_hb, eval = has_inla} compare_sae(fit_tf, fit_hb, names = c("EBLUP two-fold", "HB two-fold")) ``` `compare_sae()` reports smaller posterior variances than EBLUP MSEs. The two are not the same quantity, and here the posterior mean of $\sigma_u^2$ is also smaller than the REML estimate (both are below the true value 0.5, which is hard to separate from sampling variances of 0.5 to 2). A smaller reported uncertainty is therefore not evidence of a more accurate estimator. With simulated data, the error against the truth settles it; here the two are practically equal: ```{r truth_hb, eval = has_inla} c( eblup_two_fold = emp_mse(fit_tf$df_eblup$eblup, sampled | in_sampled_area), hb_two_fold = emp_mse(fit_hb$df_hb$hb, sampled | in_sampled_area) ) ``` For counts or proportions, use `family = "poisson"` with `exposure` or `family = "binomial"` with `trials`; with `spatial = "bym2"` and a neighbourhood graph `W` between areas (for example from `create_weights()`), the area effects are spatially correlated. ### MCMC with `engine = jax_engine()` The same model can be sampled by MCMC (NumPyro's No-U-Turn sampler on JAX) instead of INLA. It needs Python with `jax` and `numpyro` (see `?setup_numpyro_env`), so this chunk is not run here. The result has the same structure, plus the posterior draws in `draws` and the largest split R-hat in `rhat_max`; the area aggregates in `df_area` come from the MCMC draws. ```{r fit_hb_jax, eval = FALSE} fit_jax <- hb_twofold( y ~ x1, vardir = "vardir", domain = "area", subarea = "subarea", weight = "pop", data = dat, engine = jax_engine(warmup = 1000, samples = 1000, chains = 4), print_result = FALSE ) fit_jax$rhat_max head(fit_jax$df_area) ``` Both engines use the same priors by default (`jax_engine(priors = "inla")`), so their estimates differ only by INLA's approximation and the Monte Carlo error. ### Proportions from complex surveys A proportion $\hat p_{ij}$ estimated from a clustered, weighted sample is less precise than one from a simple random sample of the same size $n_{ij}$. The binomial/logit-normal model with *effective sample sizes* (Chen, Wakefield and Lumley 2014; Franco and Bell 2015) takes this into account: with the design effect $\mathrm{deff}_{ij}$, the effective sample size is $n^*_{ij} = n_{ij} / \mathrm{deff}_{ij}$ and the effective number of cases $y^*_{ij} = n^*_{ij} \hat p_{ij}$, and $$y^*_{ij} \mid p_{ij} \sim \text{Binomial}(n^*_{ij}, p_{ij}), \qquad \text{logit}(p_{ij}) = x_{ij}^\top \beta + v_i + u_{ij},$$ so that the binomial variance $p_{ij}(1 - p_{ij}) / n^*_{ij}$ matches the design variance of $\hat p_{ij}$. Neither $n^*_{ij}$ nor $y^*_{ij}$ is a whole number. `hb_twofold()` (and `hb_area()`) take the proportions as the response and the effective sample sizes as `trials`, and do not round them: INLA then uses its `"xbinomial"` likelihood, which equals the binomial likelihood on whole numbers. A direct estimate of 0 is a sampling zero and enters the model as such. ```{r neff_data} set.seed(7) prov <- rep(1:12, times = sample(4:8, 12, replace = TRUE)) svy <- data.frame( prov = prov, kab = seq_along(prov), x1 = rnorm(length(prov)), pop = round(runif(length(prov), 2e4, 2e5)) ) p_true <- plogis(-2.5 + 0.5 * svy$x1 + rnorm(12, sd = 0.4)[prov] + rnorm(nrow(svy), sd = 0.3)) n <- sample(30:200, nrow(svy), replace = TRUE) svy$deff <- runif(nrow(svy), 1.3, 2.5) # from the survey, e.g. design-based variance estimates svy$n_eff <- n / svy$deff # clustered sample: the direct estimate varies like a binomial proportion on n_eff svy$p_hat <- rbinom(nrow(svy), round(svy$n_eff), p_true) / round(svy$n_eff) sum(svy$p_hat == 0) # districts with a direct estimate of 0 ``` ```{r neff_fit, eval = has_inla} fit_neff <- hb_twofold( p_hat ~ x1, domain = "prov", subarea = "kab", weight = "pop", data = svy, family = "binomial", trials = "n_eff", print_result = FALSE ) head(fit_neff$df_hb[, c("domain", "subarea", "y", "hb", "sd", "ci_lower", "ci_upper")]) head(fit_neff$df_area) # the districts with a direct estimate of 0 fit_neff$df_hb[svy$p_hat == 0, c("subarea", "hb", "ci_lower", "ci_upper")] ``` The model is only as good as the design effects: they are not estimated by the model, so effective sample sizes that are too large give intervals that are too narrow. Design effects of small areas are themselves noisy, and are often smoothed (for example, pooled within a larger area or modeled) before use. The default priors are those of the other families (`vignette("hb-priors")`); a $N(0, 1)$ prior on the coefficients, for example, is set with `inla_engine(control.fixed = list(mean = 0, prec = 1))`. --- ## Practical notes - Give one row per sub-area. The `domain` column identifies the area and is required. - `vardir` must be finite and positive for every sampled sub-area; unsampled sub-areas can have `NA`. Covariates must be complete for all rows, sampled or not. - At least two areas need sampled sub-areas, otherwise $\sigma_v^2$ cannot be estimated. - `eblup_twofold()` has `self_benchmark = TRUE` for estimates that add up to a known total (see `?eblup_twofold`). ## References - Torabi, M., & Rao, J. N. K. (2014). On small area estimation under a sub-area level model. *Journal of Multivariate Analysis*, 127, 36–55. - Chen, C., Wakefield, J., & Lumley, T. (2014). The use of sampling weights in Bayesian hierarchical models for small area estimation. *Spatial and Spatio-temporal Epidemiology*, 11, 33–43. - Franco, C., & Bell, W. R. (2015). Borrowing information over time in binomial/logit normal models for small area estimation. *Statistics in Transition new series*, 16(4), 563–584. - Rao, J. N. K., & Molina, I. (2015). *Small Area Estimation* (2nd ed.). John Wiley & Sons. - Rue, H., Martino, S., & Chopin, N. (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. *Journal of the Royal Statistical Society: Series B*, 71(2), 319–392.