Example i.i.d. model

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

library(makemyprior)

In this vignette, we show how to use makemyprior. This is a package for easy and transparent prior construction, and uses the hierarchical decomposition framework of Fuglstad et al. (2020). We do not go into details on the framework and refer to Fuglstad et al. (2020) for details about the HD prior.

Note that you can also use the standard way of specifying priors component-wise on individual variance components, we show this below.

Prior options

Before going into the details on how to use the software, we give you a short introduction to the HD prior. We also refer to Fuglstad et al. (2020) for details.

We use the penalized complexity (PC) prior (Simpson et al. 2017) to induce shrinkage. This can make a robust prior that stabilizes the inference. We do not go into details on the PC prior here, but the following priors are available in makemyprior:

Consider a random intercept model \(y_{i,j} = a_i + \varepsilon_{i,j}\) for \(i,j = 1, \dots, 10\), where \(a_i \overset{\text{iid}}{\sim} \mathcal{N}(0, \sigma_{\mathrm{a}}^2)\) is a group effect and \(\varepsilon_i \overset{\text{iid}}{\sim} \mathcal{N}(0, \sigma_{\varepsilon}^2)\) is a residual effect. We define the variance proportion \(\omega_{\frac{\mathrm{a}}{\mathrm{a+\varepsilon}}} = \frac{\sigma_{\mathrm{a}}^2}{\sigma_{\mathrm{a}}^2 + \sigma_{\varepsilon}^2}\). Then we denote the different PC prior distributions as: * \(\sigma_{\mathrm{*}} \sim \mathrm{PC}_{\mathrm{0}}(U, \alpha)\), with \(\mathrm{Prob}(\sigma_{\mathrm{*}} > U) = \alpha\), and shrinkage towards \(\sigma_{\mathrm{*}} = 0\). * \(\omega_{\frac{\mathrm{a}}{\mathrm{a+\varepsilon}}} \sim \mathrm{PC}_{\mathrm{0}}(m)\) with \(\mathrm{Prob}(\omega_{\frac{\mathrm{a}}{\mathrm{a+\varepsilon}}} > m) = 0.5\) so that \(m\) defines the median, and shrinkage towards \(\omega_{\frac{\mathrm{a}}{\mathrm{a+\varepsilon}}} = 0\), i.e., the base model is a model with only \(\bm{\varepsilon}\). * \(\omega_{\frac{\mathrm{a}}{\mathrm{a+\varepsilon}}} \sim \mathrm{PC}_{\mathrm{1}}(m)\) with \(\mathrm{Prob}(\omega_{\frac{\mathrm{a}}{\mathrm{a+\varepsilon}}} > m) = 0.5\) so that \(m\) defines the median, and shrinkage towards \(\omega_{\frac{\mathrm{a}}{\mathrm{a+\varepsilon}}} = 1\), i.e., the base model is a model with only \(\bm{a}\). * \(\omega_{\frac{\mathrm{a}}{\mathrm{a+\varepsilon}}} \sim \mathrm{PC}_{\mathrm{M}}(m, c)\) with \(\mathrm{Prob}(\omega_{\frac{\mathrm{a}}{\mathrm{a+\varepsilon}}} > m) = 0.5\) and \(\mathrm{Prob}(\mathrm{logit}(1/4) < \mathrm{logit}(\omega_{\frac{\mathrm{a}}{\mathrm{a+\varepsilon}}}) - \mathrm{logit}(m) < \mathrm{logit}(3/4)) = c\) so that \(m\) defines the median, and \(c\) says something about how concentrated the distribution is around the median. The shrinkage is towards \(\omega_{\frac{\mathrm{a}}{\mathrm{a+\varepsilon}}} = m\), i.e., the base model is a combination of the effects \(\bm{a}\) and \(\bm{\varepsilon}\).

Note that \(\mathrm{PC}_{\mathrm{1}}(m)\) on \(\omega_{\frac{\mathrm{a}}{\mathrm{a+\varepsilon}}}\) is equivalent to \(\mathrm{PC}_{\mathrm{0}}(1-m)\) on \(1-\omega_{\frac{\mathrm{a}}{\mathrm{a+\varepsilon}}} = \omega_{\frac{\mathrm{\varepsilon}}{\mathrm{a+\varepsilon}}}\).

The priors listed above are denoted pc, pc0. pc1, and pcM in makemyprior.

Model description

Consider the hierarchical model for the \(n = m \cdot p\) observations \(y_{i,j}\), \(i = 1, \ldots p\) and \(j = 1, \ldots, m\), given by \[\begin{align*} y_{i,j}|\eta_{i,j}, \sigma_{\varepsilon}^2 &\sim \mathcal{N}(\eta_{i,j}, \sigma_{\varepsilon}^2), \\ \eta_{i,j} &= \mu + x_i \beta + a_i + b_j, \end{align*}\] where \(\mu\) is an intercept, \(x_i\) is a covariate with coefficient \(\beta\), and \(a_1, a_2, \ldots, a_p \overset{\text{iid}}{\sim} \mathcal{N}(0, \sigma_\mathrm{a}^2)\) and \(b_1, b_2, \ldots, b_m \overset{\text{iid}}{\sim} \mathcal{N}(0, \sigma_\mathrm{b}^2)\) are random effects. The residuals \(\varepsilon_1, \varepsilon_2, \dots, \varepsilon_n \sim \mathcal{N}(0, \sigma_{\varepsilon}^2)\).

