Default Priors of the Hierarchical Bayes Models

The models at a glance

The three functions fit generalised linear mixed models. Written out with their default priors, where \(s\) is the scale of the standard deviations: \(s = \mathrm{sd}(y)\) for Gaussian models (the standard deviation of the direct estimates, or of the unit responses for hb_unit()) and \(s = 1\) on the logit or log scale of the other families.

Area level, hb_area(). For area \(d\) (and period \(t\) in a panel),

\[y_d \mid \mu_d \sim f(\mu_d, \psi_d), \qquad g(\mu_d) = x_d^\top \beta + u_d \; [+\, \gamma_t + \delta_{dt}] \; [+\, \log E_d],\]

with the likelihood \(f\), its dispersion \(\psi_d\) and the link \(g\) set by family (see Priors of the likelihood), and

Unit level, hb_unit(). For unit \(j\) in domain \(d\),

\[y_{dj} \mid \mu_{dj} \sim f(\mu_{dj}), \qquad g(\mu_{dj}) = x_{dj}^\top \beta + u_d,\]

with \(f\) normal (\(e_{dj} \sim N(0, \sigma_e^2)\), \(\sigma_e \sim\) half-normal\((s)\)), Bernoulli or Poisson, \(u_d\) iid, Besag or BYM2 with the priors above, and the target the population mean of domain \(d\), computed from the population means of the covariates in Xpop.

Sub-area level, hb_twofold(). For sub-area \(j\) in area \(i\),

\[y_{ij} \mid \mu_{ij} \sim f(\mu_{ij}), \qquad g(\mu_{ij}) = x_{ij}^\top \beta + v_i + u_{ij},\]

with \(f\) normal with the known sampling variance \(V_{ij}\), binomial or Poisson, \(v_i\) iid, Besag or BYM2 and \(u_{ij}\) iid; both SDs have half-normal\((s)\) priors, with \(s\) from the sub-area direct estimates. The area targets aggregate the sub-area means with the weights in weight.

All intrinsic models (Besag, BYM, BYM2 and the random walks) are scaled (scale_model = TRUE) so that \(\sigma\) is a typical marginal standard deviation whatever the graph, which is what makes one prior scale sensible for all of them (Sørbye and Rue 2014).

Priors of the model parameters

Here \(\mathrm{sd}(y)\) is the standard deviation of the observed direct estimates (of the unit responses for hb_unit()). With benchmark_method = "wfq" the coefficients get nearly flat priors (\(N(0, 10^6)\), intercept \(N(0, 10^4)\)) instead, so that the benchmarking covariate is not shrunk. inla_engine() writes the half-normal prior as an expression: prior on the log precision \(\theta = -2 \log \sigma\), so that INLA and NumPyro use the same density.

jax_engine() has two priors that INLA does not have. With robust = TRUE the random effects are Student-\(t\) with degrees of freedom \(\nu \sim\) Uniform\((2.5, 30)\). With prior_beta = "horseshoe" the coefficients get a regularised horseshoe prior (Piironen and Vehtari 2017): local scales \(\lambda_j \sim C^+(1)\), global scale \(\tau \sim C^+(1)\) and slab \(c^2 \sim\) Inv-Gamma\((2, 8)\).

Priors of the likelihood

Why half-normal standard deviations

INLA’s penalised complexity (PC) prior on an SD, \(P(\sigma > U) = \alpha\), is a good default when \(U\) is chosen for the problem (Simpson et al. 2017). The common choice \(U = 1\), \(\alpha = 0.01\) is not: it is tied to the unit of the response. On the mys data, the REML estimate of the area SD is 1.62 and \(\mathrm{sd}(y) = 2.21\), but PC\((1, 0.01)\) gives \(P(\sigma > 1.62) = 0.0006\), and the prior shrinks the area effects far more than the data suggest. In a Fay-Herriot simulation (30 areas) the 95% intervals under PC\((1, 0.01)\) covered 94%, 85% and 38% of the true means with the same data in three units (y times 0.1, 1 and 10).

Scaling the prior to the response removes the dependence on the unit. Among the scaled priors, the half-normal covered as the EBLUP intervals, while the PC priors, which put much mass near \(\sigma = 0\), still shrank the area SD (benchmarks/validation/default_sd_prior.R; 150 replicates, known vardir):

For models with a link function the random effects are on the logit or log scale, where an SD of 1 is already large, so the scale is 1. Against PC\((1, 0.01)\) and half-normal\((2.5)\) (the default of tipsae), half-normal\((1)\) covered at least as well, and it estimated the area SD with less bias (default_sd_prior_glm.R, default_sd_prior_glm_more.R; mean coverage over the scenarios, worst scenario in brackets):

The half-normal prior is also the weakly informative choice of Gelman (2006) for hierarchical SDs. It suits NUTS better than INLA’s vague log-gamma prior on the precision, whose ridge towards \(\sigma = 0\) slows the sampler down.

The BYM2 and AR(1) priors are INLA’s defaults (Riebler et al. 2016), reproduced in JAX: jax_engine() evaluates the PC density of \(\phi\) from the eigenvalues of the scaled ICAR precision, as INLA:::inla.pc.bym.phi() does. With these priors the two engines agree on mys within 0.03 (iid) and 0.05 (BYM2) posterior SD.

Dispersion with one observation per area

In an area-level model each area contributes one direct estimate. For the gamma and negative binomial likelihoods the variance of that estimate depends on a dispersion parameter as well as on the area effect:

• gamma: \(\mathrm{CV}^2 \approx 1/\text{shape} + \sigma^2\);
• negative binomial: \(\mathrm{Var} \approx \mu + \mu^2/\text{size} + \mu^2(e^{\sigma^2} - 1)\).

