---
title: "Default Priors of the Hierarchical Bayes Models"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Default Priors of the Hierarchical Bayes Models}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", eval = FALSE)
```

`hb_area()`, `hb_unit()` and `hb_twofold()` fit their models with INLA (`inla_engine()`,
the default) or with NUTS in NumPyro (`jax_engine()`). Both engines use the same default
priors, so a model gives the same estimates whichever engine fits it. This article lists
those priors, explains why they were chosen and summarises the simulations behind them.

With few areas, or with one observation per area, the prior of a variance or dispersion
parameter is not a technicality: the data hardly inform it, and its prior decides how much
the estimates are shrunk and how wide the intervals are.

## 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](#priors-of-the-likelihood)), and

| Term | Argument | Model | Default prior |
|---|---|---|---|
| $\beta$ | `formula` | fixed effects | intercept flat; others $N(0, 1000)$ |
| $u_d$ | `spatial = "none"` | iid $N(0, \sigma_u^2)$ | $\sigma_u \sim$ half-normal$(s)$ |
| | `"besag"` | intrinsic CAR (ICAR) on the graph of `W`, scaled | $\sigma_u \sim$ half-normal$(s)$ |
| | `"bym"` | iid + ICAR, each with its own SD | both SDs half-normal$(s)$ |
| | `"bym2"` | $\sigma_u(\sqrt{1 - \phi}\, v_d + \sqrt{\phi}\, u^*_d)$, $u^*$ a scaled ICAR | $\sigma_u \sim$ half-normal$(s)$; PC prior $P(\phi < 0.5) = 0.5$ |
| | `"leroux"` | precision $\sigma_u^{-2}\{(1 - \lambda) I + \lambda R\}$, $R$ the ICAR structure | $\sigma_u \sim$ half-normal$(s)$; $\lambda \sim$ Uniform$(0, 1)$ |
| | `"sar"` | $u = (I - \rho W)^{-1} \varepsilon$, $\varepsilon_d \sim N(0, \sigma_u^2)$ | $\sigma_u \sim$ half-normal$(s)$; $\rho \sim$ Uniform$(-0.999, 0.999)$ |
| | `"slm"` | $\eta = (I - \rho W)^{-1}(X\beta + \varepsilon)$ | $\sigma_u \sim$ half-normal$(s)$; $\rho$ uniform on the range allowed by the eigenvalues of `W` |
| $\gamma_t$ | `temporal` | `"rw1"`, `"rw2"`, `"ar1"` or `"iid"` | $\sigma_\gamma \sim$ half-normal$(s)$; AR(1): $N(0, 1/0.15)$ on $\log\{(1 + \rho)/(1 - \rho)\}$ |
| $\delta_{dt}$ | `st_interaction` | domain-specific trends, a separable field or Knorr-Held types I-IV | its SD half-normal$(s)$; AR(1) correlation as for $\gamma_t$ |

**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

| Parameter | Default (`inla_engine()`, `jax_engine(priors = "inla")`) | `jax_engine(priors = "weakly_informative")` |
|---|---|---|
| Random-effect SD $\sigma$ (area, sub-area, time, area-by-time) | Half-normal$(s)$: $s = \mathrm{sd}(y)$ for Gaussian models, $s = 1$ on the link scale otherwise | Half-normal$(1)$ on the standardised response or the link scale |
| BYM2 mixing parameter $\phi$ | PC prior with $P(\phi < 0.5) = 0.5$ | Uniform$(0, 1)$ |
| AR(1) correlation $\rho$ | $N(0, 1/0.15)$ on $\log\{(1 + \rho)/(1 - \rho)\}$ | Uniform$(-1, 1)$ |
| Leroux $\lambda$, SAR $\rho$ | Uniform | Uniform |
| Coefficients $\beta$ | $N(0, 1000)$; intercept flat (INLA) or nearly flat (JAX) | $N(0, 10^2)$ on the standardised scale ($N(0, 2.5^2)$ for logit links) |

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

| Parameter | When | Prior |
|---|---|---|
| Sampling variance (Gaussian), gamma precision | `vardir` given | Fixed per area: $1/V_d$; $y_d^2 / V_d$ |
| Beta precision | `vardir` given | Fixed per area at $\hat\mu_d(1 - \hat\mu_d)/V_d - 1$, with $\hat\mu_d$ from a first fit |
| Beta precision | `trials` given (no `vardir`) | Fixed per area at $n_d - 1$ |
| Negative binomial size | `vardir` given | Fixed per area at $\hat\mu_d^2 / (V_d - \hat\mu_d)$, with $\hat\mu_d$ from a first fit |
| Gamma precision | no `vardir` | Half-normal$(10)$, both engines |
| Negative binomial size | no `vardir` | Half-normal$(10)$, both engines |
| Gaussian residual SD | no `vardir` (and the unit errors of `hb_unit()`) | Half-normal$(\mathrm{sd}(y))$, both engines |
| Beta precision $\phi$ | no `vardir`, no `trials` | Log-gamma$(1, 0.1)$ on $\log\phi$, i.e. $\phi \sim$ Exponential$(0.1)$, both engines (half-normal$(10)$ with `jax_engine(priors = "weakly_informative")`) |

## 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`):