Make data and linear predictor

First we specify our model by making a formula object (see ?mc):

formula <- y ~ x + mc(a) + mc(b)

Then we put our data in a list (a data.frame can also be used). We simulate the data here.

p <- 10
m <- 10
n <- m*p

set.seed(1)
data <- list(a = rep(1:p, each = m),
             b = rep(1:m, times = p),
             x = runif(n))

data$y <- data$x + rnorm(p, 0, 0.5)[data$a] +
  rnorm(m, 0, 0.3)[data$b] + rnorm(n, 0, 1)

Make prior

Then we make the prior object using the function make_prior. It needs the arguments we made above, formula and data, a likelihood family (Gaussian likelihood is the default), and optional priors for intercept and covariate coefficients (both have a Gaussian distribution with \(0\) mean and a standard deviation of \(1000\)). Note that the observations y are not used to create the prior, but is included in the prior object as all the information about the inference is stored there.

prior <- make_prior(formula, data, family = "gaussian",
                    intercept_prior = c(0, 1000),
                    covariate_prior = list(x = c(0, 100)))
#> Warning: Did not find a tree, using default tree structure instead.

This gives the default prior, which is a prior where all model effects are assigned an equal amount of variance through a symmetric Dirichlet distribution. The default prior on the total variance depends on the likelihood. See Section default settings for details on default settings.

We print details about the prior, plot the prior to see how the distributions look, and plot the prior tree structure:

summary(prior)
#> Model: y ~ x + mc(a) + mc(b)
#> Tree structure: a_b_eps = (a,b,eps)
#> 
#> Weight priors:
#>  (w[a/a_b_eps], w[b/a_b_eps]) ~ Dirichlet(3)
#> Total variance priors:
#>  V[a_b_eps] ~ Jeffreys'
#> 
#> Covariate priors: intercept ~ N(0, 1000^2), x ~ N(0, 100^2)

plot_prior(prior) # or plot(prior)

plot_tree_structure(prior)

Now we can use a graphical interface to choose our prior. We do not show this in the vignette, but it can be opened with the following command:

new_prior <- makemyprior_gui(prior)

The output (which we store in new_prior) is of the same class as the output from make_prior, and can be used directly for inference.

With the following command, we specify this prior:

\[\begin{equation} \omega_{\frac{\mathrm{a}}{\mathrm{a+b}}} \sim \mathrm{PC}_{\mathrm{M}}(0.7, 0.5),\, \omega_{\frac{\mathrm{a+b}}{\mathrm{a+b} + \varepsilon}} \sim \mathrm{PC}_{\mathrm{0}}(0.25),\,\text{and}\, \sigma_{\mathrm{a+b} + \varepsilon} \sim \mathrm{PC}_{\mathrm{0}}(3, 0.05). \label{eq:software:examplemodel_prior} \end{equation}\]


new_prior <- make_prior(
  formula, data,
  prior = list(
    tree = "s1 = (a, b); s2 = (s1, eps)",
    w = list(s1 = list(prior = "pcM", param = c(0.7, 0.5)),
             s2 = list(prior = "pc1", param = 0.75)),
    V = list(s2 = list(prior = "pc0", param = c(3, 0.05)))
    ),
    covariate_prior = list(x = c(0, 100))
)

summary(new_prior)
#> Model: y ~ x + mc(a) + mc(b)
#> Tree structure: a_b = (a,b); eps_a_b = (eps,a_b)
#> 
#> Weight priors:
#>  w[a/a_b] ~ PCM(0.7, 0.5)
#>  w[eps/eps_a_b] ~ PC0(0.25)
#> Total variance priors:
#>  sqrt(V)[eps_a_b] ~ PC0(3, 0.05)
#> 
#> Covariate priors: intercept ~ N(0, 1000^2), x ~ N(0, 100^2)
plot_prior(new_prior)

plot_tree_structure(new_prior)

Inference

We can carry out inference with Stan (Carpenter et al. 2017) and INLA (Rue, Martino, and Chopin 2009). Note that we in this vignette do not run the inference, as it takes time and will slow down the compilation of the vignette and thus the package download, but the code is included below and the user can carry out the inference with that.

Stan

First, we look at inference with Stan. We must start by compiling the Stan-code:

compile_stan(save = TRUE)

Then we can do the inference:

posterior1 <- inference_stan(new_prior, iter = 1e4, chains = 1, seed = 1)

We can look at the graphs of the posterior:

