#  makemyprior: Intuitive Construction of Joint Priors for Variance Parameters
in R

Ingeborg Gullikstad Hem111Corresponding author<EMAIL_ADDRESS>Department of Mathematical Sciences, NTNU, Norway Geir-Arne Fuglstad
Department of Mathematical Sciences, NTNU, Norway Andrea Riebler Department of
Mathematical Sciences, NTNU, Norway

Abstract

Priors allow us to robustify inference and to incorporate expert knowledge in
Bayesian hierarchical models. This is particularly important when there are
random effects that are hard to identify based on observed data. The challenge
lies in understanding and controlling the joint influence of the priors for
the variance parameters, and makemyprior is an R package that guides the
formulation of joint prior distributions for variance parameters. A joint
prior distribution is constructed based on a hierarchical decomposition of the
total variance in the model along a tree, and takes the entire model structure
into account. Users input their prior beliefs or express ignorance at each
level of the tree. Prior beliefs can be general ideas about reasonable ranges
of variance values and need not be detailed expert knowledge. The constructed
priors lead to robust inference and guarantee proper posteriors. A graphical
user interface facilitates construction and assessment of different choices of
priors through visualization of the tree and joint prior. The package aims to
expand the toolbox of applied researchers and make priors an active component
in their Bayesian workflow.

Keywords: Bayesian hierarchical models, robust inference, joint prior
distributions, hierarchical variance decomposition, graphical user interface,
R.

##  1 Introduction

Bayesian modelling is more available than ever through fast and easy-to-use
software programs for Bayesian inference such as Integrated Nested Laplace
Approximations (INLA, Rue et al., , 2009) in R through INLA (see
www.r-inla.org), Stan (Carpenter et al., , 2017) with the R interface rstan
(Stan Development Team, , 2020) and dependencies such as rstanarm (Goodrich et
al., , 2020), shinystan (Gabry, , 2018) and loo (Vehtari et al., , 2020),
Template Model Builder (TMB, Kristensen et al., , 2016), JAGS (Plummer, ,
2017), Bayesian Analysis Toolkit (BAT, Caldwell et al., , 2009), and more.
These programs offer many ways to construct complex models suitable for a wide
range of applications. Default settings for priors and hyperparameters are
usually given, and instructions on how to change the default settings can be
found in the software documentation. However, even though there is an
increasing focus on that prior distributions should be chosen consciously
(Zondervan-Zwijnenburg et al., , 2017; Gelman et al., , 2020; Smid and Winter,
, 2020), guides on how the priors should be chosen are missing. Our goal with
makemyprior is to close this gap and to increase the awareness of prior
choices. We empower users to actively select priors that are suitable for
their model structure and application at hand. The user is made aware of which
priors are used in the model and what they express. In particular, the
assumptions underlying the default settings are made clear. The package is
available from the Comprehensive R Archive Network (CRAN):
https://CRAN.R-project.org/package=makemyprior.

In Bayesian hierarchical models, the variation in the observations is modelled
through a combination of an observation model, a latent model, and priors for
the parameters. The observation model defines a generative model for the
observed data conditional on the latent model. The link between the latent
model and each observation acts through the transformation of a linear
predictor, which we assume is a linear combination of fixed and random
effects. The goal of the linear predictor is to explain the variation in the
true signal. We concentrate on latent Gaussian models, where the linear
predictor is composed of random model components that follow multivariate
Gaussian distributions conditional on the parameters. These parameters control
the random effects and need prior distributions. We focus on the most central
type of model parameters: the variance parameters.

Parameters controlling means and medians, such as the coefficients of fixed
effects, are close to the data and tolerate vague priors (Goel and Degroot, ,
1981; Gelman et al., , 2020). We give them Gaussian priors with zero mean and
a fixed high variance. The specification of a prior distribution for a
variance parameter is challenging (Lambert et al., , 2005; Gelman et al., ,
2017), but is at the same time a strong feature of Bayesian inference. Prior
knowledge, obtained from previous experiments or comparable investigations,
and expert knowledge can be included to make the model more robust. It is
convenient to use the default priors included in programs for inference,
however, by using defaults we are not fully utilizing the Bayesian framework.
Currently, there are no other R packages that intuitively allows the proper
inclusion of such knowledge and simple visualization of chosen priors in a
straightforward manner.

The makemyprior package applies the hierarchical decomposition (HD) prior
framework proposed by Fuglstad et al.,  (2020), where the variation in the
observed data is distributed to the random effects following a prior tree
structure. In this tree, the leaf nodes represent the random effects specified
in the linear predictor, and the top (root) node represents the sum of the
random effect variance, denoted the total variance. How much variance is
distributed from a parent node to its child nodes is determined by a split in
the tree, and this procedure continues down until the leaf nodes. This gives
us a parameterization with proportions of variance, instead of the more common
variance parameter parameterization. The priors for those variance proportion
parameters can be specified intuitively and transparently as they often
coincide with the scale on which prior or expert knowledge exists, such as in
genomic modelling (Holand et al., , 2013; Hem et al., , 2021) and disease
mapping (Wakefield, , 2006). Through the makemyprior package, we raise the
awareness for prior selection and guide the user to formulate sensible prior
distributions, that are automatically computed and visualized in an intuitive
manner.

Prior distribution can be either specified directly within R or through a
graphical user interface (GUI). In the GUI the user can inspect the prior tree
and adapt it as needed. One can click through the splits independently to
specify the beliefs for each split. A guide asking relevant questions about
the model to ease the prior specification process further is also available.
The user can be ignorant and distribute the variance equally to the child
nodes through a Dirichlet prior, or exploit expert knowledge implemented via
penalized complexity (PC, Simpson et al., , 2017) priors. After completing the
prior specification, inference can be carried out directly with the R packages
rstan and INLA.

We begin with explaining the concepts of total variance and hierarchical
variance decomposition through two motivating examples in Section 2. We then
introduce the necessary background in Section 3 before we present the
makemyprior package with general explanations on how to use it in Section 4.
Section 5 gives more detailed examples showing how to use the package in
various situations. A summary and discussion is given in Section 6.

##  2 Motivating examples

We demonstrate the core ideas of total variance and hierarchical decomposition
of the total variance through two illustrative examples. These examples are
simplified versions of examples we use in Section 5.

One of the key quantities in quantitative genetics concerns the distribution
of observed variation to genetic and environmental sources. In this setting,
good intuition exists on the ratio of genetic to phenotypic variation, also
known as the heritability, whereas it is more difficult to define suitable
priors separately on the two variance parameters.

######  Example 2.1 (Genomic models).

Consider a group of $n$ individuals, where each individual $i$ has an observed
phenotype $y_{i}$. A simple genomic model is:

$y_{i}=\mu+a_{i}+\varepsilon_{i},\quad i=1,\dots,n,$ (1)

where $\mu$ is an intercept,
$\bm{a}=(a_{1},\dots,a_{n})^{\top}\sim\mathcal{N}_{n}(\bm{0},\sigma_{\mathrm{a}}^{2}\mathbf{A})$
is an additive genetic effect, and
$\bm{\varepsilon}=(\varepsilon_{1},\dots,\varepsilon_{n})^{\top}\sim\mathcal{N}_{n}(0,\sigma_{\mathrm{\varepsilon}}^{2}\mathbf{I}_{n})$
is environmental noise. The covariance matrix $\mathbf{A}$ is calculated based
on genetic sequencing of the $n$ individuals and is scaled so that
$\sigma_{\mathrm{a}}^{2}$ is representative of the variance arising from the
genetic effect; see Selle et al.,  (2019); Hem et al.,  (2021) for details.

In this simple model, two key quantities are the total variance
$\sigma_{\mathrm{P}}^{2}=\sigma_{\mathrm{a}}^{2}+\sigma_{\varepsilon}^{2}$,
also known as the _phenotypic variance_ , and the heritability
$h^{2}=\sigma_{\mathrm{a}}^{2}/(\sigma_{\mathrm{a}}^{2}+\sigma_{\varepsilon}^{2})$,
which is the proportion of the phenotypic variance explained by the genetic
effect. When expert knowledge is available about these two quantities, this
information can be directly exploited through a joint prior assigned to
phenotypic variance and heritability. A simple and intuitive visualization of
this parameterization is given by the tree in Figure 1(a) where the phenotypic
variance
$\sigma_{\mathrm{P}}^{2}=\sigma_{\mathrm{a}}^{2}+\sigma_{\varepsilon}^{2}$ in
the top (root) node is distributed to the additive genetic variance
$\sigma_{\mathrm{a}}^{2}$ and the environmental variance
$\sigma_{\varepsilon}^{2}$ in the two leaf nodes. Note that the intercept is
treated independently with a vague Gaussian prior. See Hem et al.,  (2021) for
a detailed description.

The idea of expressing a parameterization that is given in terms of total
variance and proportions of variances through a tree, extends to more complex
models with more random effects. For example, when analysing data arising from
designed experiments.

$a+\varepsilon$$\varepsilon$$a$ (a) Genomic model.

$a+b+c+\varepsilon$$\varepsilon$$a+b+c$$a$$b$$c$ (b) Latin square model.

$a+b+c+\varepsilon$$\varepsilon$$a+b+c$$a+b$$a$$b$$c$ (c) Latin square model.

Figure 1: Prior trees for (1(a)) the genomic model in Example 2.1, and (1(b),
1(c)) the latin square model in Example 2.2.

######  Example 2.2 (Latin square design).

Based on the ideas in Fuglstad et al.,  (2020), we assume that an agricultural
field is split into rows and columns resulting in a $9\times 9$ grid, and that
one out of $9$ different strengths of fertilizer is applied in each grid cell.
Outcomes $y_{i,j}$ are observed in row $i$ and column $j$ under treatment
$k[i,j]$, and modelled through

$y_{i,j}=\alpha+\beta\cdot
k[i,j]+a_{i}+b_{j}+c_{k[i,j]}+\varepsilon_{i,j},\quad i,j=1,\dots,9,$ (2)

where $\alpha$ is an intercept and $\beta\cdot k[i,j]$ is a linear effect of
treatment. $\alpha$ and $\beta$ are assigned vague Gaussian priors
$\mathcal{N}(0,\sigma=1000)$.
$\bm{a}=(a_{1},\dots,a_{9})^{\top}\sim\mathcal{N}_{9}(\bm{0},\sigma_{\mathrm{a}}^{2}\mathbf{I}_{9})$
is a row effect,
$\bm{b}=(b_{1},\dots,b_{9})^{\top}\sim\mathcal{N}_{9}(\bm{0},\sigma_{\mathrm{b}}^{2}\mathbf{I}_{9})$
is a column effect,
$\bm{c}=(c_{1},\dots,c_{9})^{\top}\sim\mathcal{N}_{9}(\bm{0},\sigma_{\mathrm{c}}^{2}\mathbf{I}_{9})$
is a treatment effect, and the residual noise is
$\bm{\varepsilon}=(\varepsilon_{1,1},\varepsilon_{1,2}\dots,\varepsilon_{9,9})^{\top}\sim\mathcal{N}_{81}(\bm{0},\sigma_{\varepsilon}^{2}\mathbf{I}_{81})$.
$\bm{a}$, $\bm{b}$ and $\bm{c}$ have sum-to-zero constrains. The individual
variances, $\sigma_{\mathrm{a}}^{2}$, $\sigma_{\mathrm{b}}^{2}$,
$\sigma_{\mathrm{c}}^{2}$ and $\sigma_{\varepsilon}^{2}$, are nuisance
parameters, and it may be difficult to have prior knowledge about them.

