Title:	Intuitive Construction of Joint Priors for Variance Parameters
Version:	1.2.2
Description:	Tool for easy prior construction and visualization. It helps to formulates joint prior distributions for variance parameters in latent Gaussian models. The resulting prior is robust and can be created in an intuitive way. A graphical user interface (GUI) can be used to choose the joint prior, where the user can click through the model and select priors. An extensive guide is available in the GUI. The package allows for direct inference with the specified model and prior. Using a hierarchical variance decomposition, we formulate a joint variance prior that takes the whole model structure into account. In this way, existing knowledge can intuitively be incorporated at the level it applies to. Alternatively, one can use independent variance priors for each model components in the latent Gaussian model. Details can be found in the accompanying scientific paper: Hem, Fuglstad, Riebler (2024, Journal of Statistical Software, <doi:10.18637/jss.v110.i03>).
License:	GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
URL:	https://github.com/ingebogh/makemyprior
BugReports:	https://github.com/ingebogh/makemyprior/issues
Encoding:	UTF-8
LazyData:	true
Imports:	ggplot2, Matrix, methods, shiny, shinyjs, shinyBS, visNetwork, rlang, MASS
Suggests:	rstan (≥ 2.21.2), INLA (≥ 20.03.17), knitr, rmarkdown
Depends:	R (≥ 3.6.0)
Additional_repositories:	https://inla.r-inla-download.org/R/stable/
VignetteBuilder:	knitr
RoxygenNote:	7.2.3
NeedsCompilation:	no
Packaged:	2024-08-23 10:35:16 UTC; ingebogh
Author:	Ingeborg Hem [cre, aut], Geir-Arne Fuglstad [aut], Andrea Riebler [aut]
Maintainer:	Ingeborg Hem <ingeborg.hem@ntnu.no>
Repository:	CRAN
Date/Publication:	2024-08-23 11:10:06 UTC

Compile stan-model

Description

Function that compiles the stan-code for inference that is included in the model. The compiled version is stored in a .rds-file, which is by default stored in tempdir(). Can also be stored in the package (permanent = TRUE), or in a custom directory. In the latter case, this custom directory must be specified every time inference_stan is called. This will also be done by inference_stan, but this way can save some time if it is not already pre-compiled.

Usage

compile_stan(save = FALSE, permanent = FALSE, path = NULL)

Arguments

Details

Note that you will get a message saying something about integer division. The PC priors on variance proportions are represented by splines, and to evaluate them in Stan we look up values, and use integer division for this. This does not cause problems.

Value

Returns the stan-model invisibly.

Examples

if (interactive() && requireNamespace("rstan")){
  compile_stan(save = TRUE) # saving in tempdir()
}

Create a "skeleton" for custom Stan code

Description

Makes and saves files with generic code for writing custom Stan code and still use the HD prior.

Usage

create_stan_file(save = FALSE, location = "")

Arguments

Details

Must be in an interactive session to store the code. A folder called "my_stan_code" will be created in location. If the folder already exists in the specified location, you get an error. The folder contains:

main_file.stan

Main file. Can put all necessary functions here, but for a cleaner code that is easier to read, we put the functions in separate files.

jacobian.stan

Function that automatically computes the Jacobian, needed to transform from weights and total variance parameterization to log-variance parameterization.

prior_distributions.stan

Functions for computing the prior distributions.

The provided code is written so a random intercept model with an intercept, a group effect and a residual effect can be fitted:

example_custom_stan.R

R script showing how one can fit a random intercept model using the provided code.

The code can be expanded to fit the desired model. This requires some knowledge about Stan. No more documentation is given, as this is merely an offer to users who want to use other models than what are provided in the package already, and will be highly model specific.

Value

Nothing.

Examples

## Not run: 
create_stan_file(TRUE, "")

## End(Not run)

Evaluate the joint variance prior

Description

Function for evaluating the joint variance prior stored in prior_obj. To compute the joint prior, the functions needs to know the transformation from the total variance/variance proportion scale to log-variance scale. This is computed before inference, but is not stored in the mmp_prior-object. To avoid having to recompute this for every evaluation and thus improve the speed, we make a condensed data object with the function make_eval_prior_data.

