Simulate response datasets for predictive checks

Generates new response datasets for the same observations and predictors. Reuse them with check_partner_dependence() to check whether a fitted model reproduces features of the observed data. For a complete example see check_partner_dependence().

Usage

simulate_dyad_responses(model, nsim = 1000, seed = NULL)

Arguments

model

A fitted glmmTMB model.

nsim

Number of complete response datasets to simulate. Default: 1000.

seed

Optional seed for reproducible simulations, interpreted as an integer (see set.seed()). When supplied, the random-number state is restored after the function returns, or when it stops after an error.

Value

A dyadMLM_response_simulations object for use with check_partner_dependence().

The result keeps all components in fitted-row order (after missing-data exclusions):

• simulated_responses: a matrix with nsim rows and one column per fitted observation. Each row is a complete simulated dataset.

• observed_response and predicted_response: numeric vectors with one value per fitted observation.

• model_frame: the fitted model frame (variables in the model formulas only), in the same row order.

The dyadMLM attribute records the model and simulation settings, including the seed.

Supported models

Supports unweighted glmmTMB models with the following families:

• gaussian()

• poisson(), glmmTMB::compois(), glmmTMB::genpois(), glmmTMB::bell()

• glmmTMB::nbinom1(), glmmTMB::nbinom2(), glmmTMB::nbinom12()

• glmmTMB::truncated_poisson(), glmmTMB::truncated_nbinom1(), glmmTMB::truncated_nbinom2(), glmmTMB::truncated_compois(), glmmTMB::truncated_genpois()

• glmmTMB::tweedie()

• Gamma(), glmmTMB::ziGamma()

• glmmTMB::beta_family()

• glmmTMB::lognormal()

• glmmTMB::skewnormal()

• glmmTMB::t_family() with more than two degrees of freedom

• glmmTMB::ordinal() (currently only available in the development version of glmmTMB)

Smooth terms (s()) are not supported: simulations would draw new curves. Fixed-effect splines such as splines::ns() are supported.

Zero-inflated and hurdle versions are supported where available. Checks describe the combined response, including zeros, rather than each model component separately. Good agreement does not establish that the zero and response components each fit well.

Ordinal checks use category scores 1, 2, ..., K in their fitted order, matching glmmTMB's predictions. The plots compare variation and partner correlation in these scores. The scores do not measure distances on an underlying continuous scale.

The model's fitted link is used for prediction and simulation. Predictions and simulated responses must be finite.

check_partner_dependence() currently requires cross-sectional dyads.

Technical details

Each simulation draws new random effects and then new responses from the fitted model. Random effects within each block are drawn together using their fitted variances and correlations. This also applies to random effects in the zero-inflation and dispersion models, if present.

Fitted parameters and predictors, including any lagged responses, stay fixed. The model is not refitted, and uncertainty in parameter estimates is not included. This is a plug-in predictive reference. If dyads are the only grouping factor, the simulations represent hypothetical new dyads under the same study design.

predicted_response contains predicted mean responses with random effects in the conditional and zero-inflation models set to zero. For zero-inflated and hurdle models, this is the conditional response mean multiplied by one minus the zero-component probability (Brooks et al., 2017, Appendix A; doi:10.32614/RJ-2017-066 ).

By default, later checks subtract these same predictions from observed and simulated responses. Both random effects and observation-level noise still contribute to response variance. With nonlinear links, setting random effects to zero generally differs from averaging predictions over them.