Figures 1(b) and 1(c) visualize two ways in which the total variance
$\sigma_{\mathrm{a+b+c}+\varepsilon}^{2}=\sigma_{\mathrm{a}}^{2}+\sigma_{\mathrm{b}}^{2}+\sigma_{\mathrm{c}}^{2}+\sigma_{\varepsilon}^{2}$
can be distributed to the individual variances $\sigma_{\mathrm{a}}^{2}$,
$\sigma_{\mathrm{b}}^{2}$, $\sigma_{\mathrm{c}}^{2}$ and
$\sigma_{\varepsilon}^{2}$ in the leaf nodes. Using Figure 1(b), we could
envision that, in the top split, shrinkage is applied to the latent variance
$\sigma_{\mathrm{a+b+c}}^{2}=\sigma_{\mathrm{a}}^{2}+\sigma_{\mathrm{b}}^{2}+\sigma_{\mathrm{c}}^{2}$
relative to the residual variance $\sigma_{\varepsilon}^{2}$ with the goal of
reducing the risk of overfitting. Then in the second split, we could envision
that we want to express ignorance about how the latent variance is distributed
to the individual variances $\sigma_{\mathrm{a}}^{2}$,
$\sigma_{\mathrm{b}}^{2}$ and $\sigma_{\mathrm{c}}^{2}$. Alternatively, using
Figure 1(c), we may want to express ignorance about how
$\sigma_{\mathrm{a+b}}^{2}=\sigma_{\mathrm{a}}^{2}+\sigma_{\mathrm{b}}^{2}$ is
distributed to $\sigma_{\mathrm{a}}^{2}$ and $\sigma_{\mathrm{b}}^{2}$, but
apply shrinkage to the variance of the treatment variance
$\sigma_{\mathrm{c}}^{2}$ relative to $\sigma_{\mathrm{a+b}}^{2}$. The full
details can be found in Fuglstad et al.,  (2020).

Examples 2.1 and 2.2 make it clear that in some cases it is natural to think
in terms of proportions of variances. However, explicitly writing out how the
proportions are defined may obfuscate the key ideas that one want to express.
Therefore, trees such as shown in Figure 1 are critical to define such priors.

We want to emphasize that a prior can be chosen in many ways, and there are no
wrong choices. Our message is that it is a choice to use the default prior, it
is a choice to use literature-based priors, and it is a choice to use a prior
based on prior and expert knowledge. We believe it is important to communicate
this and to make prior selection an active part of the Bayesian workflow.

##  3 Background

In this section, we present the necessary theoretical and methodological
background behind the makemyprior package and introduce terminology used
throughout the paper.

###  3.1 Hierarchical variance decomposition along a prior tree

This section introduces all information needed to use the makemyprior package
through the guide provided in the graphical user interface. We include more
specific details about the different priors in Section 3.2.

####  3.1.1 General definition of a prior tree

A prior tree is a directed acyclic graph consisting of a top node, split nodes
and leaf nodes, connected by directed edges. Each random model component is
represented by a leaf node. The top node in the tree represents total
variance, i.e., the sum of the variance of the leaf nodes (random effects). We
follow Fuglstad et al.,  (2020) and consider variance that is not explained by
fixed effects, and therefore omit fixed effects from the hierarchical
decomposition (HD) prior and the prior tree. Two or more leaf nodes are
combined in a split node in a way that reflects the hierarchical model and our
prior beliefs about how the total variance is distributed. The split nodes
represent a variance proportion.

In some cases we only want to include a subset of the model components in the
same prior tree, and in that way have several trees. An example is the animal
model (e.g., Holand et al., , 2013), a mixed model that is usually used to
decompose environmental and genetic variances in animal populations. The
genetic contribution might be split into several sub-components, like
additive, dominance and mutational variance. This is an extension of the
genomic model in Example 2.1. We may have a good intuition on the absolute
magnitude of the variation that is explained by the environment and by all the
genetic contributions separately, but only have knowledge about the relative
magnitude of the genetic variation. The genetic contributions are often
confounded, and it is useful to split the variance of the genetic contribution
using a prior tree. In that way we have both variance parameters, total
variance parameters, and variance proportions.

If a variance component is assigned an individual prior, its corresponding
random effect is represented by a singleton: a node not connected to any other
nodes. The singletons can be considered as trees with only one node. Several
prior trees gives a prior forest. We refer to the forest of trees as the prior
tree structure of the model. Each tree in a tree structure is associated with
their own joint prior distribution, independent of the priors belonging to the
other trees. The tree structure is made based on prior knowledge about the
model, data and problem at hand.

####  3.1.2 Defining priors for the split nodes: Shrinkage versus ignorance

Given a tree structure describing the distribution of the variance in the
model, we can use the rest of our pre-existing knowledge to steer the variance
to the different model components by choosing suitable priors for the
different parameters belonging to the prior tree.

A good prior distribution can improve the robustness of the inference, by
helping to avoid estimating spurious effects. We apply the penalized
complexity (PC) priors of Simpson et al.,  (2017) as they shrink towards a so-
called base model and are thus robust by design. For a random effect with zero
mean and variance parameter $\sigma^{2}$, the base model is a model where
$\sigma^{2}=0$, meaning the effect does not contribute in the model. This
gives a simpler model with one effect less. When the parameter varies, we move
away from the base model, and this deviation is penalized, ensuring that we do
not overfit the model. The hyperparameters are chosen using prior knowledge on
an interpretable scale relevant for the parameter. For example, a tail
probability for a standard deviation.

Any model parameter can be assigned a PC prior. For a variance proportion
parameter, the base model can be any value in $[0,1]$. In the corner cases $0$
and $1$, we are in a similar situation as for the variance parameter, where
the base model represents a model with fewer random effects as one of them
shrinks away. The heritability $h^{2}$ in Example 2.1 is a variance
proportion. We can imagine that a geneticist says that the heritability is
around $0.4$ for this phenotype and species, but that it is unsure whether the
additive effect contributes or not. Then we can use a PC prior that shrinks
the heritability to $0$, which gives only residual effect in the base model,
and ensure that the prior has a median of $0.4$, i.e.,
$\text{Prob}(h^{2}>0.4)=0.5$.

If the base model is somewhere between $0$ and $1$, we express that both
effects involved in the variance proportion should be present in the model. If
the geneticist is certain that the additive effect is present in the model, we
can choose a PC prior that shrinks the heritability $h^{2}$ towards $0.4$, and
with median at the same value: $\text{Prob}(h^{2}>0.4)=0.5$. This corresponds
to $40$% additive and $60$% residual effect in the base model. See Section 3.2
and Fuglstad et al.,  (2020) for details.

By using a multi-split, a split with more than two child nodes, the user
expresses no strong opinions on how the variance is distributed among the
components involved in the split. Then we use a prior that assigns an equal
amount of variance to each of the model effects involved in the split through
the ignorant symmetric Dirichlet prior. This ignorant prior can also be used
on a dual split, in which case it reduces to a uniform prior on $[0,1]$, and
with that the user expresses that no prior information is used in that split.

####  3.1.3 Defining priors for top nodes and singletons

Appropriate prior distributions for variance parameters varies with the
likelihood. We describe the train of thought when specifying prior
distributions for variance parameters when the likelihood is Gaussian,
binomial and Poisson.

In a model with a Gaussian likelihood, the total variance is usually easy to
identify and does not need an informative prior. This is a parameter close to
the data in the model and settles with a vague prior (Goel and Degroot, ,
1981; Gelman et al., , 2020). For the total variance in the top nodes,
Fuglstad et al.,  (2020) recommend the scale-invariant, improper Jeffreys’
prior when all model components are involved in one single tree. This prior
does not require any hyperparameters and is straight-forward to use.

In cases where the user has specific knowledge about the total variance, a
proper prior can be used to include this knowledge in the prior. In cases
where the nodes are not all involved in the same tree and we have a prior
forest, a scale-invariant prior is not meaningful, and an improper prior may
lead to an improper posterior. This implies that singletons always must be
assigned a proper prior. Due to its desirable shrinking properties Fuglstad et
al.,  (2020) recommend the PC prior (Simpson et al., , 2017) for the
variances, where the hyperparameters can be selected using a tail probability
of the standard deviation: $\mathrm{Prob}(\sigma>U)=\alpha$. We denote this
prior $\mathrm{PC}_{\mathrm{0}}(U,\alpha)$, where the subscript indicates that
the shrinkage is towards $0$. For the genomic model in Example 2.1, assume we
have prior or expert knowledge saying it is unlikely that the total variance
in the observed data is greater than $4$. We want to use this knowledge, and
choose $U=\sqrt{4}$. The value we choose for $\alpha$ says something about how
certain we are in the value of $U$. Using $\alpha=0.05$ is a suitable choice,
so that $\mathrm{Prob}(\sigma>\sqrt{4})=0.05$ and
$\sigma\sim\mathrm{PC}_{\mathrm{0}}(\sqrt{4},0.05)$.

Other likelihoods require proper priors on all variances, also the total
variance, as scale-invariance is not meaningful for data that are not
Gaussian. Again we follow the recommendation of Fuglstad et al.,  (2020) and
suggest PC priors. However, instead of choosing a prior using the upper tail
probability of the standard deviation, we think on a different scale than for
Gaussian data, and choose a credible interval for the variance parameter on a
suitable scale. Then we transform the interval into an upper tail probability
$\mathrm{Prob}(\sigma>U)=\alpha$. For both binomial likelihood with logit link
and Poisson likelihood with log link, an exponential scale is appropriate.
This will correspond to thinking on odds-ratio scale for binomial data, and on
the scale of the data for Poisson data.

Assume we have a model with linear predictor $\eta_{i}=\mu+a_{i}+b_{i}$ and
binomial likelihood with logit link function. An intuitive way of choosing a
prior for the total variance of $a_{i}+b_{i}$ is to choose an equal-tailed
credible interval for the effect of the random effects on the odds-ratio,
$\exp(a_{i}+b_{i})$, i.e., $\mathrm{Prob}(l<\exp(a_{i}+b_{i})<u)=p$ (Fong et
al., , 2010). For example, we can say we want
$\mathrm{Prob}(0.1<\exp(a_{i}+b_{i})<10)=0.9$. This corresponds to a 90%
credible interval $[0.1,10]$ for $\exp(a_{i}+b_{i})$. The idea is the same for
a Poisson likelihood with a log link: We can think on the effect of the random
effects on the relative risk.

###  3.2 Priors for variance proportions

Details around the prior distributions discussed in the previous section are
presented here. While the makemyprior package can be used without knowing
these details, they are needed to fully understand the hierarchical
decomposition (HD) prior framework.

####  3.2.1 Penalized complexity

The penalized complexity (PC) priors of Simpson et al.,  (2017) are based on
principles independent of model and application. They use the distance
$d(\cdot)$ between the base model and a flexible extension of the base model,
measured using the Kullback-Leibler divergence. In the base model the
parameter of interest $\theta$ is fixed to $\theta_{0}$. In the flexible model
$\theta$ is allowed to vary, and the PC prior induces shrinkage towards the
base model, which gives a robust prior that aids to avoid overfitting. As in
Simpson et al.,  (2017), we use an exponential prior for the distance
$d(\cdot)$ between the two models, and transform this to a prior for the
desired parameter. This means that the PC prior always is an exponential
distribution on the distance, while for the parameter in question $\theta$ the
distribution varies with parameterization, choice of base model and covariance
matrix structures. The PC prior does in general not have an analytical
expression. The parameter $\theta$ can be a variance or standard deviation
(Simpson et al., , 2017), a variance proportion (Fuglstad et al., , 2020; Hem
et al., , 2021), or, for example, a correlation parameter (Guo et al., ,
2017). In makemyprior we consider standard deviations $\sigma$ (and variances
$\sigma^{2}$) together with variance proportions $\omega$. For a standard
deviation the distance is simply $d(\sigma)=\sigma$ (Simpson et al., , 2017).
For a variance proportion parameter the distance will be a function of the
covariance matrices of the random effects involved in the split, see Fuglstad
et al.,  (2020) for details.

