Neonatal mortality

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: \[ \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.

Prior 1

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)

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

Neonatal mortality

Ingeborg Gullikstad Hem (ingeborg.hem@ntnu.no)

Model

Prior 1

Prior 2

References