Usage

eval_joint_prior(theta, prior_data)

make_eval_prior_data(prior_obj)

Arguments

Details

Note that a Jeffreys' prior is improper and sampling with the prior only will not work when it is used. For sampling from the prior (for example for debugging), use a proper prior for all parameters instead.

The make_eval_prior_data function is used to create a condensed version of the prior object from make_prior, that only contains what is needed to compute the joint prior. Since the HD prior is chosen on total variances and variance proportions, some additional information is needed to compute the Jacobian for the joint prior. To improve the speed, we do this once before evaluating the prior.

Expert option: make_eval_prior_data can also be used to extract the prior to be used with 'regular' INLA. See examples for how this can be done.

Value

Logarithm of the prior density.

Examples

ex_model <- makemyprior_example_model()
get_parameter_order(ex_model) # a, b, eps
prior_data <- make_eval_prior_data(ex_model)
eval_joint_prior(c(0, 0, 0), prior_data)
eval_joint_prior(c(-1, 0, 1), prior_data)

# a model with only 2 variance parameters
if (interactive()){

  data <- list(
    a = rep(1:10, each = 10)
  )
  set.seed(1); data$y <- rnorm(10, 0, 0.4)[data$a] + rnorm(100, 0, 1)

  # random intercept model
  ex_model2 <- make_prior(y ~ mc(a), data, family = "gaussian",
                          prior = list(tree = "s2 = (a, eps)",
                                       w = list(s2 = list(prior = "pc0", param = 0.25)),
                                       V = list(s2 = list(prior = "pc", param = c(3, 0.05)))),
                          intercept_prior = c(0, 1000))

  prior_data2 <- make_eval_prior_data(ex_model2)
  # evaluating the prior in a grid
  theta_a <- seq(-8, 4, 0.1)
  theta_eps <- seq(-8, 4, 0.1)
  res <- matrix(nrow = 0, ncol = 3)
  for (ind in 1:length(theta_a)){
    for (jnd in 1:length(theta_eps)){
      res <- rbind(res, c(theta_a[ind], theta_eps[jnd],
                          eval_joint_prior(c(theta_a[ind], theta_eps[jnd]), prior_data2)))
    }
  }

  # graph showing the prior
  if (requireNamespace("ggplot2")){
    res2 <- as.data.frame(res)
    names(res2) <- c("x", "y", "z")
    # Note from the "exp(z)" that we use the posterior, and not log posterior, in this plot
    ggplot(res2, aes(x = x, y = y, z = exp(z), fill = exp(z))) +
      geom_raster() +
      geom_contour(color = "black") +
      scale_fill_viridis_c(option = "E") +
      xlab("Log variance of 'a'") +
      ylab("Log residual variance") +
      labs(fill = "Density") +
      theme_bw()
  }

}

## Not run: 

# How an HD prior can be computed with \code{makemyprior}, and then sent to regular INLA 
# (expert option).
# Note the use of the hidden \code{make_jpr}-function.
# Also note that the order of the parameters must be the same as in the call to \code{make_prior}.
# The residual variance is put in the correct place by \code{make_jpr}.
data <- list(
  a = rep(1:10, each = 100),
  b = rep(1:100, times = 10)
)
set.seed(1); data$y <- rnorm(100, 0, 0.4)[data$a] + rnorm(100, 0, 0.6)[data$b] + rnorm(1000, 0, 1)
prior <- make_prior(y ~ mc(a) + mc(b), data, family = "gaussian",
                    prior = list(tree = "s1 = (a, b); s2 = (s1, eps)",
                                 w = list(s2 = list(prior = "pc0", param = 0.25)),
                                 V = list(s2 = list(prior = "pc", param = c(3, 0.05)))),
                    intercept_prior = c(0, 1000))
jpr_dat <- make_eval_prior_data(prior)
res <- inla(y ~ f(a) + f(b),
            data = data,
            control.fixed = list(prec.intercept = 1/1000^2),
            control.expert = list(jp = makemyprior:::make_jpr(jpr_dat)))

## End(Not run)

Evaluate PC prior for variance proportion

Description

Evaluate PC prior for a variance proportion.

Usage