####  3.2.2 Shrinkage priors for variance proportions

We construct the joint hierarchical decomposition (HD) prior using a bottom-up
approach following the prior tree, and the prior will thus be dependent on
prior tree structure. The distance measure $d(\cdot)$ for a variance
proportion depends on the covariance matrices of the effects of the child
nodes in a split (Fuglstad et al., , 2020), and the covariance matrix of a
split node will be a function of the variance proportions(s) involved in this
split. This means we must condition on the variance proportions associated
with splits lower in the tree (if any), and that each prior depends on choices
and covariance matrices at that and lower levels. We omit the dependence of
tree structure, covariance matrices and prior choices for other splits in the
notation of the PC prior for readability. 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{\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\sim\mathrm{PC}_{\mathrm{0}}(m)$ with $\mathrm{Prob}(\omega>m)=0.5$ so
that $m$ defines the median, and shrinkage towards $\omega=0$, i.e., the base
model is a model with only $\bm{\varepsilon}$.

  * •

$\omega\sim\mathrm{PC}_{\mathrm{1}}(m)$ with $\mathrm{Prob}(\omega>m)=0.5$ so
that $m$ defines the median, and shrinkage towards $\omega=1$, i.e., the base
model is a model with only $\bm{a}$.

  * •

$\omega\sim\mathrm{PC}_{\mathrm{M}}(m,c)$ with $\mathrm{Prob}(\omega>m)=0.5$
and
$\mathrm{Prob}(\mathrm{logit}(1/4)<\mathrm{logit}(\omega)-\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=m$,
i.e., the base model is a combination of the effects $\bm{a}$ and
$\bm{\varepsilon}$.

Note that $\mathrm{PC}_{\mathrm{1}}(m)$ for $\omega$ is equivalent to
$\mathrm{PC}_{\mathrm{0}}(1-m)$ for
$1-\omega=\frac{\sigma_{\varepsilon}^{2}}{\sigma_{\mathrm{a}}^{2}+\sigma_{\varepsilon}^{2}}$.
Since the PC prior is a prior put on the distance between two models, and then
transformed to the parameter of interest, we do not distinguish the notation
between the PC prior on a standard deviation and variance parameter, as it
will result in the same prior.

In Figure 2 we show examples for the four different priors. Code producing
these graphs can be found in Appendix A. The shape of the PC prior on a
standard deviation, shown in Figure 2(a), is independent of the
hyperparameters and other hyperparameters will simply give a rescaling of the
axes.

Figure 2(b) shows a prior where we have shrinkage towards $\omega=0$, which
corresponds to a model with only residual variance. In Figure 2(c) we show a
prior where with shrinkage towards $\omega=1$, giving a model with only group
effect. This base model has a singular covariance matrix, giving a distance
measure that is infinite. In such cases we cannot have a median that is
further than $0.25$ from the base model (see Fuglstad et al.,  (2020, Theorem
1) for details). In this case that means that the median cannot be smaller
than $0.75$.

The covariance matrix of the base model is always non-singular for a
$\mathrm{PC}_{\mathrm{M}}$ prior (shown in Figure 2(d)). The concentration
parameter $c$ in $\mathrm{PC}_{\mathrm{M}}(m,c)$ measures how certain we are
about the prior median. A concentration of less than $0.5$ will indicate that
a prior with median at $\omega=0.5$ is less concentrated around $0.5$ than a
uniform distribution on $[0,1]$ would be, and we set the lower limit of $c$ to
$0.5$. As this parameter says something about how much of the distribution
mass is in an interval that is smaller than the parameter space
($0\leq\omega\leq 1$), the distribution will not change much when $c$
approaches $1$.

(a) $\sigma\sim\mathrm{PC}_{\mathrm{0}}(3,0.05)$.

(b) $\omega\sim\mathrm{PC}_{\mathrm{0}}(0.25)$.

(c) $\omega\sim\mathrm{PC}_{\mathrm{1}}(0.75)$. Note that the base model has a
singular covariance matrix.

(d) $\omega\sim\mathrm{PC}_{\mathrm{M}}(0.25,0.85)$.

Figure 2: Examples of the different PC priors for a random intercept model
$y_{i,j}=a_{i}+\varepsilon_{i,j}$ for $i,j=1,\dots,10$.
$\omega=\sigma_{\mathrm{a}}^{2}/(\sigma_{\mathrm{a}}^{2}+\sigma_{\varepsilon}^{2})$.

The difference between the
$\mathrm{PC}_{\mathrm{0}}$/$\mathrm{PC}_{\mathrm{1}}$ and
$\mathrm{PC}_{\mathrm{M}}$ priors is how certain the user is in the prior
knowledge. In the genomic model in Example 2.1, a geneticist told us that the
heritability $h^{2}$ is around $0.4$. If the contribution of the additive
effect is unclear, for example due to a small data sample, we can use
$\mathrm{PC}_{\mathrm{0}}(0.4)$ for $h^{2}$, which has shrinkage towards $0$
(only residual effect) and median $0.4$. If we are certain that the additive
genetic effect is contributing to the variation in the observed phenotype, we
use a $\mathrm{PC}_{\mathrm{M}}$ prior with median $m=0.4$, and we choose the
concentration parameter $c$ based on how strongly we believe in this value.
This gives a prior with shrinkage towards a heritability of $0.4$.

####  3.2.3 Ignorance priors for variance proportions

If we want to express ignorance about the variance distribution in a split, we
use a symmetric Dirichlet prior distribution. For a split with $p\geq 2$
children, this prior is given by:

$\pi(\bm{\omega};p)=\mathrm{Dirichlet}(p)=\frac{\Gamma\left(p\alpha\right)}{\Gamma\left(\alpha\right)^{p}}\left(\prod_{i=1}^{p}\omega_{i}\right)^{\alpha-1},$

where $\bm{\omega}$ is the vector of variance proportions involved in the
split with $0<\omega_{i}<1$ for $i=1,\ldots,p$ and
$\sum_{i=1}^{p}\omega_{i}=1$. Further, $\Gamma(\cdot)$ denotes the gamma
function and $\alpha>0$ is chosen so
$\mathrm{P}(\mathrm{logit}(1/4)<\mathrm{logit}(w_{i})-\mathrm{logit}(1/p)<\mathrm{logit}(3/4))=1/2$
for $i=1,\dots,p$ (if this is achieved for one $i$, it is by symmetry achieved
for all $i$). This prior assigns equal amount of variance to each model
component. For a split consisting of $p$ nodes we denote this as
$\bm{\omega}\sim\text{Dirichlet}(p)$ for each proportion in the split. A
multi-split is by itself a way of showing ignorance, and all multi-splits are
given this Dirichlet prior.

####  3.2.4 Shrinkage priors for multi-splits

There may be situations where the user wants to assign unequal amounts of
variance to the components in a multi-split. In that case, the user has
opinions about the variance decomposition in the split, and an ignorant multi-
split with a Dirichlet prior is not suitable. Instead we transform the multi-
split to several dual splits and assign a PC prior to each of the dual splits.
In Example 2.2, assume that we want a (20, 30, 50) division of the total
variance of $\bm{\alpha}$, $\bm{\beta}$ and $\bm{\gamma}$. We achieve this by
first splitting the variance 50/50 between $\bm{\gamma}$ and $\bm{\alpha}$,
$\bm{\beta}$ with a $\mathrm{PC}_{\mathrm{M}}(0.5,c)$ prior, and then dividing
the variance of $\bm{\alpha}$ and $\bm{\beta}$ 40/60 with a
$\mathrm{PC}_{\mathrm{M}}(0.4,c)$ prior for some suitable value of $c$. The
corresponding tree structure is shown in Figure 1(c). We could also have
chosen another effect to split off first, but Fuglstad et al.,  (2020) show
that this order does not have much impact on the resulting prior.

####  3.2.5 Pre-computing the marginal priors

To ease the computation of the joint PC prior, we follow Fuglstad et al.,
(2020) and condition on the base models of the splits below in the tree,
instead of on the parameters themselves. In this way we can pre-compute the
marginal priors for each split in the tree. The base model for the Dirichlet
distribution is equal variance to each component, i.e., for a split with $p$
components the base model is $1/p$ of the variance in the split node to each
child node. A split with a Dirichlet prior does not use information from lower
levels in the tree.

##  4 Software

Throughout this section, the use of the R package makemyprior is exemplified
by the following model.

######  Model 1 (Example model).

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

$\displaystyle y_{i,j}|\eta_{i,j},\sigma_{\varepsilon}^{2}$
$\displaystyle\sim\mathcal{N}(\eta_{i,j},\sigma_{\varepsilon}^{2}),$
$\displaystyle\eta_{i,j}$ $\displaystyle=\mu+x_{i}\beta+a_{i}+b_{j},$

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.
$\varepsilon_{1},\varepsilon_{2},\dots,\varepsilon_{n}\overset{\text{iid}}{\sim}\mathcal{N}(0,\sigma_{\varepsilon}^{2})$
are residuals.

The data for this model could be some measures made for year $i$ and age group
$j$ with a time-specific covariate $x_{i}$. We do not interpret the meaning of
the data in this section.

We include four different tree structures for this model in Table 1. Several
likelihoods, latent models and prior distributions are available in
makemyprior, which can be listed with the function:

makemyprior_models(type = c("prior", "latent", "likelihood"), select = NULL)

###  4.1 Specifying the linear predictor and preparing the data object

First, the linear predictor is specified using a formula object with a syntax
similar to e.g. lm() from stats and inla() from INLA. Covariates are included
directly by name in the formula, and mc() is used to include information about
each random effect:

mc(label, model = "iid", ...)The main arguments are label (name of the effect)
and model (the type of latent model effect, i.i.d. is the default). The other
arguments depend on the choice of latent model, see documentation for details.

For Model 1 where both $\bm{a}$ and $\bm{b}$ are i.i.d., the formula is:

R> formula <\- y ~ x + mc(a) + mc(b) The intercept $\mu$ is included by
default, but can be removed with -1. The residual effect for a Gaussian
likelihood should not be specified in the formula.

The second step is to gather data and create a data object as either a
data.frame or list with names corresponding to the elements of the formula.
For Model 1, we need the observed response y, the covariate x, and indexes for
a and b (these must be specified as integers). A simple simulated dataset is:

R> p <\- 10R> m <\- 10R> n <\- m*pR>R> set.seed(1)R> data <\- list(a =
rep(1:p, each = m),\+ b = rep(1:m, times = p),\+ x = runif(n))R> data$y <\-
data$x + rnorm(p, 0, 0.5)[data$a] +\+ rnorm(m, 0, 0.3)[data$b] + rnorm(n, 0,
1)We recommend the use of short names for the input data because these names
need to be used to refer to model components in later steps. An example is
"rain" instead of "rainfall_august_2020". Note that the observations $y_{i,j}$
are not used to make the prior.

###  4.2 Exploring and selecting the prior graphically

| Tree structure | Text string
---|---|---
Prior 1 |  |  $a$$b$$\varepsilon$
---
| "(a); (b); (eps)"
---
Prior 2 |  |  $\varepsilon$$a+b$$a$$b$
---
| "s1 = (a, b); (eps)"
---
Prior 3 |  |  $a+b+\varepsilon$$\varepsilon$$a$$b$
---
| "s1 = (a, b, eps)"
---
Prior 4 |  |  $a+b+\varepsilon$$\varepsilon$$a+b$$a$$b$
---
| "s1 = (a, b); s2 = (s1, eps)"
---
Table 1:  Four tree structures for Model 1 with text strings specifying them.
The names s1 and s2 on the splits (variance proportions) are chosen by the
user in the initial specification of the prior, and are used in a nested
formulation to specify the prior tree structure.

We provide a graphical user interface (GUI) where the user can construct the
prior in an interactive way. The GUI is implemented as a shiny app (Chang et
al., , 2020) running locally. It allows the user to first define the tree
structure, and then be guided sequentially through the steps of selecting
priors for each split, singleton and top node. The shiny apps are useful to,
for example, display and investigate results from analyses. For example, the
user can customize graphs and tables in a simple way. Such apps are used by
Depaoli et al.,  (2020) to show why prior sensitivity analysis is important,
and by Smid and Winter,  (2020) to let users explore the impact of prior
distributions in inference. However, to the extent of our knowledge, as of
today there are no packages or apps that allows the user to use a shiny app
(or similar) to specify priors for custom models and data and directly carry
out inference.

The first step is to initialize a prior object. This is done with the function
make_prior():

make_prior(formula, data, family = "gaussian", prior = list(), intercept_prior
= c(), covariate_prior = list())formula and data are the objects created in
Section 4.1, and family is the likelihood ("gaussian" is the default;
"binomial" and "poisson" are also available). The prior argument is not
relevant when the GUI is used, and its description is deferred to Section 4.3.
intercept_prior and covariate_prior specifies the mean and standard deviation
for the Gaussian priors on the intercept and covariate coefficients. The
default in makemyprior is $\mathcal{N}(0,\sigma=1000)$ for both, and the
coefficient priors are specified as a named list with names corresponding to
the covariate names.

For Model 1, we can create a prior object with the following command:

R> prior <\- make_prior(formula, data, family = "gaussian",\+ intercept_prior
= c(0, 1000),\+ covariate_prior = list(x = c(0, 100)))

Warning message:Did not find a tree, using default tree structure instead.
This gives a prior with a single tree with one split as shown in Prior 3 in
Table 1, which is the default setting. A warning makes the user aware that the
default prior is chosen. For Gaussian likelihoods, Jeffreys’ prior is by
default set for the total variance. Let
$\sigma_{\mathrm{a+b}+\varepsilon}^{2}$ denote total variance, and
$\bm{\omega}=\big{(}\omega_{\frac{\mathrm{a}}{\mathrm{a+b}+\varepsilon}},\omega_{\frac{\mathrm{b}}{\mathrm{a+b}+\varepsilon}},1-\omega_{\frac{\mathrm{a}}{\mathrm{a+b}+\varepsilon}}-\omega_{\frac{\mathrm{b}}{\mathrm{a+b}+\varepsilon}}\big{)}$
describe the attribution of variance to the three different sources, then the
initial choice of priors is:

$\bm{\omega}\sim\mathrm{Dirichlet}(3)\quad\text{and}\quad\sigma_{\mathrm{a+b}+\varepsilon}^{2}\sim\mathrm{Jeffreys^{\prime}}.$
(3)

The intercept has a $\mathcal{N}(0,\sigma=1000)$ prior, and the covariate
$\beta$ a $\mathcal{N}(0,\sigma=100)$ prior.

Figure 3: Screenshot of the GUI in makemyprior with the default prior
(Equation 3).

The function makemyprior_gui() allows the user to select the desired prior
tree structure and choose prior distributions interactively:

makemyprior_gui(prior, guide = FALSE, no_pc = FALSE)This function takes the
arguments prior, created with make_prior() earlier, guide which specifies
whether or not the guide should automatically start (the guide can be started
at any time), and no_pc. The PC prior is computed using the covariance matrix
structure of the model components, which may be slow for large models. For a
better user experience, the user can turn off the computation of the PC prior
in the GUI using no_pc = TRUE. The prior will be computed upon closing, and
this will only affect the plotting of the prior in the GUI.

For Model 1, we start the GUI by running:

R> new_prior <\- makemyprior_gui(prior) This saves the new prior chosen to the
variable new_prior when we close the GUI. Figure 3 shows a screenshot of the
GUI for Model 1 for $m=p=10$ and the default prior. By selecting nodes and
using the buttons located in the panels to the left, the user can create the
desired tree structure and choose the desired priors for this tree.

(a) Initial window.

(b) Question for prior elicitation.

(c) Question for prior elicitation.

Figure 4: Screenshots from the guide in the makemyprior GUI.

The user can click ”Begin guide” to the lower left to be guided through the
prior construction. Figure 4(a) shows the initial window of the guide, which
consists of two main steps: First, the user is helped to create the desired
tree structure. Note that every time the tree structure is modified, both when
using the guide and otherwise, all splits are set to have the default
Dirichlet prior. After deciding on a tree structure, the user will be taken
through the prior tree structure in a step-wise manner. For each split, top
node and singleton they will be asked simple questions about the existing
prior knowledge about the corresponding parameter. This eases the process of
using the knowledge to create the prior, as the user does not need to think in
terms of prior distributions directly, but only needs to answer questions such
as shown in Figure 4(b) and 4(c). After modifying the prior to have the
following distributions:

$\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),$
(4)

the GUI looks like the screenshot in Figure 5. The chosen prior distributions
and the connection between the parameterization and model variances are easily
seen to the left.

Figure 5: Screenshot of the GUI in makemyprior when the prior in Equation 4 is
chosen.

The GUI is intuitive and contains a thorough description of the options, and
we do not explain the features in detail here. We instead recommend using the
guide in the GUI to get familiar with the package. A summary of the prior
object can be printed with:

R> summary(new_prior)

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] ~ PC1(0.75)Total variance priors:sqrt(V)[eps_a_b]
~ PC0(3, 0.05)Covariate priors: intercept ~ N(0, 1000^2), x ~ N(0, 100^2)y ~ x
+ mc(a) + mc(b)

