Neonatal mortality

library(makemyprior)

Model

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: ηi, j = logit(pi, j) = μ + xi, jβ + ui + vi + νi, j, i = 1, …, n, j = 1, …, mi, and use yi, j|bi, j, pi, j ∼ Binomial(bi, j, pi, 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).

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.

Prior 1

We prefer coarser over finer unstructured effects, and unstructured over structured effects. That means that we prefer v over u and v + u over ν 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 (vi + ui + ν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)
plot_prior(prior1) # or plot(prior1)
plot_tree_structure(prior1)

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.

Prior 2

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)
plot_prior(prior2)
plot_tree_structure(prior2)

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.4.2 (2024-10-31)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] makemyprior_1.2.2 rmarkdown_2.29   
#> 
#> loaded via a namespace (and not attached):
#>  [1] sass_0.4.9        utf8_1.2.4        generics_0.1.3    lattice_0.22-6   
#>  [5] digest_0.6.37     magrittr_2.0.3    evaluate_1.0.1    grid_4.4.2       
#>  [9] fastmap_1.2.0     jsonlite_1.8.9    Matrix_1.7-1      promises_1.3.0   
#> [13] fansi_1.0.6       scales_1.3.0      jquerylib_0.1.4   cli_3.6.3        
#> [17] shiny_1.9.1       rlang_1.1.4       visNetwork_2.1.2  munsell_0.5.1    
#> [21] splines_4.4.2     withr_3.0.2       cachem_1.1.0      yaml_2.3.10      
#> [25] tools_4.4.2       dplyr_1.1.4       colorspace_2.1-1  ggplot2_3.5.1    
#> [29] httpuv_1.6.15     buildtools_1.0.0  vctrs_0.6.5       R6_2.5.1         
#> [33] mime_0.12         lifecycle_1.0.4   htmlwidgets_1.6.4 MASS_7.3-61      
#> [37] shinyjs_2.1.0     pkgconfig_2.0.3   pillar_1.9.0      bslib_0.8.0      
#> [41] later_1.3.2       gtable_0.3.6      glue_1.8.0        Rcpp_1.0.13-1    
#> [45] xfun_0.49         tibble_3.2.1      tidyselect_1.2.1  sys_3.4.3        
#> [49] knitr_1.49        farver_2.1.2      xtable_1.8-4      htmltools_0.5.8.1
#> [53] maketools_1.3.1   labeling_0.4.3    compiler_4.4.2

References

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.