We carry out a similar study of neonatal mortality in Kenya as one by Fuglstad et al. (2020). We model the neonatal mortality, defined as the number of deaths if infants the first month of live per birth. We use the linear predictor: \[ \eta_{i,j} = \mathrm{logit}(p_{i,j}) = \mu + x_{i,j} \beta + u_i + v_i + \nu_{i,j}, \ i = 1, \dots, n, \ j = 1, \dots, m_i, \] and use \(y_{i,j} | b_{i,j}, p_{i,j} \sim \mathrm{Binomial}(b_{i,j}, p_{i,j})\), for cluster \(j\) in county \(i\). We have between $m_i {6, 7, 8} clusters in each of the \(n = 47\) counties (see e.g. Fuglstad et al. (2020) for a map of the counties).

- \(b_{i,j}\) is the number of live births,
- \(y_{i,j}\) is the number of neonatal deaths,
- \(\mu\) is an intercept with a \(N(0, 1000^2)\) prior,
- \(\beta\) is a coefficient with a \(N(0, 1000^2)\) prior for \(x_{i,j}\), which is an indicator classifying cluster \(j\) in county \(i\) as urban (\(x_{i,j} = 1\)) or rural (\(x_{i,j} = 0\)),
- \(\nu_i \sim N_n(0, \sigma_{\nu}^2)\) is an i.i.d. random effect for cluster
- \(v_i \sim N_n(0, \sigma_v^2)\) is an i.i.d. random effect for county, and
- \(\mathbf{u}\) is a Besag effect on county with variance \(\sigma_u^2\) and a sum-to-zero constraint.

We need a neighborhood graph for the counties, which is found in
`makemyprior`

. We scale the Besag effect to have a
generalized variance equal to \(1\).

```
# neighborhood graph
graph_path <- paste0(path.package("makemyprior"), "/neonatal.graph")
formula <- y ~ urban + mc(nu) + mc(v) +
mc(u, model = "besag", graph = graph_path, scale.model = TRUE)
```

We use the dataset `neonatal_mortality`

in
`makemyprior`

, and present three priors. We do not carry out
inference, as it takes time and will slow down the compilation of the
vignettes by a lot, but include code so the user can run the inference
themselves.

We prefer coarser over finer unstructured effects, and unstructured
over structured effects. That means that we prefer \(\mathbf{v}\) over \(\mathbf{u}\) and \(\mathbf{v} + \mathbf{u}\) over \(\mathbf{\nu}\) in the prior. We achieve
this with a prior that distributes the county variance with shrinkage
towards the unstructured county effect, and the total variance towards
the county effects. Following (**fuglstad?**), we induce
shrinkage on the total variance such that we have a 90% credible
interval of \((0.1, 10)\) for the
effect of \(\exp(v_i + u_i +
\nu_{i,j})\). We use the function
`find_pc_prior_param`

in `makemyprior`

to find the
parameters for the PC prior:

```
set.seed(1)
find_pc_prior_param(lower = 0.1, upper = 10, prob = 0.9, N = 2e5)
#> U = 3.353132
#> Prob(0.09866969 < exp(eta) < 9.892902) = 0.9
```

```
prior1 <- make_prior(
formula, neonatal_data, family = "binomial",
prior = list(tree = "s1 = (u, v); s2 = (s1, nu)",
w = list(s1 = list(prior = "pc0", param = 0.25),
s2 = list(prior = "pc1", param = 0.75)),
V = list(s2 = list(prior = "pc",
param = c(3.35, 0.05)))))
prior1
#> Model: y ~ urban + mc(nu) + mc(v) + mc(u, model = "besag", graph = graph_path,
#> scale.model = TRUE)
#> Tree structure: v_u = (v,u); nu_v_u = (nu,v_u)
#>
#> Weight priors:
#> w[v/v_u] ~ PC1(0.75)
#> w[nu/nu_v_u] ~ PC0(0.25)
#> Total variance priors:
#> sqrt(V)[nu_v_u] ~ PC0(3.35, 0.05)
```

Inference can be carried out by running:

```
posterior1 <- inference_stan(prior1, iter = 15000, warmup = 5000,
seed = 1, init = "0", chains = 1)
plot_posterior_stan(posterior1, param = "prior", plot_prior = TRUE)
```

For inference with INLA:

```
posterior1_inla <- inference_inla(prior1, Ntrials = neonatal_data$Ntrials)
plot_posterior_stdev(posterior1_inla)
```

Note the `Ntrials`

argument fed to
`inference_inla`

.

We use a prior without any knowledge, and use the default prior:

```
prior2 <- make_prior(formula, neonatal_data, family = "binomial")
#> Warning: Did not find a tree, using default tree structure instead.
prior2
#> Model: y ~ urban + mc(nu) + mc(v) + mc(u, model = "besag", graph = graph_path,
#> scale.model = TRUE)
#> Tree structure: nu_v_u = (nu,v,u)
#>
#> Weight priors:
#> (w[nu/nu_v_u], w[v/nu_v_u]) ~ Dirichlet(3)
#> Total variance priors:
#> sqrt(V)[nu_v_u] ~ PC0(1.6, 0.05)
```

Inference can be carried out by running:

```
posterior2 <- inference_stan(prior2, iter = 15000, warmup = 5000,
seed = 1, init = "0", chains = 1)
plot_posterior_stan(posterior2, param = "prior", plot_prior = TRUE)
```

```
sessionInfo()
#> R version 4.3.2 (2023-10-31)
#> Platform: aarch64-apple-darwin20 (64-bit)
#> Running under: macOS Sonoma 14.1.1
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
#>
#> locale:
#> [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> time zone: Europe/Oslo
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] makemyprior_1.2.2
#>
#> loaded via a namespace (and not attached):
#> [1] Matrix_1.6-1.1 gtable_0.3.4 jsonlite_1.8.8 highr_0.10
#> [5] dplyr_1.1.4 compiler_4.3.2 promises_1.2.1 tidyselect_1.2.0
#> [9] Rcpp_1.0.12 later_1.3.2 jquerylib_0.1.4 splines_4.3.2
#> [13] scales_1.3.0 yaml_2.3.8 fastmap_1.1.1 mime_0.12
#> [17] lattice_0.21-9 ggplot2_3.4.4 R6_2.5.1 labeling_0.4.3
#> [21] shinyjs_2.1.0 generics_0.1.3 knitr_1.45 htmlwidgets_1.6.4
#> [25] visNetwork_2.1.2 MASS_7.3-60 tibble_3.2.1 munsell_0.5.0
#> [29] shiny_1.8.0 bslib_0.6.1 pillar_1.9.0 rlang_1.1.3
#> [33] utf8_1.2.4 cachem_1.0.8 httpuv_1.6.14 xfun_0.42
#> [37] sass_0.4.8 cli_3.6.2 withr_3.0.0 magrittr_2.0.3
#> [41] digest_0.6.34 grid_4.3.2 xtable_1.8-4 rstudioapi_0.14
#> [45] lifecycle_1.0.4 vctrs_0.6.5 evaluate_0.23 glue_1.7.0
#> [49] farver_2.1.1 fansi_1.0.6 colorspace_2.1-0 rmarkdown_2.25
#> [53] ellipsis_0.3.2 tools_4.3.2 pkgconfig_2.0.3 htmltools_0.5.7
```

Fuglstad, Geir-Arne, Ingeborg Gullikstad Hem, Alexander Knight, Håvard
Rue, and Andrea Riebler. 2020. “Intuitive
Joint Priors for Variance Parameters.” *Bayesian
Anal.* 15 (4): 1109–37.