###  4.3 Selecting the prior non-graphically

If the prior has been constructed with the GUI, this section can be skipped.
Here we describe an alternative way to specify the prior without the use of
the GUI. This is done with same function we use to make the default prior,
make_prior(), by specifying the prior argument. prior is a named list with the
following arguments:

tree

The tree structure as a string. A split is specified as s1 = (a, b), where s1
represents a split node and can be any name except names of the input data in
data and the reserved eps, which is used for residuals for a Gaussian
likelihood. Short names are recommended. Note that these split names are just
used in the initial specification. The child nodes for each split are included
in parentheses separated by commas, and each split is separated by semicolons.
Singletons are included as (a). Examples of strings for different tree
structures for Model 1 are shown in Table 1.

V

A named list with information on the priors on each top node and singleton,
i.e., all variances. Options are "pc", "jeffreys", "invgam", and "hc" (Half-
Cauchy). The names in the list are the top node and singleton names from the
tree argument.

w

A named list with information on the priors on each split, i.e., all variance
proportions. The names in the list are the split node names from the tree
argument. Options are "pc0", "pc1", "pcM" and "dirichlet".

V and w must have the following structure for each element in the list:
list(prior = prior_name, param = parameter_vector), except for Jeffreys’ and
Dirichlet priors, where param is not specified (as the distributions do not
have hyperparameters). See Section 5 for how to specify priors in specific
examples, and makemyprior_models() for more details about the different prior
distributions.

The prior in Equation 4 for Model 1 can be specified as:

R> 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 = "pc0", param = 0.25)),\+ V = list(s2 = list(prior =
"pc0", param = c(3, 0.05)))\+ ),\+ covariate_prior = list(x = c(0, 100)))R>
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] ~ PC1(0.75)Total
variance priors:sqrt(V)[eps_a_b] ~ PC0(3, 0.05)Covariate priors: intercept ~
N(0, 1000^2), x ~ N(0, 100^2) The names s1 and s2 are chosen by the user in
the initial specification of the prior, and can be any names that are not used
for the data or eps (reserved for residuals). Note that s1 and s2 are only
used as a link between the splits and priors, and have been automatically
changed to a_b and eps_a_b by make_prior(). The order we list the children for
each split node in the tree argument decides which way we have shrinkage with
the PC priors. For s1 = (a, b), $\mathrm{PC}_{\mathrm{0}}(m)$ shrinks effect
$\bm{a}$ ($\omega_{\frac{\mathrm{a}}{\mathrm{a+b}}}=0$ as base model),
$\mathrm{PC}_{\mathrm{1}}(m)$ shrinks effect $\bm{b}$ and
$\mathrm{PC}_{\mathrm{M}}(m,c)$ gives shrinkage towards $m\bm{a}+(1-m)\bm{b}$.
All three has median at $\omega_{\frac{\mathrm{a}}{\mathrm{a+b}}}=m$. Note
that $\mathrm{PC}_{\mathrm{0}}(m)$ on
$\omega_{\frac{\mathrm{a}}{\mathrm{a+b}}}$ is equivalent to
$\mathrm{PC}_{\mathrm{1}}(1-m)$ on $\omega_{\frac{\mathrm{b}}{\mathrm{a+b}}}$.

All top nodes, singletons and split nodes without a specified prior will get
the default prior. The default settings in makemyprior are chosen based on the
findings of Fuglstad et al.,  (2020) to ensure robust inference:

  * •

If no prior is specified (neither tree structure nor priors), the prior will
be a joint prior where all latent components (including a possible residual
effect) get an equal amount of the total variance in the prior through the
symmetric Dirichlet prior, and the default total variance prior.

  * •

The default prior on the total variance (top node) varies with likelihood:

    * –

Jeffreys’ prior for Gaussian likelihood for a tree structure with one tree,
$\mathrm{PC}_{\mathrm{0}}(3,0.05)$ otherwise.

    * –

$\mathrm{PC}_{\mathrm{0}}(1.6,0.05)$ for binomial likelihood.

    * –

$\mathrm{PC}_{\mathrm{0}}(1.6,0.05)$ for Poisson likelihood.

  * •

The default prior on an individual variance (singleton) varies with
likelihood:

    * –

$\mathrm{PC}_{\mathrm{0}}(3,0.05)$ for Gaussian likelihood.

    * –

$\mathrm{PC}_{\mathrm{0}}(1.6,0.05)$ for binomial likelihood.

    * –

$\mathrm{PC}_{\mathrm{0}}(1.6,0.05)$ for Poisson likelihood.

  * •

The default prior on a variance proportion (split node) is a Dirichlet prior
assigning equal amount of variance to each of the model components involved in
the split.

The reasoning behind these choices are as follows. For the variance
proportions, there is no way of knowing what behaviour is desired in general.
Therefore we use the ignorant Dirichlet prior as the default. A standard
Gaussian distribution with mean $0$ and variance $1$ will have close to all
the density mass between $-3$ and $3$, and $\mathrm{PC}_{\mathrm{0}}(3,0.05)$
is a vague prior for the standard deviation in such a distribution. We use
Jeffreys’ prior when applicable (see Section 3.1). To choose the default
variance (standard deviation) prior for other likelihoods, we follow the idea
of Fong et al.,  (2010) and use a credible interval on a suitable scale. The
default prior for all variance parameters (both for top nodes and singletons)
in makemyprior for binomial and Poisson likelihoods is a PC prior with a 95%
credible interval between $0.2$ and $5$ for the multiplicative effect on the
odds ratio and risk, respectively. This is obtained with a
$\mathrm{PC}_{\mathrm{0}}(1.6,0.05)$ prior. We want to emphasize that before
selecting the default prior, both when using makemyprior and otherwise, you
should stop and think about whether or not it is suitable for your model and
data.

###  4.4 Performing inference

We include functions for inference that are compatible with the prior object
obtained from make_prior() (and makemyprior_gui()). Both Stan (Carpenter et
al., , 2017) through rstan (Stan Development Team, , 2020) and Integrated
Nested Laplace Approximations (INLA, Rue et al., , 2009) through INLA (see
www.r-inla.org) can be used for the inference.

Stan is a probabilistic programming language, where Hamiltonian Monte Carlo
(HMC) is used to sample from the posterior distribution (Carpenter et al., ,
2017). Stan implements HMC using the No U-Turn Sampler (NUTS, Hoffman and
Gelman, , 2014). NUTS reduces the need for tuning of the sampler, making it
easy to use as no manual settings are needed for the algorithm to run, and the
user only needs to provide the joint prior and likelihood model, implemented
in a programming language similar to C++. We provide Stan code for fitting
latent Gaussian models with certain likelihoods and latent effects. The
internal parameterization in the provided Stan code is log-variance, and the
prior is transformed from the parameterization given by the prior tree
structure to log-variances.