With one observation per area the data identify only the sum, so the dispersion follows its prior. INLA’s own defaults favour little sampling noise: log-gamma\((1, 0.01)\) on the gamma precision and pc.mgamma\((7)\) on the size, which shrinks towards Poisson. They put the estimate at 76-84 (gamma) and 110-180 (size) whatever the truth, and the intervals then missed the true means when the sampling noise was real. Both engines therefore use a half-normal\((10)\) prior on the gamma precision and on the size. Coverage of 95% intervals, sampled areas, area SD 0.3 / 0.8 (gamma_shape_prior.R, nb_size_prior.R; 30 areas, 100 replicates):

The half-normal\((10)\) prior has a cost: when the sampling noise is small (shape or size 100) its intervals are too wide and its errors larger. hb_area() therefore warns when a gamma or negative binomial model is fitted without vardir. The sampling variances fix the dispersion of each area and give 92-94% (gamma) and 91-95% (negative binomial) coverage (dispersion_vardir.R). A Poisson model is not an alternative for overdispersed counts: it ignores the overdispersion and covered 68-75%.

The Gaussian and beta likelihoods have the same problem without vardir. For the Gaussian residual SD, INLA’s own default, log-gamma\((1, 5 \times 10^{-5})\) on the precision, is disastrous here: it put the residual SD near zero (precision about \(1.5 \times 10^4\)), the area means collapsed onto the direct estimates and their intervals covered almost none of the true means. The half-normal prior with the scale of the response, the prior of the area SD, covered nominally in every unit of the response (likelihood_priors.R; 30 areas, 60 replicates, coverage of 95% intervals for sampled areas, \(y\) in its unit / times 10):

For the beta precision, INLA’s default log-gamma\((1, 0.1)\) covered as well as half-normal\((10)\) or better (area SD 0.3 / 0.8):

so both engines use it. Neither prior covers well when the sampling noise is small (\(\phi = 200\)), and the areas without a direct estimate then covered only 73-80%; give vardir or trials when they are known.

Checking the priors of a fit

The priors that INLA actually used, including any given in inla_engine(), are listed by INLA::inla.priors.used(). The half-normal priors appear as expression: priors on the log precision, with their scale \(s\) written out:

fit <- hb_area(y ~ x1 + x2 + x3, data = mys, vardir = "vardir")
INLA::inla.priors.used(fit$fit)
#> section=[random]
#>   tag=[..domain_id..] component=[..domain_id..]
#>     theta1:
#>       parameter=[log precision]
#>       prior=[expression: s = 2.20822220299071; sigma = exp(-theta/2); ...]

print(fit$engine_spec) shows the engine settings, and for jax_engine() the prior set ("inla" or "weakly_informative").

When to change a prior

The defaults are meant to be safe without tuning, not to encode knowledge of a particular survey. Consider changing them when

• you have real prior knowledge of the variation between areas. A PC prior \(P(\sigma_u > U) = \alpha\) with \(U\) the largest plausible area SD in the unit of the response (or on the link scale) states this directly; see the first example below.
• the areas differ a lot on the link scale. Half-normal\((1)\) gives \(P(\sigma_u > 2) = 0.046\); if logit- or log-scale area effects with an SD above 2 are plausible, raise the scale.
• a dispersion parameter follows its prior. Without vardir (or trials) the dispersion of the gamma, negative binomial and beta likelihoods is weakly identified, and the intervals depend on its prior. Giving the sampling variances is better than tuning that prior.

With few areas, refit with a second prior and compare the estimates in df_hb: if they move by a sizeable fraction of their posterior SD, report the sensitivity.

Changing the priors

inla_engine() takes INLA prior specifications for every hyperparameter. A prior given there replaces the default.

library(fastsae)
data("mys")

# PC prior on the area SD instead of the half-normal default
fit <- hb_area(y ~ x1 + x2 + x3,
  data = mys, vardir = "vardir",
  engine = inla_engine(prior_prec = list(prior = "pc.prec", param = c(1, 0.01)))
)

# INLA's own prior on the negative binomial size (d: counts y, covariate x, exposure E)
fit_nb <- hb_area(y ~ x,
  data = d, family = "nbinomial", exposure = "E",
  engine = inla_engine(control.family = list(hyper = list(theta = list(prior = "pc.mgamma", param = 7))))
)

The other arguments are prior_phi (BYM2), prior_rho (Leroux, SAR, SLM), prior_prec_time and prior_rho_time (temporal effect) and prior_prec_subarea (hb_twofold()). jax_engine() offers the two prior sets of the first table:

fit_jax <- hb_area(y ~ x1 + x2 + x3,
  data = mys, vardir = "vardir",
  engine = jax_engine(priors = "weakly_informative")
)

References

• Gelman, A. (2006). Prior distributions for variance parameters in hierarchical models. Bayesian Analysis, 1(3), 515-534.
• Janicki, R. (2020). Properties of the beta regression model for small area estimation of proportions and application to estimation of poverty rates. Communications in Statistics - Theory and Methods, 49(9), 2264-2284.
• Piironen, J. and Vehtari, A. (2017). Sparsity information and regularization in the horseshoe and other shrinkage priors. Electronic Journal of Statistics, 11(2), 5018-5051.
• Riebler, A., Sørbye, S. H., Simpson, D. and Rue, H. (2016). An intuitive Bayesian spatial model for disease mapping that accounts for scaling. Statistical Methods in Medical Research, 25(4), 1145-1165.
• Simpson, D., Rue, H., Riebler, A., Martins, T. G. and Sørbye, S. H. (2017). Penalising model component complexity: a principled approach to constructing priors. Statistical Science, 32(1), 1-28.
• Sørbye, S. H. and Rue, H. (2014). Scaling intrinsic Gaussian Markov random field priors in spatial modelling. Spatial Statistics, 8, 39-51.