eval_pc_prior(x, obj, param, logitscale = FALSE)

Arguments

Value

Returns density for the given variance proportion.

Examples

ex_prior <- makemyprior_example_model()
eval_pc_prior(seq(0, 1, 0.01), ex_prior, "eps/eps_a_b")
# or:
eval_pc_prior(seq(0, 1, 0.01), ex_prior, "w[eps/eps_a_b]")

expit

Description

Calculates inverse logit, exp(x)/(1+exp(x)) = 1/(1+exp(-x))

Usage

expit(x)

Arguments

Value

expit value

Examples

expit(2)
expit(seq(-5, 5, 1))

Extract the posterior of a random effect

Description

Extract the posterior of a random effect in the model for inference done with Stan

Usage

extract_posterior_effect(obj, effname)

Arguments

Value

Returns a matrix with the posterior samples of the chosen effect

Examples

if (interactive() && requireNamespace("rstan")){
  ex_prior <- makemyprior_example_model()
  res_stan <- inference_stan(ex_prior, iter = 100)
  # Note: For reliable results, increase the number of iterations (e.g., 'iter = 2000')
  extract_posterior_effect(res_stan, "a")
}

Extract the posterior parameter estimate

Description

Extract the posterior parameter estimate of a random effect variance or fixed effect coefficient when inference is done with Stan

Usage

extract_posterior_parameter(obj, param)

Arguments

Value

Returns a vector with the posterior samples of the chosen parameter, on variance scale for variances parameters and original (the common) scale for fixed effect coefficients

Examples

if (interactive() && requireNamespace("rstan")){
  ex_prior <- makemyprior_example_model()
  res_stan <- inference_stan(ex_prior, iter = 100)
  # Note: For reliable results, increase the number of iterations (e.g., 'iter = 2000')
  extract_posterior_parameter(res_stan, "intercept")
  extract_posterior_parameter(res_stan, "a")
}

Find suitable PC prior parameters

Description

Returns the U value in P(U > sigma) = alpha for a PC prior on standard deviation given an equal-tailed credible interval P(lower < func(x) < upper) = prob where x is a Gaussian variable with zero mean standard deviation sigma. Note that this function uses sampling.

Usage

find_pc_prior_param(
  lower,
  upper,
  alpha = 0.05,
  func = exp,
  N = 10000,
  prob = 0.95
)

Arguments

Value

The U-value to pass to the PC prior. NB! Store result to avoid rerunning this function, as it uses sampling. The function also prints (sampled) quantiles for the U-value that is returned.

Examples

find_pc_prior_param(0.1, 10)

Internal variance parameter order

Description

Returns the internal order of the variance parameters related to each random effect in the model.

Usage

get_parameter_order(prior_obj)

Arguments

Value

Names of the random effects in the model in the order the prior object reads them.

Examples

ex_model <- makemyprior_example_model()
get_parameter_order(ex_model)

Run inference

Description

This function helps you run inference with INLA using a prior object from make_prior. You must have INLA installed to run this. INLA can be installed with: install.packages("INLA", repos = c(getOption("repos"), INLA = "https://inla.r-inla-download.org/R/stable"), dep = TRUE). Also see r-inla.org.

Usage

inference_inla(prior_obj, use_likelihood = TRUE, print_prior = TRUE, ...)

Arguments

Details

Jeffreys' prior is improper. If use_likelihood = FALSE and Jeffreys' prior is used for the total variance, the prior will be changed to a Gaussian(0,1) prior on the log total variance. This means that it does not make sense to look at the variances/standard deviations/precisions, but the variance proportions will be correct. Note that this is only an issue when sampling from the prior (i.e., not using the likelihood).

Value

A named list with a prior object (prior), an inla-object (inla) and some data inla requires (inla_data).

Examples

## Not run: 

vignette("make_prior", package = "makemyprior")

## End(Not run)

ex_prior <- makemyprior_example_model()
if (interactive() && requireNamespace("INLA")){
  posterior <- inference_inla(ex_prior)
  plot(posterior)
}

Run inference

Description

This function helps you run inference with rstan using a prior object from make_prior. Note that you must install Stan: install.packages("rstan"), see mc-stan.org.