INLA is a non-sampling based method for doing fast and efficient Bayesian
inference on latent Gaussian models (Rue et al., , 2009), utilizing Gaussian
Markov random fields (GMRFs) with sparse precision matrices, which gives
computational benefits through the Markov property. The INLA method
approximates the posteriors by a mixture of Gaussian distributions and applies
a skewness correction to the marginals (Rue et al., , 2017). It is easy and
straight-forward to use for inference, and can fit models with a broad range
of latent effects. The internal parameterization of the model parameters in
INLA is log-precision, and in the same way as for the provided Stan code the
parameterization following the prior tree is transformed to fit INLA.

Some common latent models are included in the code for the package: i.i.d.
("iid"), random walk of first ("rw1") and second ("rw2") order, Besag
("besag"), and effects with structured covariance matrices ("generic0"). The
likelihood family was specified in make_prior(). In Section 5 we show how the
inference can be performed. Here we describe the functions that can be used
for inference.

####  4.4.1 Inference with Stan

Stan is a flexible tool for inference, however, it requires the user to write
custom code for the model that is to be fitted. makemyprior contains pre-
written Stan-code that can be used to do inference on latent Gaussian models
with a selection of latent models. We recommend to compile the Stan-code
before doing inference with Stan. This can be done with the following
function:

compile_stan(save = FALSE, path = NULL)where save indicates whether or not to
save the compiled object (must be set this to TRUE to avoid recompiling the
code every time inference is performed). path is only necessary if the default
location for saving the compiled object is not possible to use (see
?compile_stan for details). For inference with Stan we use the following
function:

inference_stan(prior_obj, use_likelihood = TRUE, print_prior = TRUE, ...)The
first argument is the prior object from make_prior() or makemyprior_gui(). The
user can specify whether to include the likelihood (use_likelihood = TRUE) or
not (FALSE). In the latter case we sample from the prior distribution.
print_prior (TRUE by default) prints details about the chosen prior.
Additional arguments that is sent directly to the rstan function sampling()
can be specified for the inference. Useful arguments include iter (total
number of iterations for each chain), warmup (number of iterations for the
warm-up), chains (number of chains), seed (for reproducibility), and control
(for specifying algorithm tuning parameters).

The internal parameterization in the Stan-code included in the package is log-
variance, however, since Stan works with samples we can look at any
parameterization we want by transforming the log-variances. For using other
latent models or more complex models than the ones provided in the included
Stan-code (see above), the user must write customized Stan-code, see Section
4.6.

####  4.4.2 Inference with INLA

For inference with INLA we use the following function:

inference_inla(prior_obj, use_likelihood = TRUE, print_prior = TRUE, ...)The
first three arguments are the same as in inference_stan(). Additional
arguments can be fed to the INLA function inla(). Useful arguments include
Ntrials for the binomial likelihood, used to specify the amount of trials,
where the response is the number of successes.

###  4.5 Visualizing priors and posteriors

We offer several functions to visualize the prior and posterior distributions.
The prior distributions for the random effects on the tree structure
parameterization can be plotted with plot_prior(obj) which take an object from
make_prior(), makemyprior_gui(), inference_stan() or inference_inla() as
input. The posterior distributions can be displayed with

plot_posterior_variance(obj)plot_posterior_stdev(obj)plot_posterior_precision(obj)obj
is an object from inference_stan() or inference_inla().

The posterior distributions of random effects from inference with Stan can be
plotted with:

plot_posterior_stan( obj, param = c("prior", "variance", "stdev",
"precision"), prior = FALSE)Here, obj is an object from inference_stan(),
param specifies which parameterization the plots should have where param =
"prior" gives the posterior on the same parameterization as the prior. prior
indicates whether or not to plot the prior together with the posterior for
param = "prior". The total variance prior will only be plotted if it is not
Jeffreys’ prior. Fixed effect posteriors can be plotted with
plot_posterior_fixed(obj). More details about visualization can be found with
?makemyprior_plotting.

###  4.6 More complex models in Stan

Latent models may have parameters that are not variances, such as
correlations. These non-variance parameters are handled independently in the
HD prior (Fuglstad et al., , 2020). The Stan code included in makemyprior is
only applicable for certain latent models and likelihoods (see Section 4.4).
We provide a “skeleton” code and a description on how the user can write
custom Stan-code and include the joint prior created with make_prior(). This
can be accessed with:

create_stan_file(location = "")location is a string to a path where a folder
with necessary files will be stored. The user can edit the code and include
custom latent components etc. We do not include details on this, as it will be
highly model specific and is merely an offer to the users who want to apply
the HD prior in more advanced models.

##  5 Using makemyprior: Examples

In this section we provide three examples where we use the makemyprior package
to construct priors and run inference. Two examples are with Gaussian
responses, and one is with Binomial responses. We have used Stan for the
inference (with inference_stan()), but the procedure is the same for inference
with INLA (using inference_inla() instead).

###  5.1 Gaussian responses

####  5.1.1 Genomic selection in wheat breeding

This is an extended version of the model in Example 2.1. In addition to the
additive genetic effect $\bm{a}$, we now also include two nonadditive effects:
Dominance $\bm{d}$ and additive-by-additive epistasis $\bm{x}$. This example
is taken from Hem et al.,  (2021). The response $y_{i}$ is grain yield for
individual $i$. We utilize the expert knowledge elicited from experts in the
field, and create a prior distribution reflecting this knowledge. We model the
response as:

$y_{i}=\mu+a_{i}+d_{i}+x_{i}+\varepsilon_{i},\ i=1,\dots,100,$ (5)

where $\mu$ is an intercept with default $\mathcal{N}(0,1000^{2})$ prior and
$\varepsilon_{i}$ is the residual effect, representing environmental noise.
Further, $a_{i}$, $d_{i}$ and $x_{i}$ are additive, dominance and epistasis
(additive-by-additive epistasis) effects, respectively. These three add up to
the genetic effect $g_{i}=a_{i}+d_{i}+x_{i}$. We assume that
$\bm{a}=(a_{1},\dots,a_{100})\sim\mathcal{N}_{100}(\bm{0},\sigma_{\mathrm{a}}^{2}\bm{A})$,
$\bm{d}=(d_{1},\dots,d_{100})\sim\mathcal{N}_{100}(\bm{0},\sigma_{\mathrm{a}}^{2}\bm{D})$,
and
$\bm{x}=(x_{1},\dots,x_{100})\sim\mathcal{N}_{100}(\bm{0},\sigma_{\mathrm{a}}^{2}\bm{X})$,
and we use a sum-to-zero constraint on all genetic effects. The covariance
matrices $\bm{A}$, $\bm{D}$ and $\bm{X}$ are computed from the single
nucleotide polymorphism (SNP) matrix with thousands of genetic markers, see
Hem et al.,  (2021) for details. This model has structured covariance
matrices, and we use the "generic0" latent model. "generic0" requires the
argument Cmatrix, which is the precision (inverse covariance) matrix
$\bm{Q}_{*}$ for the effect. With these data, we get the following formula:

R> formula <\- y ~\+ mc(a, model = "generic0", Cmatrix = Q_a, constr = TRUE)
+\+ mc(d, model = "generic0", Cmatrix = Q_d, constr = TRUE) +\+ mc(x, model =
"generic0", Cmatrix = Q_x, constr = TRUE) We go through the reasoning behind a
prior where we use all available prior knowledge, following the tree structure
in Table 2.

Tree structure | Parameters, priors
---|---
|  $a+d+x+\varepsilon$$\varepsilon$$a+d+x$$a$$d+x$$d$$x$
---
| $\sigma_{\mathrm{a+d+x}+\varepsilon}^{2}\sim\mathrm{Jeffreys^{\prime}}$
---
$\omega_{\frac{\mathrm{g}}{\mathrm{g}+\varepsilon}}\sim\mathrm{PC}_{\mathrm{0}}(0.25)$
$\omega_{\frac{\mathrm{a}}{\mathrm{g}}}\sim\mathrm{PC}_{\mathrm{M}}(0.85,0.8)$
$\omega_{\frac{\mathrm{d}}{\mathrm{d+x}}}\sim\mathrm{PC}_{\mathrm{M}}(0.67,0.8)$
Table 2:  Tree structures and the corresponding parameters for the genomic
example: $g_{i}=a_{i}+d_{i}+x_{i}$.

The expert in genetics has information on the heritability, which is the
amount of total variance attributed to the genetic effects and on the
distribution of the genetic effect $\bm{g}$ to the additive, dominance and
epistasis effects $\bm{a}$, $\bm{d}$ and $\bm{x}$. The expert says the
heritability $\omega_{\frac{\mathrm{g}}{\mathrm{g}+\varepsilon}}$ is around
0.25, and that we want to avoid overfitting, which leads us to
$\omega_{\frac{\mathrm{g}}{\mathrm{g}+\varepsilon}}\sim\mathrm{PC}_{\mathrm{0}}(0.25)$.
The additive, dominance and epistasis effects have according to the expert a
division of the genetic variance that is around (85, 10, 5)%, respectively. To
achieve this, we must use two dual-splits to decompose the genetic variation,
and do this by splitting off the additive effect first, with a
$\mathrm{PC}_{\mathrm{M}}(0.85,0.8)$ prior on
$\omega_{\frac{\mathrm{a}}{\mathrm{g}}}$ (the amount of genetic variance that
is additive). Then we attribute the remaining 15% of the genetic variance to
$\bm{d}$ and $\bm{x}$ with 67% to $\bm{d}$ with
$\mathrm{PC}_{\mathrm{M}}(0.67,0.8)$ on
$\omega_{\frac{\mathrm{d}}{\mathrm{d+x}}}$. We choose a concentration
parameter value of $0.8$ because the expert is quite sure about the (85, 10,
5)% division. This corresponds to having 75% of the density mass in the
interval $[\mathrm{logit}(m)-1,\mathrm{logit}(m)+1]$. The expert does not want
to use expert knowledge for the total variance
$\sigma_{\mathrm{a+d+x}+\varepsilon}^{2}$, so we use Jeffreys’ prior.

We have simulated a dataset following the description in Hem et al.,  (2021)
(see also Gaynor et al., , 2017; Selle et al., , 2019), using the R package
AlphaSimR (Faux et al., , 2016; Gaynor, , 2019). The source code for
simulating the dataset is available in the Supplemental Materials in Hem et
al.,  (2021) (Hem et al., , 2020). The dataset is included as wheat_data in
makemyprior. To incorporate the expert knowledge in a unified way, we first
scale the covariance matrices have typical variance equal to 1 (see Sørbye and
Rue,  (2017) for details), using the function scale_precmat() in makemyprior:

R> wheat_data_scaled <\- wheat_dataR> wheat_data_scaled$Q_a <\-
scale_precmat(wheat_data$Q_a)R> wheat_data_scaled$Q_d <\-
scale_precmat(wheat_data$Q_d)R> wheat_data_scaled$Q_x <\-
scale_precmat(wheat_data$Q_x) This model is implemented as follows:

R> prior <\- make_prior(formula, wheat_data_scaled, prior = list(\+ tree = "s1
= (d, x); s2 = (a, s1); s3 = (s2, eps)",\+ w = list(s1 = list(prior = "pcM",
param = c(0.67, 0.8)),\+ s2 = list(prior = "pcM", param = c(0.85, 0.8)),\+ s3
= list(prior = "pc0", param = 0.25)))) Note that we omit the specification of
the total variance prior, as we choose the default Jeffreys’ prior. It can be
specified with adding V = list(s3 = list(prior = "jeffreys")) to the list
provided to the prior argument. We now do inference on this model and plot the
results:

