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.
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 yi,j=ai+εi,j for i,j=1,…,10, where aiiid∼N(0,σ2a) is a group effect and εiiid∼N(0,σ2ε) is a residual effect. We define the variance proportion ωaa+ε=σ2aσ2a+σ2ε. Then we denote the different PC prior distributions as: * σ∗∼PC0(U,α), with Prob(σ∗>U)=α, and shrinkage towards σ∗=0. * ωaa+ε∼PC0(m) with Prob(ωaa+ε>m)=0.5 so that m defines the median, and shrinkage towards ωaa+ε=0, i.e., the base model is a model with only εε. * ωaa+ε∼PC1(m) with Prob(ωaa+ε>m)=0.5 so that m defines the median, and shrinkage towards ωaa+ε=1, i.e., the base model is a model with only aa. * ωaa+ε∼PCM(m,c) with Prob(ωaa+ε>m)=0.5 and Prob(logit(1/4)<logit(ωaa+ε)−logit(m)<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 ωaa+ε=m, i.e., the base model is a combination of the effects aa and εε.
Note that PC1(m) on ωaa+ε is equivalent to PC0(1−m) on 1−ωaa+ε=ωεa+ε.
The priors listed above are denoted pc
,
pc0
. pc1
, and pcM
in
makemyprior
.
Consider the hierarchical model for the n=m⋅p observations yi,j, i=1,…p and j=1,…,m, given by yi,j|ηi,j,σ2ε∼N(ηi,j,σ2ε),ηi,j=μ+xiβ+ai+bj, where μ is an intercept, xi is a covariate with coefficient β, and a1,a2,…,apiid∼N(0,σ2a) and b1,b2,…,bmiid∼N(0,σ2b) are random effects. The residuals ε1,ε2,…,εn∼N(0,σ2ε).
First we specify our model by making a formula object (see
?mc
):
Then we put our data in a list
(a
data.frame
can also be used). We simulate the data
here.
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)
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:
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:
ωaa+b∼PCM(0.7,0.5),ωa+ba+b+ε∼PC0(0.25),andσa+b+ε∼PC0(3,0.05).
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)
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.
First, we look at inference with Stan. We must start by compiling the Stan-code:
Then we can do the inference:
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:
Inference with INLA is carried out in a similar way:
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.
See vignette("plotting", package = "makemyprior")
for
more details on functions for plotting.
See:
for details on default settings.
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)
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)
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