Usage

inference_stan(
  prior_obj,
  use_likelihood = TRUE,
  print_prior = TRUE,
  path = NULL,
  ...
)

Arguments

Details

We cannot sample from a Jeffreys' prior since it is improper. If use_likelihood = FALSE and Jeffreys' prior is used for the total variance, the prior will be changed to a Gaussian(0,1) prior on the log total variance. This means that it does not make sense to look at the variances/standard deviations/precisions, but the variance proportions will be correct. Note that this is only an issue when sampling from the prior (i.e., not using the likelihood).

Value

A named list with a prior object (prior), a stan-object (stan) and some data stan requires (stan_data).

Examples

## Not run: 

vignette("make_prior", package = "makemyprior")

## End(Not run)

ex_prior <- makemyprior_example_model()
if (interactive() && requireNamespace("rstan")){
  posterior <- inference_stan(ex_prior, iter = 100)
  # Note: For reliable results, increase the number of iterations (e.g., 'iter = 2000')
  plot(posterior)
}

## Not run: 

posterior <- inference_stan(ex_prior, use_likelihood = TRUE, iter = 1e4, chains = 1, seed = 1)
plot(posterior)

## End(Not run)

Latin square experiment data

Description

Simulated dataset for latin square experiment with 81 observations.

Usage

latin_data

Format

A list with the following variables:

y

Response

lin

Covariate for linear effect of treatment

row

Row indexes

col

Column indexes

treat_iid, treat_rw2

Treatment indexes

Examples

## Not run: 

vignette("latin_square", package = "makemyprior")

## End(Not run)

if (interactive() && requireNamespace("rstan")){

  formula <- y ~ lin + mc(row) + mc(col) + mc(iid) +
    mc(rw2, model = "rw2", constr = TRUE, lin_constr = TRUE)

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

  posterior <- inference_stan(prior, iter = 150, warmup = 50,
                              seed = 1, init = "0", chains = 1)
  # Note: For reliable results, increase the number of iterations

  plot(prior)
  plot_tree_structure(prior)
  plot_posterior_fixed(posterior)
  plot_posterior_stan(posterior, param = "prior", prior = TRUE)
}

## Not run: 

posterior <- inference_stan(prior, iter = 15000, warmup = 5000,
                            seed = 1, init = "0", chains = 1)

plot(prior)
plot_tree_structure(prior)
plot_posterior_fixed(posterior)
plot_posterior_stan(posterior, param = "prior", prior = TRUE)

## End(Not run)

logit

Description

Calculates logit, log(x/(1-x)) = log(x) - log(1-x)

Usage

logit(x)

Arguments

Value

logit value

Examples

logit(0.5)
logit(seq(0, 1, 0.1))

Making a prior object

Description

Make a prior object with all necessary information about the prior and model. The object can either be sent to makemyprior_gui or used directly for inference with Stan (inference_stan) or INLA (inference_inla). eval_joint_prior can be used to evaluate the prior.

Usage

make_prior(
  formula,
  data,
  family = "gaussian",
  prior = list(),
  intercept_prior = c(),
  covariate_prior = list()
)

Arguments

Details

See makemyprior_models for details on available priors and likelihoods.

Value

Prior object.

Examples

## Not run: 

vignette("make_prior", package = "makemyprior")

## End(Not run)

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)

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

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

Returning a simple example prior object

Description

Creating a simple prior object using make_prior. Used in examples of other functions in the package.

Usage

makemyprior_example_model(seed = 1)

Arguments

Details

See the example for what model is made.

Value

An object of class mmp_prior.

Examples

ex_model <- makemyprior_example_model()

## Not run: 
# The function corresponds to the following model call:

set.seed(1)

data <- list(
  a = rep(1:10, each = 10),
  b = rep(1:10, times = 10)
)
data$y <- rnorm(10, 0, 0.4)[data$a] + rnorm(10, 0, 0.6)[data$b] + rnorm(100, 0, 1)

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

prior <- make_prior(formula, data, family = "gaussian",
                    prior = list(tree = "s1 = (a, b); s2 = (s1, eps)",
                                 w = list(s2 = list(prior = "pc0", param = 0.25)),
                                 V = list(s2 = list(prior = "pc", param = c(3, 0.05)))),
                    intercept_prior = c(0, 1000))