R> posterior <\- inference_stan(prior, iter = 15000, warmup = 5000,\+ seed =
1, init = "0", chains = 1)R> plot_posterior_stan(posterior, param = "prior",
prior = TRUE)

Tree structure: d_x = (d,x); a_d_x = (a,d_x); eps_a_d_x = (eps,a_d_x)Weight
priors:w[d/d_x] ~ PCM(0.67, 0.8)w[a/a_d_x] ~ PCM(0.85, 0.8)w[eps/eps_a_d_x] ~
PC1(0.75)Total variance priors:V[eps_a_d_x] ~ Jeffreys’SAMPLING FOR MODEL
’full_file’ NOW (CHAIN 1).Chain 1:Chain 1: Gradient evaluation took 0.000463
secondsChain 1: 1000 transitions using 10 leapfrog steps per transition
wouldtake 4.63 seconds.Chain 1: Adjust your expectations accordingly!Chain
1:Chain 1:Chain 1: Iteration: 1 / 15000 [ 0%] (Warmup)Chain 1: Iteration: 1500
/ 15000 [ 10%] (Warmup)Chain 1: Iteration: 3000 / 15000 [ 20%] (Warmup)Chain
1: Iteration: 4500 / 15000 [ 30%] (Warmup)Chain 1: Iteration: 5001 / 15000 [
33%] (Sampling)Chain 1: Iteration: 6500 / 15000 [ 43%] (Sampling)Chain 1:
Iteration: 8000 / 15000 [ 53%] (Sampling)Chain 1: Iteration: 9500 / 15000 [
63%] (Sampling)Chain 1: Iteration: 11000 / 15000 [ 73%] (Sampling)Chain 1:
Iteration: 12500 / 15000 [ 83%] (Sampling)Chain 1: Iteration: 14000 / 15000 [
93%] (Sampling)Chain 1: Iteration: 15000 / 15000 [100%] (Sampling)Chain
1:Chain 1: Elapsed Time: 30.2285 seconds (Warm-up)Chain 1: 74.8169 seconds
(Sampling)Chain 1: 105.045 seconds (Total)Chain 1: Figure 6 shows the prior
and posterior together on the parameterization of the prior (see Table 2).

Figure 6: Prior and posterior distribution of the random effect parameters for
the genomic selection example.

We see that we do not have enough data to estimate the variance proportion for
the additive and nonadditive genetic effects: The posterior distribution is
almost identical to the prior distribution. Hem et al.,  (2021) have conducted
an extensive simulation study on this and similar models. They saw a strong
need for robust prior distributions, which we also see in Figure 6, because
the nonadditive effects $\bm{d}$ and $\bm{x}$ are strongly confounded with the
environmental effect $\bm{\varepsilon}$, and the number of observations is
small compared to the number of genetic markers that needs to be estimated
(Sorensen and Gianola, , 2007).

####  5.1.2 Latin square experiment

We consider the latin square experiment in Example 2.2. In line with Fuglstad
et al.,  (2020, Section 5.2), we expand the model and assume the treatment
effect now consists of a smooth signal
$\bm{c}^{(1)}=(c_{1}^{(1)},\dots,c_{9}^{(1)})\sim(\bm{0},\sigma_{\mathrm{c}^{(1)}}^{2}\mathbf{Q}_{\mathrm{RW2}}^{-1})$
where $\sigma_{\mathrm{c}^{(1)}}^{2}$ is the variance and
$\mathbf{Q}_{\mathrm{RW2}}^{-1}$ is the covariance matrix describing the
intrinsic second-order random walk (Rue and Held, , 2005, Chapter 3), and
random noise
$\bm{c}^{(2)}=(c_{1}^{(2)},\dots,c_{9}^{(2)})\sim\mathcal{N}_{9}(\bm{0},\sigma_{\mathrm{c}^{(2)}}^{2}\mathbf{I}_{9})$.
We remove implicit intercept and linear effect by requiring
$\sum_{i=1}^{9}c_{i}^{(1)}=0$ and $\sum_{i=1}^{9}ic_{i}^{(1)}=0$. To simplify
the notation, we use $f_{i,j}=a_{i}+b_{j}+c_{k[i,j]}^{(1)}+c_{k[i,j]}^{(2)}$.

We show how to create the prior distributions in Table 3. We want to avoid
overfitting of the model, and use a prior with shrinkage towards the residuals
in the top split with a median giving 75% residual effect. We do not have any
preference for the attribution of the row, column and treatment effects, and
use an ignorant Dirichlet prior for the middle split. In the bottom split we
again we want to avoid overfitting, and use a prior with shrinkage towards the
unstructured treatment effect and a median corresponding to 75% unstructured
treatment effect. At last we do not want to say anything about the scale of
the total variance, and use the default Jeffreys’ prior.

Tree structure | Parameters, priors
---|---
|
$a+b+c^{(1)}+c^{(2)}+\varepsilon$$\varepsilon$$a+b+c^{(1)}+c^{(2)}$$a$$b$$c^{(1)}+c^{(2)}$$c^{(1)}$$c^{(2)}$
---
|
$\sigma_{\mathrm{a+b}+\mathrm{c}^{(1)}+\mathrm{c}^{(2)}+\varepsilon}^{2}\sim\mathrm{Jeffreys^{\prime}}$
---
$\omega_{\frac{f_{i}}{f_{i}+\varepsilon}}\sim\mathrm{PC}_{\mathrm{0}}(0.25)$
$\big{(}\omega_{\frac{\mathrm{a}}{f_{i}}},\omega_{\frac{\mathrm{b}}{f_{i}}},1-\omega_{\frac{\mathrm{a}}{f_{i}}}-\omega_{\frac{\mathrm{b}}{f_{i}}}\big{)}\sim\mathrm{Dirichlet}(3)$
$\omega_{\frac{\mathrm{c}^{(1)}}{\mathrm{c}^{(1)}+\mathrm{c}^{(2)}}}\sim\mathrm{PC}_{\mathrm{0}}(0.25)$
Table 3:  Tree structures and the corresponding parameters for the prior used
in the latin square model:
$f_{i,j}=a_{i}+b_{j}+c_{k[i,j]}^{(1)}+c_{k[i,j]}^{(2)}$.

The dataset used in this model is included in makemyprior as latin_data. It is
a simulated dataset, following the description of Fuglstad et al.,  (2020,
Section 5.2), where we have used
$\sigma_{\mathrm{a}}=\sigma_{\mathrm{b}}=\sigma_{\mathrm{c}^{(2)}}=\sigma_{\varepsilon}=0.1$
and true treatment effect $c_{i}^{(1)}=0.02\cdot((i-5)^{2}-20/3)$. lin in the
formula below is the linear effect of treatment $k[i,j]$ and both the
intercept $\mu$ and the coefficient $\beta$ has a default
$\mathcal{N}(0,\sigma=1000)$ prior. The following will fit this model and
produce the plots in Figure 7:

R> formula <\- y ~ lin + mc(row) + mc(col) + mc(iid) +\+ mc(rw2, model =
"rw2", constr = TRUE, lin_constr = TRUE)R> prior <\- make_prior(\+ formula,
latin_data,\+ prior = list(tree = "s1 = (rw2, iid);\+ s2 = (row, col, s1); s3
= (s2, eps)",\+ w = list(s1 = list(prior = "pc0", param = 0.25),\+ s2 =
list(prior = "dirichlet"),\+ s3 = list(prior = "pc0", param = 0.25))))R>
posterior <\- inference_stan(prior, iter = 15000, warmup = 5000,\+ seed = 1,
init = "0", chains = 1,\+ control = list(adapt_delta = 0.9))R>
plot_posterior_stan(posterior, param = "prior", prior = TRUE)

Tree structure: iid_rw2 = (iid,rw2); row_col_iid_rw2 =
(row,col,iid_rw2);eps_row_col_iid_rw2 = (eps,row_col_iid_rw2)Weight
priors:w[iid/iid_rw2] ~ PC1(0.75)(w[row/row_col_iid_rw2],
w[col/row_col_iid_rw2]) ~ Dirichlet(3)w[eps/eps_row_col_iid_rw2] ~
PC1(0.75)Total variance priors:V[eps_row_col_iid_rw2] ~ Jeffreys’SAMPLING FOR
MODEL ’full_file’ NOW (CHAIN 1).Chain 1:Chain 1: Gradient evaluation took
0.000167 secondsChain 1: 1000 transitions using 10 leapfrog steps per
transition wouldtake 1.67 seconds.Chain 1: Adjust your expectations
accordingly!Chain 1:Chain 1:Chain 1: Iteration: 1 / 15000 [ 0%] (Warmup)Chain
1: Iteration: 1500 / 15000 [ 10%] (Warmup)Chain 1: Iteration: 3000 / 15000 [
20%] (Warmup)Chain 1: Iteration: 4500 / 15000 [ 30%] (Warmup)Chain 1:
Iteration: 5001 / 15000 [ 33%] (Sampling)Chain 1: Iteration: 6500 / 15000 [
43%] (Sampling)Chain 1: Iteration: 8000 / 15000 [ 53%] (Sampling)Chain 1:
Iteration: 9500 / 15000 [ 63%] (Sampling)Chain 1: Iteration: 11000 / 15000 [
73%] (Sampling)Chain 1: Iteration: 12500 / 15000 [ 83%] (Sampling)Chain 1:
Iteration: 14000 / 15000 [ 93%] (Sampling)Chain 1: Iteration: 15000 / 15000
[100%] (Sampling)Chain 1:Chain 1: Elapsed Time: 30.845 seconds (Warm-up)Chain
1: 67.4428 seconds (Sampling)Chain 1: 98.2878 seconds (Total)Chain 1: Figure 7
shows the prior and posterior together on the parameterization of the prior.
The posterior distribution of the bottom split, attributing the treatment
effect to the random noise and smooth signal, is only slightly different from
the prior, indicating that there is no strong signal about the smooth
treatment effect in the data. By using a prior with shrinkage towards only
random noise treatment effect, we avoid overfitting. The model has learned
about the three other variance proportions, and we see that even though the
prior on the amount of total variance going to the residual effect has
shrinkage towards $1$, the model is not restricted by this (top right plot).

Note that for inference with INLA, we must implement the linear constraint
with extraconstr = list(A = matrix(1:9, 1, 9), e = matrix(0, 1, 1)) in mc(rw2,
...) in the formula.

Figure 7: Prior and posterior distribution of the random effect parameters for
the latin square example.

###  5.2 Binomial responses: Neonatal mortality

This example is based on a study carried out by Fuglstad et al.,  (2020).
Child mortality is an important indicator of health and well-being in a
country. We define neonatal mortality as the number of deaths of infants the
first month of life per live birth, which can be estimated using national
household surveys from Demographic and Health Surveys (Kenya National Bureau
of Statistics et al., , 2015). From such surveys we can extract the number of
live births $b_{i,j}$ and the number of neonatal deaths $y_{i,j}$ in cluster
$j$ in county $i$, and use an indicator $x_{i,j}$ for classifying cluster $j$
in county $i$ as rural ($x_{i,j}=0)$ or urban ($1$). We model
$y_{i,j}|b_{i,j},p_{i,j}\sim\mathrm{Binomial}(b_{i,j},p_{i,j})$ with 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},$ (6)

where $v_{i}\sim\mathcal{N}(0,\sigma_{v}^{2})$ and
$\nu_{i,j}\sim\mathcal{N}(0,\sigma_{\nu}^{2})$ are i.i.d. random effects with
sum-to-zero constraints for counties and clusters, respectively, and $\bm{u}$
is a Besag effect on county with variance $\sigma_{u}^{2}$ and a sum-to-zero
constraint. In the Besag model, the spatial effect of each county depends on
the effects in the neighboring regions (see e.g. Besag et al.,  (1991) for
details), and when combining it with an i.i.d. effect on the same level in the
hierarchy, we get a BYM (Besag, York and Mollié) model (Besag et al., , 1991).
We want to investigate whether or not there is a spatial effect present.