| Areas | Area SD | EBLUP | Half-normal$(\mathrm{sd}\,y)$ | PC, $\alpha = 0.01$ | PC, $\alpha = 0.05$ |
|---|---|---|---|---|---|
| 15 | 0.4 | 0.960 | 0.966 | 0.917 | 0.942 |
| 30 | 0.4 | 0.940 | 0.958 | 0.905 | 0.929 |
| 60 | 0.4 | 0.927 | 0.943 | 0.898 | 0.918 |
| 15 | 0.8 | 0.930 | 0.932 | 0.853 | 0.890 |
| 30 | 0.8 | 0.934 | 0.929 | 0.888 | 0.909 |
| 60 | 0.8 | 0.942 | 0.939 | 0.919 | 0.928 |

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):

| Response | Half-normal$(1)$ | Half-normal$(2.5)$ | PC$(1, 0.01)$ |
|---|---|---|---|
| Binomial, $n_d$ 20-200 | 0.927-0.950 | 0.927-0.952 | 0.918-0.944 |
| Binomial, $n_d$ 5-20 | 0.950 (0.932) | 0.953 (0.938) | 0.919 (0.864) |
| Beta, $n_d$ 5-20 | 0.938 (0.925) | 0.940 (0.929) | 0.908 (0.873) |
| Poisson | 0.947 (0.944) | 0.948 (0.944) | 0.943 (0.934) |
| Binomial BYM2 | 0.947 (0.941) | 0.947 (0.942) | 0.942 (0.934) |

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):

| Likelihood | True dispersion | Half-normal$(10)$ | INLA's default |
|---|---|---|---|
| Gamma | shape 5 | 0.930 / 0.925 | 0.593 / 0.522 |
| Gamma | shape 20 | 0.959 / 0.986 | 0.862 / 0.826 |
| Gamma | shape 100 | 0.950 / 0.999 | 0.979 / 0.993 |
| Negative binomial | size 3 | 0.928 / 0.896 | 0.816 / 0.756 |
| Negative binomial | size 10 | 0.949 / 0.965 | 0.913 / 0.915 |
| Negative binomial | size 100 | 0.939 / 0.990 | 0.965 / 0.980 |

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):

| Area SD | Residual SD | Half-normal$(\mathrm{sd}\,y)$ | INLA's default |
|---|---|---|---|
| 0.5 | 0.5 | 0.943 / 0.943 | 0.038 / 0.006 |
| 0.5 | 1 | 0.949 / 0.947 | 0.024 / 0.001 |
| 1 | 0.5 | 0.955 / 0.968 | 0.033 / 0.003 |
| 1 | 1 | 0.953 / 0.938 | 0.019 / 0.002 |

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):

| True $\phi$ | Log-gamma$(1, 0.1)$ | Half-normal$(10)$ |
|---|---|---|
| 10 | 0.960 / 0.937 | 0.965 / 0.938 |
| 50 | 0.931 / 0.903 | 0.913 / 0.869 |
| 200 | 0.875 / 0.925 | 0.841 / 0.871 |

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.

### Beta precision from `vardir`

The beta variance is $\mu_d(1 - \mu_d)/(\phi_d + 1)$, so a known sampling variance fixes
$\phi_d = \mu_d(1 - \mu_d)/V_d - 1$ at the unknown true proportion. Putting the direct
estimate $y_d$ in place of $\mu_d$ (Janicki 2020) gives the model a sampling variance of
$V_d\,\mu_d(1 - \mu_d)/\{y_d(1 - y_d)\}$ instead of $V_d$, too small wherever $y_d$ happens
to fall far from $\mu_d$. `hb_area()` therefore fits the model twice and takes $\phi_d$ at
the posterior mean of the first fit. tipsae (`type_disp = "var"`) instead lets $\phi_d$
follow the model mean inside the MCMC; that keeps the variance at $V_d$ exactly, but a low
mean then also lowers the precision, which gives the likelihood a second mode near 0. Coverage
of 95% intervals for sampled areas, area SD 0.3 / 0.8 (`beta_vardir_phi.R`; 30 areas,
60 replicates, $V_d$ the true sampling variance):

| Proportions near | Direct estimate | Two stages | tipsae | True $\phi_d$ |
|---|---|---|---|---|
| 0.08 | 0.844 / 0.812 | 0.938 / 0.896 | 0.938 / 0.931 | 0.953 / 0.918 |
| 0.23 | 0.922 / 0.931 | 0.948 / 0.933 | 0.951 / 0.954 | 0.950 / 0.944 |

tipsae covered as well, but 11 of its 240 fits did not converge (R-hat up to 8.7).

## 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:

```{r}
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.

```{r}
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:

```{r}
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.