## End(Not run)

Graphical prior construction

Description

This functions opens a shiny app where the specified prior can be seen, and changed.

Usage

makemyprior_gui(prior, guide = FALSE, no_pc = FALSE)

Arguments

Value

Returns an object that can be sent to inference_stan or inference_inla. Can also be sent to makemyprior_gui again.

Examples

## Not run: 

vignette("make_prior", package = "makemyprior")


## End(Not run)

if (interactive()){
   ex_prior <- makemyprior_example_model()
   new_prior <- makemyprior_gui(ex_prior)
}

List available priors, latent models and likelihoods

Description

Prints available priors, latent models and likelihoods to use with make_prior.

Usage

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

Arguments

Value

None.

Examples

makemyprior_models("prior", c("pc0", "pc1"))
makemyprior_models("latent")
makemyprior_models("likelihood", "all")

List of available plotting functions

Description

Directs to this help page.

Usage

makemyprior_plotting()

Details

The available plotting functions are:

• plot_marginal_prior

• plot_posterior_fixed

• plot_posterior_precision

• plot_posterior_stan

• plot_posterior_stdev

• plot_posterior_variance

• plot_prior

• plot_several_posterior_stan

• plot_tree_structure

In addition the following functions can be used to extract posterior samples of effects and parameters:

• extract_posterior_effect

• extract_posterior_parameter

eval_pc_prior can be used to evaluate a PC prior for a weight parameter, and eval_joint_prior to evaluate the whole joint prior.

Value

None.

Define latent component

Description

Function for defining a latent component for the HD prior package. All model components must be specified, and if an HD prior is used, this is specified later. See make_prior for more details and examples.

Usage

mc(
  label,
  model = "iid",
  constr = NULL,
  lin_constr = FALSE,
  Cmatrix = NULL,
  graph = NULL,
  ...
)

Arguments

Details

The graph argument is a path to a file describing neighbouring relationship on the following form: First row: number of elements The rest: First number is the index of this element, second is the number of neighbours, the rest is the index numbers of all neighbours for this element. If element 1 and 4 are neighbours, 1 should have 4 in its neighbour list, and 4 should have 1.

Value

For specifying details on this latent component.

Examples

## Not run: 

vignette("make_prior", package = "makemyprior")

## End(Not run)

Neonatal mortality data

Description

Simulated neonatal mortality data with 323 observations.

Usage

neonatal_data

Format

A list with the following variables:

y

Response

Ntrials

Number of trials for each cluster

urban

Covariate indicating if cluster is urban (1) or rural (0)

nu

Cluster effect indexes

v

County effect indexes for iid effect

u

County effect indexes for Besag effect

Examples

## Not run: 

vignette("neonatal_mortality", package = "makemyprior")

## End(Not run)

if (interactive() && requireNamespace("rstan")){

  graph_path <- paste0(path.package("makemyprior"), "/neonatal.graph")

  formula <- y ~ urban + mc(nu) + mc(v) +
    mc(u, model = "besag", graph = graph_path, scale.model = TRUE)

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

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

  posterior <- inference_stan(prior, iter = 150, warmup = 50,
                              seed = 1, init = "0", chains = 1)
  # Note: For reliable results, increase the number of iterations

  plot(prior)
  plot_tree_structure(prior)
  plot_posterior_fixed(posterior)
  plot_posterior_stan(posterior, param = "prior", prior = TRUE)
}

## Not run: 

posterior <- inference_stan(prior, iter = 15000, warmup = 5000,
                            seed = 1, init = "0", chains = 1)

plot(prior)
plot_tree_structure(prior)
plot_posterior_fixed(posterior)
plot_posterior_stan(posterior, param = "prior", prior = TRUE)

## End(Not run)

Plotting

Description

Plotting

Usage

## S3 method for class 'mmp_prior'
plot(x, ...)

## S3 method for class 'mmp_inla'
plot(x, ...)

## S3 method for class 'mmp_stan'
plot(x, ...)

Arguments

Details