We simulated a dataset following the description in Fuglstad et al.,  (2020,
Section 6.2) with the 47 counties in Kenya (see Figure 9 for a map). We used
6, 7 or 8 clusters in each county which gave in total $327$ clusters, and thus
$327$ observations, $b_{i,j}=25$ live births in each cluster, and parameters
$\mu=-4$, $\beta=0.1$, $\sigma_{\nu}^{2}=0.2$, $\sigma_{\mathrm{v}}^{2}=0.1$,
and $\sigma_{\mathrm{u}}^{2}=0.5$.

This dataset is available in makemyprior as neonatal_data, as well as other
necessary files for fitting the model.

Tree structure | Parameters, priors
---|---
|  $u+v+\nu$$u+v$$u$$v$$\nu$
---
| $\sigma_{\mathrm{u+v}+\nu}^{2}\sim\mathrm{PC}_{0}(3.35,0.05)$
---
$\omega_{\frac{\mathrm{u+v}}{\mathrm{u+v}+\nu}}\sim\mathrm{PC}_{\mathrm{1}}(0.75)$
$\omega_{\frac{\mathrm{u}}{\mathrm{u+v}}}\sim\mathrm{PC}_{\mathrm{0}}(0.25)$
Table 4:  Tree structures and the corresponding parameters for the neonatal
mortality model.

We prefer coarser over finer unstructured effects, and unstructured over
structured effects. That means that we prefer $\bm{v}$ over $\bm{u}$ and
$\bm{v}+\bm{u}$ over $\bm{\nu}$ in the prior. The BYM model is intuitively
represented with a dual split in the prior tree, where one leaf node
represents a Besag effect and the other represents an i.i.d. effect. We
achieve this with a prior that distributes the between-county variance with
shrinkage towards the unstructured county effect, which gives the BYM2 model
of Riebler et al.,  (2016), and we shrink the total variance towards the
county effects. Following Fuglstad et al.,  (2020), 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:

R> set.seed(1)R> find_pc_prior_param(lower = 0.1, upper = 10, prob = 0.9, N =
2e5)

U = 3.353132Prob(0.09866969 < exp(eta) < 9.892902) = 0.9 This gives a
$\mathrm{PC}_{\mathrm{0}}(3.35,0.05)$ prior. The tree structure and a summary
of the prior distributions can be found in Table 4. We fit the model with
Stan:

R> graph_path <\- paste0(path.package("makemyprior"), "/neonatal.graph")R>
formula <\- y ~ mc(nu) + mc(v) +\+ mc(u, model = "besag", graph = graph_path,
scale.model = TRUE)R> prior <\- 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)))))R>
posterior <\- inference_stan(prior, iter = 15000, warmup = 5000,\+ seed = 1,
init = "0", chains = 1,\+ control = list(adapt_delta = 0.85))

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)SAMPLING FOR MODEL ’full_file’ NOW (CHAIN 1).Chain 1:Chain 1:
Gradient evaluation took 0.000511 secondsChain 1: 1000 transitions using 10
leapfrog steps per transitionwould take 5.11 seconds.Chain 1: Adjust your
expectations accordingly!Chain 1:Chain 1:Chain 1: Iteration: 1 / 15000 [ 0%]
(Warmup)Chain 1: Iteration: 1500 / 15000 [ 10%] (Warmup)Chain 1: Iteration:
3000 / 15000 [ 20%] (Warmup)Chain 1: Iteration: 4500 / 15000 [ 30%]
(Warmup)Chain 1: Iteration: 5001 / 15000 [ 33%] (Sampling)Chain 1: Iteration:
6500 / 15000 [ 43%] (Sampling)Chain 1: Iteration: 8000 / 15000 [ 53%]
(Sampling)Chain 1: Iteration: 9500 / 15000 [ 63%] (Sampling)Chain 1:
Iteration: 11000 / 15000 [ 73%] (Sampling)Chain 1: Iteration: 12500 / 15000 [
83%] (Sampling)Chain 1: Iteration: 14000 / 15000 [ 93%] (Sampling)Chain 1:
Iteration: 15000 / 15000 [100%] (Sampling)Chain 1:Chain 1: Elapsed Time:
73.7369 seconds (Warm-up)Chain 1: 145.026 seconds (Sampling)Chain 1: 218.763
seconds (Total)Chain 1: For inference with INLA, the Ntrials argument must be
provided to inference_inla(). The following produce the plots in Figure 8 and
gives some key information on the posterior:

R> plot_posterior_fixed(posterior)R> plot_posterior_stan(posterior, param =
"prior", prior = TRUE)R> posterior

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)Inference
done with Stan. Param. mean median sd V[nu_v_u] 0.662 0.630 0.244 w[v/v_u]
0.655 0.715 0.263 w[nu/nu_v_u] 0.333 0.330 0.190 intercept -4.157 -4.153 0.135
urban 0.469 0.469 0.169 Figure 9 shows the posterior spatial effect
$e^{u_{i}}$ plotted in a map. We see a spatial variation between the counties.
The necessary data for creating the spatial map are not included, but can be
obtained from https://gadm.org/. The samples for the effects can easily be
extracted with extract_posterior_effect() which take the arguments obj from
inference_stan() and the name of the effect:

R> u <\- extract_posterior_effect(posterior, "u")

(a) Posterior of the effect of urban/rural and the intercept.

(b) Random effect parameters.

Figure 8: Prior and posterior distribution of 8(a)) coefficients of the fixed
effects and 8(b)) total variance and variance proportions of the random
effects for the neonatal mortality example.

The fixed effects in Figure 8(a) show that the intercept is not contributing
much to the linear predictor, while the effect of urban/rural shows that there
is a higher mortality in urban areas (which is the case also for real data,
see Kenya National Bureau of Statistics et al.,  (2015)). From Figure 8(b) we
see that the model has learned about the total variance from the data and
about the amount of total (latent) variance to the cluster effect
($\bm{\nu}$), but there is not enough information in the data about the amount
of county variance to the structured county effect ($\bm{u}$). The following
fits the model without the likelihood (sampling from the prior) and produces
the plots of the prior and posterior on standard deviation scale in Figure 10:

Figure 9: Posterior median of $e^{u_{i}}$ for each county in Kenya. Note that
this is based on simulated data. Figure 10: Prior and posterior distribution
for the neonatal mortality example on standard deviation scale.

R> prior_samps <\- inference_stan(prior, use_likelihood = F, print_prior =
F,\+ iter = 15000, warmup = 5000,\+ seed = 1, init = "0", chains = 1)R>
plot_several_posterior_stan(list(Prior = prior_samps,\+ Posterior =
posterior), "stdev")

SAMPLING FOR MODEL ’full_file’ NOW (CHAIN 1).Chain 1:Chain 1: Gradient
evaluation took 0.000119 secondsChain 1: 1000 transitions using 10 leapfrog
steps per transition wouldtake 1.19 seconds.Chain 1: Adjust your expectations
accordingly!Chain 1:Chain 1:Chain 1: Iteration: 1 / 15000 [ 0%] (Warmup)Chain
1: Iteration: 1500 / 15000 [ 10%] (Warmup)Chain 1: Iteration: 3000 / 15000 [
20%] (Warmup)Chain 1: Iteration: 4500 / 15000 [ 30%] (Warmup)Chain 1:
Iteration: 5001 / 15000 [ 33%] (Sampling)Chain 1: Iteration: 6500 / 15000 [
43%] (Sampling)Chain 1: Iteration: 8000 / 15000 [ 53%] (Sampling)Chain 1:
Iteration: 9500 / 15000 [ 63%] (Sampling)Chain 1: Iteration: 11000 / 15000 [
73%] (Sampling)Chain 1: Iteration: 12500 / 15000 [ 83%] (Sampling)Chain 1:
Iteration: 14000 / 15000 [ 93%] (Sampling)Chain 1: Iteration: 15000 / 15000
[100%] (Sampling)Chain 1:Chain 1: Elapsed Time: 6.10027 seconds (Warm-up)Chain
1: 18.8748 seconds (Sampling)Chain 1: 24.9751 seconds (Total)Chain 1: From
these graphs we see that the posterior of the standard deviations are clearly
different from the prior. We saw in Figure 8(b) that the model did not learn
much about the amount of county variation accounted for by the Besag effect
($\bm{u}$), but we cannot see this from plots of the posterior standard
deviations, and they do not show the whole picture. This is another advantage
of the HD prior; it is easy to see that even though we get the impression that
the model has learned from the data, that knowledge is not necessarily about
the whole model. This shows, as Fuglstad et al.,  (2020) points out, that one
should be careful before drawing conclusions on first impressions about the
results, and more investigation should be done.

##  6 Summary and discussion

The makemyprior package offers an intuitive and transparent way of choosing
and visualizing prior distributions for Bayesian hiearchical models. It is
easy to utilize expert knowledge, and clear what prior distributions are used,
also when the default settings are chosen. The package works with the flexible
and widely used latent Gaussian models, and offers Gaussian, binomial and
Poisson likelihoods. This, together with the latent models included, makes the
package applicable in a range of problems and applications.

The package can be used to investigate the prior choices, and the user can
simulate from the prior with Stan and look at the prior distributions on
different parameterizations. In this way, crucial misunderstandings of what
prior distributions are used can be discovered and corrected, and thus
increase the understanding and meaning of the prior. The usage of the
hierarchical decomposition (HD) prior is not limited to inference carried out
with rstan or INLA, and the prior can be computed with makemyprior and fed
into other software programs for Bayesian inference.

To see how the individual fixed effects contribute to the total data variation
would be interesting, and could be done with the HD prior framework. However,
fixed effects are often correlated, and the variance that is explained by each
single fixed effect is not well defined. The perhaps most intuitive way to
include fixed effects directly in the tree structure is to assign one variance
parameter to each effect, but this can quickly increase the amount of variance
parameters to a level where inference become computationally hard. Gelman and
Hill,  (2007) and Zhang et al.,  (2020) have proposed prior distributions
related to the coefficient of determination, $R^{2}$, which measures the
amount of variance explained by the model. The generalized $R^{2}$ proposed by
Gelman and Hill,  (2007) measures this at each level in the hierarchical
model. To include fixed effects directly in the joint prior in the HD prior
framework is discussed by Fuglstad et al.,  (2020). We have chosen to give
them independent and vague priors.

Other exiting additions to the HD prior framework include extending it to
models outside the class of latent Gaussian models, or to models where the
hyperparameters of the priors will get prior distributions. This will require
further development of the framework, and may be highly computational
expensive, as the penalized complexity (PC) prior cannot be pre-computed in
the same way as we do now with the conditioning on the hyperparameters.

To open for using Dirichlet distributions with custom hyperparameter values
can be a natural next step. This addition to the HD prior framework itself
will open for even easier integration of the HD prior in other software (such
as Template Model Builder (TMB, Kristensen et al., , 2016)) as expert
knowledge can be included without having to compute the PC prior. However,
this will complicate the intuition behind the prior, and it will be more
difficult to use prior and expert knowledge in a transparent way. It will
require more thoughtful prior choices, and we lose one of the big advantages
with the easy-to-use and intuitive way of making priors with makemyprior, in
addition to the shrinkage properties of the PC prior. Allowing the user to
specify custom prior distributions for variance parameters will introduce even
more flexibility, however, the included variance prior distributions cover the
most popular choices.

Including more latent models will further increase the amount of applications
the package can be used for without specifying custom Stan-code. This includes
handling parameters such as correlations, which can be done by giving them
independent priors and conditioning on a representative value (e.g. the mean
or median of the chosen prior) for this parameter when the joint prior is
computed. To open for easy integration into other software for inference, such
as TMB, can be useful for models that are very complex and will be highly time
consuming and difficult to fit with rstan or INLA.

