Simulate responses from a formula and parameters

Builds the design from formula and data exactly as frm() would, sets the parameters from newparams (and/or draws them from prior), and simulates responses. Random effects are redrawn from their covariance for every simulation unless newparams$b supplies them.

Usage

frm_simulate(
  formula,
  data,
  family = NULL,
  newparams = NULL,
  prior = NULL,
  nsim = 1,
  seed = NULL,
  data2 = list()
)

Arguments

formula: A bf() formula (with a family attached) or a plain formula plus family.
data: A data frame of model data, including a dummy response column.
family: Family, when formula does not carry one.
newparams: Named list of parameters, in either spelling (see Details). par_template() discovers the internal spelling for a formula and data. Optional when prior pins everything.
prior: A set_prior() specification to draw parameters from, once per simulation, or a brmsprior object brms built, which is translated row by row. The argument takes brms's spelling, prior; the priors of releases before 0.43 is gone rather than aliased.
nsim, seed: As in simulate().
data2: Structural objects, as in frm().

Value

A data frame with nsim columns of simulated responses, carrying the drawn parameters in attr(., "pars") when prior is used.

Details

data must contain a response column with values that are valid for the family (any dummy values do; they only anchor the design).

Draws come back in the response's own type, exactly as simulate()'s do: an ordered factor for an ordinal family, an unordered one for a categorical family, and a matrix column for a matrix response (multinomial(), frmtmb.latent's mixture_mvn()).

The structured families draw here through the same implementation simulate() uses (see its Structured draws section): mixture(groups = ~g) takes one class per group, a multivariate normal mixture uses its class covariances, a residual correlation term (ar(), ma(), ...) contributes one correlated residual per group, and a family from an extension package draws through whatever its structure declares. The de novo frame carries those structures, so nothing is lost relative to a fit; ar() and friends need their thetaac entry in the internal newparams spelling, since a correlation parameter has no natural-scale name here.

Two spellings for `newparams`

Natural scale (recommended): brms's parameter names, the strings variables() returns, one number each.

coefficients as b_Intercept, b_x, b_sigma_Intercept, on the LINK scale;
sigma (and any other distributional parameter nobody wrote a formula for: shape, phi, zi, ...) on the RESPONSE scale, so sigma = 0.7 is a residual SD of 0.7;
a mixture's weights with no theta formula as brms's simplex, all of theta1 ... thetaK together, each in (0, 1) and summing to one: theta1 = 0.9, theta2 = 0.1;
sd_<group>__<coef> for group-level standard deviations, cor_<group>__<c1>__<c2> for their correlations, and sds_<label>_1 for a smooth's smoothing SD (sds_sx_1 for s(x)). Unset correlations are 0.

A bare name (x, Intercept, sigma_Intercept) is refused with its brms spelling. Priors are different: set_prior(coef = "x") takes the bare coefficient, as in brms.

Internal scale: named after the parameter components - beta, betad, theta, and optionally b - each a full-length vector, on the internal parameterization (theta holds log SDs and Cholesky correlation parameters, betad holds dispersion dpars on their link scale). par_template() returns that layout for a formula and data without fitting anything, filled with the default values and carrying the name of every entry; edit it and pass it back as newparams.

The two spellings cannot be mixed: newparams is read as internal when every name is a par_template component (or b), and as natural otherwise. The natural spelling covers us, diag, homdiag, smooth, gr_cov and gr_prec random-effect structures; other structures (ar1, cs, toep, GPs, ...) have no inverse map from SDs and correlations and need the internal spelling.

Prior-predictive simulation

With prior, each of the nsim simulations draws a fresh parameter vector from the set_prior() specification and simulates from it - the analog of brms's sample_prior = "only" followed by posterior_predict(). Draws are taken on the scale each class defines: class "b"/"Intercept" on the link scale, class "sd" on the natural SD scale (mapped to theta afterwards), class "theta" on the internal scale. lb/ub truncate by rejection. The drawn values come back as attr(result, "pars"), one row per simulation, so a prior-predictive check can relate parameters to outcomes.

That table reports each parameter on the scale newparams names it on, which for b_Intercept is the intercept at ZERO. A class = "Intercept" prior is a density on the intercept at the MEAN of the predictors (see the Where an intercept prior lands section of set_prior()), so on a design with uncentered predictors the number drawn and the number reported differ by colMeans(X) times the slopes. Reporting the written value is what lets a pars row be handed straight back as newparams. A prior on a distributional parameter's own brms class is likewise drawn on that parameter and reported there, not on its link scale.

The column name says which scale a dpar's draw is on, as brms's names do: sigma holds sigma, drawn under set_prior(class = "sigma") on a model that gives sigma no formula, and b_sigma_Intercept holds log sigma, drawn under set_prior(class = "Intercept", dpar = "sigma") on a model that does.

Parameters without a prior keep their newparams value. Whenever prior are used, or newparams uses the natural spelling, every fixed coefficient and every random-effect SD must be pinned by one or the other; an unpinned parameter is an error rather than a silent zero effect or unit SD.

Examples

# power analysis: simulate from a design with chosen parameters
dd <- data.frame(x = rnorm(60), g = factor(rep(1:6, 10)), y = 0)
sims <- frm_simulate(bf(y ~ x + (1 | g)) + gaussian(), dd,
                     newparams = list(b_Intercept = 1, b_x = 0.5,
                                      sigma = 0.7,
                                      sd_g__Intercept = 0.5),
                     nsim = 3, seed = 1)
head(sims)
#>       sim_1      sim_2      sim_3
#> 1 0.8647286 0.36178434  0.6460370
#> 2 1.9956516 1.46591169  0.5879778
#> 3 1.3777358 0.59158649  0.2057305
#> 4 1.9654916 2.67194124  1.4370064
#> 5 2.3704051 1.07485483 -0.1685644
#> 6 0.2364781 0.01962257 -0.2147892
# the same thing on the internal scale
par_template(bf(y ~ x + (1 | g)) + gaussian(), dd)   # the layout
#> <frmtmb parameter template> starting values
#> $beta  (2)
#> (Intercept)           x 
#>           0           0 
#> $betad  (1)
#> sigma_(Intercept) 
#>                 0 
#> $b  (6)
#> b_1 b_2 b_3 b_4 b_5 b_6 
#>   0   0   0   0   0   0 
#> $theta  (1)
#> theta_1 
#>       0 
#> Edit and pass back as frm(start = ) or frm_simulate(newparams = ). Parentheses are optional in names you supply.
sims2 <- frm_simulate(bf(y ~ x + (1 | g)) + gaussian(), dd,
                      newparams = list(beta = c(1, 0.5),
                                       betad = log(0.7),
                                       theta = log(0.5)),
                      nsim = 3, seed = 1)
# prior-predictive draws
pp <- frm_simulate(bf(y ~ x + (1 | g)) + gaussian(), dd,
                   prior = set_prior("normal(0, 1)", class = "b") +
                     set_prior("normal(0, 2)", class = "Intercept") +
                     set_prior("exponential(1)", class = "sd") +
                     set_prior("exponential(1)", class = "sigma"),
                   nsim = 4, seed = 1)
head(attr(pp, "pars"))
#>          b_x b_Intercept sd_g__Intercept     sigma
#> 1 -0.6264538   3.1939242       0.1457067 0.1397953
#> 2  0.4755095  -2.5098192       1.5652413 0.8145358
#> 3  1.1765833  -0.9333763       3.3072809 1.4953203
#> 4 -0.1643758  -0.7996112       0.3260219 1.5316271