See plot_prior (objects of class mmp_prior), plot_posterior_stan (objects of class mmp_stan), and plot_posterior_variance (objects of class mmp_inla),

Value

None.

Examples

pri <- makemyprior_example_model()
plot(pri)

if (interactive() && requireNamespace("rstan")){
  res_stan <- inference_stan(ex_prior, iter = 100)
  # Note: For reliable results, increase the number of iterations (e.g., 'iter = 2000')
  plot(res_stan)
}

if (interactive() && requireNamespace("INLA")){
  res_inla <- inference_inla(pri)
 plot(res_inla)
}

Plotting prior for a single parameter (weight or variance (not standard deviation))

Description

Following the parameterization of the prior.

Usage

plot_marginal_prior(x, obj, param, sd = FALSE)

Arguments

Value

A ggplot with the posterior distribution. See also makemyprior_plotting.

Examples

ex_prior <- makemyprior_example_model()
plot_marginal_prior(seq(0, 1, 0.001), ex_prior, "w[a/a_b]")
plot_marginal_prior(seq(0, 1, 0.001), ex_prior, "w[eps/eps_a_b]")
plot_marginal_prior(seq(0, 5, 0.01), ex_prior, "V[eps_a_b]")

Plotting posterior distributions

Description

Function for plotting the posterior distributions of the coefficients of the fixed effects

Usage

plot_posterior_fixed(obj)

Arguments

Value

A ggplot with the posterior distributions. See also makemyprior_plotting.

Examples

if (interactive() && requireNamespace("rstan")){
  ex_prior <- makemyprior_example_model()
  res_stan <- inference_stan(ex_prior, iter = 100)
  # Note: For reliable results, increase the number of iterations (e.g., 'iter = 2000')
  plot_posterior_fixed(res_stan)
}

Plotting posterior distributions

Description

Function for plotting the posterior distributions of the random effect variances on the scale of the tree parameterization.

Usage

plot_posterior_stan(
  obj,
  param = c("prior", "variance", "stdev", "precision"),
  prior = FALSE
)

Arguments

Value

A ggplot with the posterior distributions. See also makemyprior_plotting.

Examples

if (interactive() && requireNamespace("rstan")){
  ex_prior <- makemyprior_example_model()
  res_stan <- inference_stan(ex_prior, iter = 100)
  # Note: For reliable results, increase the number of iterations (e.g., 'iter = 2000')
  plot_posterior_stan(res_stan)
}

Plotting posterior variances, standard deviations or precisions

Description

Plotting posterior variances, standard deviations or precisions

Usage

plot_posterior_variance(obj)

plot_posterior_stdev(obj)

plot_posterior_precision(obj)

Arguments

Value

A ggplot object with the plot See also makemyprior_plotting.

Examples

if (interactive() && requireNamespace("rstan")){
  ex_prior <- makemyprior_example_model()
  res_stan <- inference_stan(ex_prior, iter = 100)
  # Note: For reliable results, increase the number of iterations (e.g., 'iter = 2000')
  plot_posterior_variance(res_stan)
}

if (interactive() && requireNamespace("INLA")){
  ex_prior <- makemyprior_example_model()
  res_inla <- inference_inla(ex_prior)
  plot_posterior_variance(res_inla)
}

Plotting prior distributions

Description

Function for plotting the prior distributions of the random effects on the scale of the parameters chosen

Usage

plot_prior(obj)

Arguments

Value

A ggplot with the prior distributions. See also makemyprior_plotting.

Examples

ex_prior <- makemyprior_example_model()
plot_prior(ex_prior)

Plotting several posterior distributions

Description

Function for plotting the posterior distributions of the random effect variances on the scale of the tree parameterization.

Usage

plot_several_posterior_stan(
  objs,
  param = c("prior", "variance", "stdev", "precision", "logvariance")
)

Arguments

Details

We cannot sample from a Jeffreys' prior since it is improper. If Jeffreys' prior is used for the total variance, the prior will be changed to a Gaussian(0,1) prior on the log total variance. This means that it does not make sense to look at the variances/standard deviations/precisions, but the variance proportions will be correct. See also makemyprior_plotting.

Value

A ggplot with the posterior distributions.

Examples