In conclusion, makemyprior offers something not offered by the range of
packages that can be used to carry out inference for Bayesian hierarchical
models. It makes it easy to include prior knowledge in an intuitive and
transparent way, can be used to verify prior choices, and allows direct
inference in a simple way. makemyprior makes users aware of what priors are
used and makes prior selection a concious choice, which is important when
doing inference to ensure that the model fitted is indeed the intended one.

## Computational details

The results in this paper were obtained with R version 4.0.2 on platform
x86_64-apple-darwin17.0 (64-bit) running under macOS Catalina 10.15.7.

Package versions: makemyprior 1.0.0, ggplot2 3.3.2, Matrix 1.2.18, knitr 1.33,
shiny 1.5.0, shinyjs 2.0.0, shinyBS 0.61, visNetwork 2.0.9, rmarkdown 2.3,
rlang 0.4.11, MASS 7.3.51.6, testthat 2.3.2, rstan 2.21.2, and INLA 20.3.17.

## Acknowledgments

Hem, Fuglstad and Riebler were supported by project number 240873 from the
Research Council of Norway.

## References

  * Besag et al., (1991) Besag, J., York, J., and Mollié, A. (1991).  Bayesian Image Restoration, With Two Applications in Spatial Statistics.  Annals of the Institute of Statistical Mathematics, 43(1):1–20.
  * Caldwell et al., (2009) Caldwell, A., Kollar, D., and Kröninger, K. (2009).  BAT–The Bayesian analysis toolkit.  Computer Physics Communications, 180(11):2197–2209.
  * Carpenter et al., (2017) Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P., and Riddell, A. (2017).  Stan: A probabilistic programming language.  Journal of Statistical Software, 76(1).
  * Chang et al., (2020) Chang, W., Cheng, J., Allaire, J., Xie, Y., and McPherson, J. (2020).  shiny: Web Application Framework for R.  R package version 1.5.0.
  * Depaoli et al., (2020) Depaoli, S., Winter, S. D., and Visser, M. (2020).  The Importance of Prior Sensitivity Analysis in Bayesian Statistics: Demonstrations Using an Interactive Shiny app.  Frontiers in Psychology, 11.
  * Faux et al., (2016) Faux, A.-M., Gorjanc, G., Gaynor, R. C., Battagin, M., Edwards, S. M., Wilson, D. L., Hearne, S. J., Gonen, S., and Hickey, J. M. (2016).  AlphaSim: Software for Breeding Program Simulation.  The Plant Genome, 9(3):1–14.
  * Fong et al., (2010) Fong, Y., Rue, H., and Wakefield, J. (2010).  Bayesian Inference for Generalized Linear Mixed Models.  Biostatistics, 11(3):397–412.
  * Fuglstad et al., (2020) Fuglstad, G.-A., Hem, I. G., Knight, A., Rue, H., and Riebler, A. (2020).  Intuitive Joint Priors for Variance Parameters.  Bayesian Analysis, 15(4):1109–1137.
  * Gabry, (2018) Gabry, J. (2018).  shinystan: Interactive Visual and Numerical Diagnostics and Posterior Analysis for Bayesian Models.  R package version 2.5.0.
  * Gaynor, (2019) Gaynor, C. (2019).  AlphaSimR: Breeding Program Simulations.  R package version 0.10.0.
  * Gaynor et al., (2017) Gaynor, R. C., Gorjanc, G., Bentley, A. R., Ober, E. S., Howell, P., Jackson, R., Mackay, I. J., and Hickey, J. M. (2017).  A Two-Part Strategy for Using Genomic Selection to Develop Inbred Lines.  Crop Science, 57(5):2372–2386.
  * Gelman and Hill, (2007) Gelman, A. and Hill, J. (2007).  Data Analysis Using Regression and Multilevel/Hierarchical Models, volume 1.  Cambridge University Press, New York, New York.
  * Gelman et al., (2017) Gelman, A., Simpson, D., and Betancourt, M. (2017).  The Prior Can Often Only Be Understood in the Context of the Likelihood.  Entropy, 19(10):555.
  * Gelman et al., (2020) Gelman, A., Vehtari, A., Simpson, D., Margossian, C. C., Carpenter, B., Yao, Y., Kennedy, L., Gabry, J., Bürkner, P.-C., and Modrák, M. (2020).  Bayesian Workflow.  arXiv preprint arXiv:2011.01808 [stat.ME].
  * Goel and Degroot, (1981) Goel, P. K. and Degroot, M. H. (1981).  Information About Hyperparameters in Hierarchical Models.  Journal of the American Statistical Association, 76(373):140–147.
  * Goodrich et al., (2020) Goodrich, B., Gabry, J., Ali, I., and Brilleman, S. (2020).  rstanarm: Bayesian Applied Regression Modeling via Stan.  R package version 2.21.1.
  * Guo et al., (2017) Guo, J., Riebler, A., and Rue, H. (2017).  Bayesian bivariate meta-analysis of diagnostic test studies with interpretable priors.  Statistics in Medicine, 36(19):3039–3058.
  * Hem et al., (2020) Hem, I. G., Selle, M., Gorjanc, G., Fuglstad, G.-A., and Riebler, A. (2020).  Supplemental Material for ”Robust Modelling of Additive and Non-Additive Variation With Intuitive Inclusion of Expert Knowledge”.
  * Hem et al., (2021) Hem, I. G., Selle, M. L., Gorjanc, G., Fuglstad, G.-A., and Riebler, A. (2021).  Robust Modeling of Additive and Nonadditive Variation With Intuitive Inclusion of Expert Knowledge.  Genetics, 217(3).  iyab002.
  * Hoffman and Gelman, (2014) Hoffman, M. D. and Gelman, A. (2014).  The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo.  Journal of Machine Learning Research, 15(1):1593–1623.
  * Holand et al., (2013) Holand, A. M., Steinsland, I., Martino, S., and Jensen, H. (2013).  Animal Models and Integrated Nested Laplace Approximations.  G3: Genes, Genomes, Genetics, 3(8):1241–1251.
  * Kenya National Bureau of Statistics et al., (2015) Kenya National Bureau of Statistics, Ministry of Health/Kenya, National AIDS Control Council/Kenya, Kenya Medical Research Institute, and National Council for Population and Development/Kenya (2015).  Kenya Demographic and Health Survey 2014.  Rockville, MD, USA.
  * Kristensen et al., (2016) Kristensen, K., Nielsen, A., Berg, C. W., Skaug, H., and Bell, B. M. (2016).  TMB: Automatic Differentiation and Laplace Approximation.  Journal of Statistical Software, 70(5):1–21.
  * Lambert et al., (2005) Lambert, P. C., Sutton, A. J., Burton, P. R., Abrams, K. R., and Jones, D. R. (2005).  How Vague is Vague? A Simulation Study of the Impact of the Use of Vague Prior Distributions in MCMC using WinBUGS.  Statistics in Medicine, 24(15):2401–2428.
  * Plummer, (2017) Plummer, M. (2017).  JAGS Version 4.3. 0 User Manual [Computer Software Manual].  Retrieved from: sourceforge.net/projects/mcmc-jags/files/Manuals/4.x.
  * Riebler et al., (2016) Riebler, A., Sørbye, S. H., Simpson, D., and Rue, H. (2016).  An Intuitive Bayesian Spatial Model for Disease Mapping That Accounts for Scaling.  Statistical Methods in Medical Research, 25(4):1145–1165.
  * Rue and Held, (2005) Rue, H. and Held, L. (2005).  Gaussian Markov Random Fields: Theory and Applications.  CRC press, Boca Raton, Florida.
  * Rue et al., (2009) Rue, H., Martino, S., and Chopin, N. (2009).  Approximate Bayesian Inference for Latent Gaussian Models by Using Integrated Nested Laplace Approximations.  Journal of the Royal Statistical Society B (Statistical Methodology), 71(2):319–392.
  * Rue et al., (2017) Rue, H., Riebler, A., Sørbye, S. H., Illian, J. B., Simpson, D. P., and Lindgren, F. K. (2017).  Bayesian Computing With INLA: A Review.  Annual Review of Statistics and Its Application, 4:395–421.
  * Selle et al., (2019) Selle, M. L., Steinsland, I., Hickey, J. M., and Gorjanc, G. (2019).  Flexible Modelling of Spatial Variation in Agricultural Field Trials with the R Package INLA.  Theoretical and Applied Genetics, 132(12):3277–3293.
  * Simpson et al., (2017) Simpson, D., Rue, H., Riebler, A., Martins, T. G., and Sørbye, S. H. (2017).  Penalising Model Component Complexity: A Principled, Practical Approach to Constructing Priors.  Statistical Science, 32(1):1–28.
  * Smid and Winter, (2020) Smid, S. C. and Winter, S. D. (2020).  Dangers of the Defaults: A Tutorial on the Impact of Default Priors When Using Bayesian SEM With Small Samples.  Frontiers in Psychology, 11:3536.
  * Sørbye and Rue, (2017) Sørbye, S. H. and Rue, H. (2017).  Penalised Complexity Priors for Stationary Autoregressive Processes.  Journal of Time Series Analysis, 38(6):923–935.
  * Sorensen and Gianola, (2007) Sorensen, D. and Gianola, D. (2007).  Likelihood, Bayesian, and MCMC methods in Quantitative Genetics.  Springer Science & Business Media.
  * Stan Development Team, (2020) Stan Development Team (2020).  RStan: the R interface to Stan.  R package version 2.21.2.
  * Vehtari et al., (2020) Vehtari, A., Gabry, J., Magnusson, M., Yao, Y., Bürkner, P.-C., Paananen, T., and Gelman, A. (2020).  loo: Efficient leave-one-out cross-validation and WAIC for Bayesian models.  R package version 2.4.1.
  * Wakefield, (2006) Wakefield, J. (2006).  Disease Mapping and Spatial Regression With Count Data.  Biostatistics, 8(2):158–183.
  * Zhang et al., (2020) Zhang, Y. D., Naughton, B. P., Bondell, H. D., and Reich, B. J. (2020).  Bayesian Regression Using a Prior on the Model Fit: The R2-D2 Shrinkage Prior.  Journal of the American Statistical Association, 0(0):1–13.
  * Zondervan-Zwijnenburg et al., (2017) Zondervan-Zwijnenburg, M., Peeters, M., Depaoli, S., and de Schoot, R. V. (2017).  Where Do Priors Come From? Applying Guidelines to Construct Informative Priors in Small Sample Research.  Research in Human Development, 14(4):305–320.

##  Appendix A Code for Figure 2 in Section 3

formula <\- y ~ -1 + mc(a)data <\- list( a = rep(1:10, 10), y = rep(0,
100))prior1 <\- make_prior( formula, data, prior = list(tree = "(a);
(eps)"))prior2 <\- make_prior( formula, data, prior = list(tree = "s1 = (a,
eps)", w = list(s1 = list(prior = "pc0", param = 0.25))))prior3 <\-
make_prior( formula, data, prior = list(tree = "s1 = (a, eps)", w = list(s1 =
list(prior = "pc1", param = 0.75))))prior4 <\- make_prior( formula, data,
prior = list(tree = "s1 = (a, eps)", w = list(s1 = list(prior = "pcM", param =
c(0.25, 0.85)))))plot_marginal_prior(seq(0, 5, 0.050), prior1, "sigma^2[a]",
sd = T)plot_marginal_prior(seq(0, 1, 0.001), prior2,
"w[a/a_eps]")plot_marginal_prior(seq(0, 1, 0.001), prior3,
"w[a/a_eps]")plot_marginal_prior(seq(0, 1, 0.001), prior4, "w[a/a_eps]")