Example i.i.d. model

In makemyprior: Intuitive Construction of Joint Priors for Variance Parameters

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
options(rmarkdown.html_vignette.check_title = FALSE)

library(makemyprior)

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

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

Prior options

Before going into the details on how to use the software, we give you a short introduction to the HD prior. We also refer to @fuglstad2020 for details.

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

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

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

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

Model description

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

Make data and linear predictor

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

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

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

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

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

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

Make prior

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

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

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

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

summary(prior)

plot_prior(prior) # or plot(prior)
plot_tree_structure(prior)

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

new_prior <- makemyprior_gui(prior)

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

With the following command, we specify this prior:

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

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

summary(new_prior)
plot_prior(new_prior)
plot_tree_structure(new_prior)

Inference

We can carry out inference with Stan [@carpenter2017] and INLA [@rue2009]. Note that we in this vignette do not run the inference, as it takes time and will slow down the compilation of the vignette and thus the package download, but the code is included below and the user can carry out the inference with that.

Stan

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

compile_stan(save = TRUE)

Then we can do the inference:

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

We can look at the graphs of the posterior:

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

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

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

INLA

Inference with INLA is carried out in a similar way:

posterior2 <- inference_inla(new_prior)

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

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

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

Default settings {#default}

• 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 variance in the prior.
• The prior on the total variance (top nodes) 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 individual variance (singletons) 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.

See:

?makemyprior_models

for details on default settings.

Additional examples

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

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

plot_prior(prior2)
plot_tree_structure(prior2)

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

plot_prior(prior3)
plot_tree_structure(prior3)

sessionInfo()

References

makemyprior documentation built on Sept. 11, 2024, 6:56 p.m.