plot_posterior_stan(posterior1, param = "prior", prior = TRUE) # on the scale of the prior, together with the prior
plot_posterior_stan(posterior1, param = "variance") # on variance scale
plot_fixed_posterior(posterior1) # fixed effects

We can also sample from the prior and compare on variance scale:

prior1 <- inference_stan(new_prior, use_likelihood = FALSE, iter = 1e4, chains = 1, seed = 1)
plot_several_posterior_stan(list(Prior = prior1, Posterior = posterior1))

INLA

Inference with INLA is carried out in a similar way:

posterior2 <- inference_inla(new_prior)

And we can look at some posterior diagnostics. Note that we can only look at the posteriors on variance/precision/standard deviation scale when doing inference with INLA.

plot_posterior_variance(posterior2) # on variance scale
plot_fixed_posterior(posterior1)

See vignette("plotting", package = "makemyprior") for more details on functions for plotting.

Default settings

See:

?makemyprior_models

for details on default settings.

Additional examples

We include some additional examples on how to create various prior distributions. We still use the same model and data, and change the joint prior on the variances. We do not run inference. Note that the values of the priors are NOT based on knowledge about the model, but chosen to show the different options of the package. See vignette("wheat_breeding", package = "makemyprior"), vignette("latin_square", package = "makemyprior"), and vignette("neonatal_mortality", package = "makemyprior") for examples where we discuss how expert knowledge can be used to set the priors.


prior2 <- make_prior(formula = formula, data = data,
                     prior = list(tree = "(a); (b); (eps)",
                                  V = list(
                                    a = list(prior = "pc", param = c(1, 0.05)),
                                    b = list(prior = "pc", param = c(2, 0.05)),
                                    eps = list(prior = "pc", param = c(3, 0.05))
                                  )))

plot_prior(prior2)

plot_tree_structure(prior2)

prior3 <- make_prior(formula = formula, data = data,
                     prior = list(tree = "s1 = (a, b); (eps)",
                                  V = list(
                                    s1 = list(prior = "pc", param = c(3, 0.05)),
                                    eps = list(prior = "pc", param = c(3, 0.05))),
                                  w = list(
                                    s1 = list(prior = "pcM", param = c(0.5, 0.8))
                                  )
                                  ))

plot_prior(prior3)

plot_tree_structure(prior3)
sessionInfo()
#> R version 4.1.0 (2021-05-18)
#> Platform: aarch64-apple-darwin20 (64-bit)
#> Running under: macOS Big Sur 11.3.1
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/lib/libRblas.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/lib/libRlapack.dylib
#> 
#> locale:
#> [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] makemyprior_1.1.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] tidyselect_1.1.1  xfun_0.30         bslib_0.3.1       shinyjs_2.1.0    
#>  [5] purrr_0.3.4       splines_4.1.0     lattice_0.20-44   colorspace_2.0-3 
#>  [9] vctrs_0.3.8       generics_0.1.1    htmltools_0.5.2   yaml_2.3.5       
#> [13] utf8_1.2.2        rlang_1.0.2       jquerylib_0.1.4   pillar_1.7.0     
#> [17] later_1.3.0       glue_1.6.2        lifecycle_1.0.1   stringr_1.4.0    
#> [21] munsell_0.5.0     gtable_0.3.0      htmlwidgets_1.5.4 visNetwork_2.1.0 
#> [25] evaluate_0.15     labeling_0.4.2    knitr_1.37        fastmap_1.1.0    
#> [29] httpuv_1.6.5      fansi_1.0.2       highr_0.9         Rcpp_1.0.8.2     
#> [33] xtable_1.8-4      scales_1.1.1      promises_1.2.0.1  jsonlite_1.8.0   
#> [37] farver_2.1.0      mime_0.12         ggplot2_3.3.5     digest_0.6.29    
#> [41] stringi_1.7.6     dplyr_1.0.7       shiny_1.7.1       grid_4.1.0       
#> [45] cli_3.2.0         tools_4.1.0       magrittr_2.0.2    sass_0.4.0       
#> [49] tibble_3.1.6      crayon_1.5.0      pkgconfig_2.0.3   ellipsis_0.3.2   
#> [53] MASS_7.3-54       Matrix_1.3-3      rmarkdown_2.13    rstudioapi_0.13  
#> [57] R6_2.5.1          compiler_4.1.0

References

Carpenter, Bob, Andrew Gelman, Matthew D Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. 2017. “"Stan: A probabilistic programming language".” "Journal of Statistical Software" 76 (1).
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.
Rue, Håvard, Sara Martino, and Nicolas Chopin. 2009. “"Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations".” "Journal of the Royal Statistical Society: Series b (Statistical Methodology)" 71 (2): 319–92.
Simpson, Daniel, Håvard Rue, Andrea Riebler, Thiago G. Martins, and Sigrunn H. Sørbye. 2017. “Penalising model component complexity: a principled, practical approach to constructing priors.” Statistical Science 32 (1): 1–28.