Prior Predictive Simulation

Simulate numerator and denominator responses from the prior predictive distribution: parameters are drawn from priors, the two linear predictors are assembled on the design formula and data imply, and a response is drawn from family at each draw. Comparing the resulting ratios against what is plausible for the system under study is the check on whether the priors are reasonable before any data are fitted.

Random-effect terms contribute at a standard deviation drawn from priors$sigma, one per intercept and one per slope column. Correlation between slopes within a term is not drawn, so a correlated term simulates at its marginal scales.

Usage

prior_predict(formula, family, data, priors = NULL, n = 100)

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

Arguments

formula

A model formula, using the combined num | denom ~ x syntax or a single response with formula_denom supplied through ratiod_formula().

family

A ratiod_family object.

data

Data frame supplying the predictors, the grouping factors and, for the binomial families, the trial counts. Response columns need not be present; where one is missing a placeholder is used, since nothing about the prior predictive distribution depends on the observed response.

priors

A ratiod_priors object. Defaults to ratiod_priors().

n

Number of prior predictive draws.

x

A ratiod_priorpred object

...

Ignored

Value

An object of class ratiod_priorpred: a list with the n x N matrices y_num, y_denom, ratio, mu_num and mu_denom, the drawn coefficients beta_num and beta_denom, the drawn random-effect scales sigma, the drawn dispersions phi_num and phi_denom, and the family and data the draws were taken on.

x, invisibly

See also

sim_ratiod() for simulating from a known parameterisation, ratiod_priors() for the prior specification

Examples

set.seed(1)
df <- data.frame(
  count = rep(0L, 40),
  total = rep(0L, 40),
  x = rnorm(40),
  site = factor(rep(1:8, each = 5))
)

pp <- prior_predict(
  count | total ~ x + (1 | site),
  family = ratiod_negbin_negbin(),
  data = df,
  n = 50
)
dim(pp$ratio)
#> [1] 50 40
summary(as.vector(pp$ratio))
#>       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
#>  0.000e+00  0.000e+00  3.000e+00 7.360e+164  3.138e+07 1.395e+168 

# A tighter coefficient prior narrows the predicted ratios
pp_tight <- prior_predict(
  count | total ~ x + (1 | site),
  family = ratiod_negbin_negbin(),
  data = df,
  priors = ratiod_priors(beta = prior_normal(0, 0.5)),
  n = 50
)