if (interactive() && requireNamespace("rstan")){
  ex_prior1 <- makemyprior_example_model(seed = 1)
  ex_prior2 <- makemyprior_example_model(seed = 2)
  # Note: For reliable results, increase the number of iterations (e.g., 'iter = 2000')
  res_stan1 <- inference_stan(ex_prior1, iter = 100)
  res_stan2 <- inference_stan(ex_prior2, iter = 100)
  plot_several_posterior_stan(list(One = res_stan1, Two = res_stan2))
}

Plotting the prior tree structure graph

Description

Can only be used for visualization in R.

Usage

plot_tree_structure(obj, nodenames)

Arguments

Value

A visNetwork with the tree graph See also makemyprior_plotting.

Examples

ex_prior <- makemyprior_example_model()
plot_tree_structure(ex_prior)

Print

Description

Print

Usage

## S3 method for class 'mmp_prior'
print(x, ...)

## S3 method for class 'mmp_inla'
print(x, ...)

## S3 method for class 'mmp_stan'
print(x, ...)

Arguments

Value

Returns input object invisible.

Examples

pri <- makemyprior_example_model()
pri # or print(pri)

if (interactive() && requireNamespace("rstan")){
  res_stan <- inference_stan(ex_prior, iter = 100)
  # Note: For reliable results, increase the number of iterations (e.g., 'iter = 2000')
  res_stan # or print(res_stan)
}

if (interactive() && requireNamespace("INLA")){
  res_inla <- inference_inla(pri)
  res_inla # or print(res_inla)
}

Scaling precision matrix

Description

Scaling a precision matrix so the corresponding covariance matrix has typical variance (geometric mean) equal to 1.

Usage

scale_precmat(Q)

Arguments

Value

Precision matrix which is now scaled to have typical variance 1.

Examples

scale_precmat(diag(10))

Short summary

Description

Short summary

Usage

## S3 method for class 'mmp_prior'
summary(object, ...)

## S3 method for class 'mmp_inla'
summary(object, ...)

## S3 method for class 'mmp_stan'
summary(object, ...)

Arguments

Value

Returns summary invisible.

Examples

pri <- makemyprior_example_model()
summary(pri)

if (interactive() && requireNamespace("rstan")){
  res_stan <- inference_stan(ex_prior, iter = 100)
  # Note: For reliable results, increase the number of iterations (e.g., 'iter = 2000')
  summary(res_stan)
}

if (interactive() && requireNamespace("INLA")){
  res_inla <- inference_inla(pri)
  summary(res_inla)
}

Compute the typical variance

Description

Computing the typical variance (geometric mean) of a matrix.

Usage

typical_variance(mat)

Arguments

Value

Typical variance.

Examples

typical_variance(diag(10))

Genomic wheat breeding model data

Description

Simulated wheat yield data with 100 observations.

Usage

wheat_data

Format

A list with the following variables

y

Response

a, b, x

Indexes for the additive, dominance and epistasis genetic effects, respectively

Q_a, Q_d, Q_x

Precision matrices for the genetic effects

Examples

## Not run: 

vignette("wheat_breeding", package = "makemyprior")

## End(Not run)

if (interactive() && requireNamespace("rstan")){

  wheat_data_scaled <- wheat_data
  wheat_data_scaled$Q_a <- scale_precmat(wheat_data$Q_a)
  wheat_data_scaled$Q_d <- scale_precmat(wheat_data$Q_d)
  wheat_data_scaled$Q_x <- scale_precmat(wheat_data$Q_x)

  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)

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

  posterior <- inference_stan(prior, iter = 150, warmup = 50,
                              chains = 1, seed = 1)
  # Note: For reliable results, increase the number of iterations

  plot(prior)
  plot_tree_structure(prior)
  plot_posterior_fixed(posterior)
  plot_posterior_stan(posterior, param = "prior", prior = TRUE)

}

## Not run: 

posterior <- inference_stan(prior, iter = 150, warmup = 50,
                            chains = 1, seed = 1)

plot(prior)
plot_tree_structure(prior)
plot_posterior_fixed(posterior)
plot_posterior_stan(posterior, param = "prior", prior = TRUE)

## End(Not run)