Package 'dyadMLM'

Description

A dyadic model should reproduce how much responses vary and how strongly partners' responses are related. This check computes variances and correlations of simulated responses based on the model and compares them to the observed response variances and correlations from the data. This helps identify mismatches in the model's assumptions. The check currently supports cross-sectional dyads only (one response per partner).

Usage

check_partner_dependence(
  simulations,
  dyad,
  role = NULL,
  plot = TRUE,
  response = c("model-centred", "raw"),
  ask = NULL,
  panels = TRUE,
  data = NULL
)

Arguments

Value

The comparison plots (shown by default) are the main output. The function invisibly returns a dyadMLM_partner_check object containing the compositions table with pair counts and a statistics tibble for each composition. Compositions with fewer than three complete dyads are skipped and have NULL statistics. Each tibble has one observed row followed by one row per simulation, identified by dataset. summary compares each observed statistic with the middle 95% of its simulations (lower, upper) and marks those outside it. The object includes omission counts and settings, and can be saved and plotted later.

Reading the plots

Histograms show simulated summaries. Red lines mark observed values. Dashed lines enclose the middle 95% of simulations (no formal confidence intervals). The heading shows the dyad composition, its number of usable dyads, and the total across all compositions. The top row shows member SDs and partner correlation. The bottom row shows the same information using dyad averages and partner differences.

An observed value far from most simulated values may indicate that the model does not reproduce that feature of the data well. Agreement does not establish that omitted dependence is negligible, especially with few dyads. The partner-dependence study illustrates how sample size affects detection when residual partner correlation is omitted.

Checking each composition can reveal differences hidden by pooling. Flags (observed values outside the middle 95%, marked * when printed) can occur by chance, especially when checking several summaries. They invite investigation, not formal rejection of the model. The covariance-pooling study illustrates detection and false alarms when checking each composition; 10 to 29% of correctly pooled models had at least one flag.

The top row compares:

• Response SDs: one for each role, or one common SD for exchangeable members.

• Partner correlation: how strongly partners' responses are related.

The bottom row shows:

• Dyad-average SD: how much dyads differ in their average response.

• Half-difference SD or RMS: each partner difference is divided by two. With distinct roles, the SD shows how much signed differences vary across dyads. For exchangeable members, the RMS (root mean square) shows their typical size, regardless of partner order.

• Dyad-average/role-difference correlation: shown only for distinct roles, because it depends on how partners are ordered. Positive values indicate greater variance for the first named role. Negative values indicate greater variance for the second.

These checks are particularly useful when a less restricted model does not converge and so cannot be used for model comparison. They show how well the simpler model reproduces the observed variances and partner correlations.

A flexible covariance model will usually reproduce features it estimated from the same data. Agreement alone therefore does not establish good fit. For example, pooled summaries (role = NULL) of a Gaussian model with an unconstrained exchangeable covariance, as below, agree almost by construction. Use a suitable model comparison to formally test a specific covariance restriction when both models can be fitted (see compare_nested_models()).

Rows with missing IDs or roles and incomplete dyads are omitted with a warning.

Technical details

After any centring, paired responses a and b are used to compute dyad averages M = (a + b) / 2 and half-differences D = (a - b) / 2. Roles follow factor levels or sorted values.

For exchangeable members, common member variance is var(M) + mean(D^2) and partner covariance is var(M) - mean(D^2). Partner correlation is covariance divided by variance. Half-difference RMS is sqrt(mean(D^2)).

The variance and covariance calculations for exchangeable dyads follow Woody and Sadler (2005).

References

Woody, E., & Sadler, P. (2005). Structural equation models for interchangeable dyads: Being the same makes a difference. Psychological Methods, 10(2), 139-158. doi:10.1037/1082-989X.10.2.139.

Gelman, A., Meng, X.-L., & Stern, H. S. (1996). Posterior predictive assessment of model fitness via realized discrepancies. Statistica Sinica, 6, 733-807.

Examples

# Data contains three compositions: female-female, female-male, and male-male.
example_data <- prepare_dyad_data(
  dyads_cross,
  dyad = coupleID,
  member = personID,
  model_types = "none",
  seed = 123
)

# This model pools all compositions and treats every dyad as exchangeable.
model <- glmmTMB::glmmTMB(
  closeness ~ 1 +
    us(1 | coupleID) +
    us(0 + .member_contrast_arbitrary | coupleID),
  dispformula = ~ 0,
  data = example_data
)

# Fewer simulations for a quick example (the default is 1000).
simulations <- simulate_dyad_responses(
  model,
  nsim = 100,
  seed = 123
)

check_partner_dependence(
  simulations,
  dyad = coupleID,
  role = gender,
  # Supply the fitting data because gender is not in the model formula.
  data = example_data
)

# Optionally, suppress plot, store object, and plot later.
check <- check_partner_dependence(
  simulations,
  dyad = coupleID,
  role = gender,
  data = example_data,
  plot = FALSE
)

# Check how well the model reproduces variances and partner correlations
# within each dyad composition.
plot(check, ask = FALSE, panels = TRUE)
# Composition checks flag differences that the pooled check misses.
print(check)

Description

Performs a likelihood-ratio test for two nested glmmTMB models. The models may use ordinary data frames or dyadMLM_data objects, and their calls do not need to refer to the same R object. Models may be supplied in either order. The model with fewer estimated parameters is shown first in the result.

Usage

compare_nested_models(model1, model2, alpha = 0.05)

Arguments

Details

Both model calls must use named data-frame objects that remain available when the models are compared. The checks assume these objects have not been modified since fitting. All ordinary data columns must be identical, including their types and attributes. For dyadMLM_data, package-generated columns may differ, but the original columns must be identical. Ordinary and prepared data may be compared with each other. Dyad metadata are checked when both models use dyadMLM_data. The function also checks fitted rows, outcomes, weights and offsets, model family and link, maximum-likelihood estimation, and model convergence. Each model must use the same untransformed response column.

These checks establish that the models use equivalent observations. They cannot establish that one model is mathematically nested within the other. The caller remains responsible for supplying genuinely nested models. The usual chi-squared reference distribution may also be inappropriate when tested variance parameters are on the boundary. alpha affects only the printed conclusion, not the test results.

Value

An anova-style data frame containing model degrees of freedom, information criteria, log-likelihoods, the likelihood-ratio statistic, and its chi-squared p-value. When printed, a short conclusion interprets the test at the selected significance level.

Examples

if (requireNamespace("glmmTMB", quietly = TRUE)) {
  restricted_data <- prepare_dyad_data(
    dyads_cross,
    dyad = coupleID,
    member = personID,
    role = gender,
    # All three observed compositions in `dyads_cross` are detected and retained
    # by default. This example focuses on `female-male` dyads, so we restrict the
    # analysis here.
    keep_compositions = "female-male",
    include_arbitrary_member_contrast = TRUE,
    seed = 123
  )
  full_data <- restricted_data

  # Compare an exchangeable model with a distinguishable model.
  restricted_model <- glmmTMB::glmmTMB(
    closeness ~ 1 +
      us(1 | coupleID) +
      us(0 + .member_contrast_arbitrary | coupleID),
    dispformula = ~ 0,
    data = restricted_data
  )

  full_model <- glmmTMB::glmmTMB(
    closeness ~ 0 + gender +
      us(0 + gender | coupleID),
    dispformula = ~ 0,
    data = full_data
  )

  # Test equal means and residual variances jointly.
  compare_nested_models(restricted_model, full_model)
}

Description

A simulated long-format dataset containing distinguishable female-male dyads and exchangeable female-female and male-male dyads. Each dyad has two members and each member has one row. closeness and provided_support are the member's averages across the 14 observations in dyads_ild.

Usage

dyads_cross

Format

A data frame with 720 rows and 5 variables:

personID

Unique person identifier.

coupleID

Dyad identifier.

gender

Gender role, with levels female and male.

closeness

Mean simulated Gaussian closeness score across 14 days.

provided_support

Mean simulated provided-support score across 14 days.

Details

Prepare with dyad = coupleID, member = personID, and role = gender. These data contain three compositions. Use keep_compositions in prepare_dyad_data() when an analysis should retain only selected compositions; omit it when all supplied compositions should remain.

Source

Simulated for dyadMLM; design adapted from https://github.com/Pascal-Kueng/05DyadicDataAnalysis. See Küng, P. M. (2026). Distinguishable and Exchangeable Dyads: Bayesian Multilevel Modelling. Zenodo. doi:10.5281/zenodo.17400655.

Description

A simulated long-format dataset containing distinguishable female-male dyads and exchangeable female-female and male-male dyads. Each dyad has two members observed on 14 diary days.

Usage

dyads_ild

Format

A data frame with 10,080 rows and 6 variables:

personID

Unique person identifier.

coupleID

Dyad identifier.

diaryday

Measurement day, from 0 through 13.

gender

Gender role, with levels female and male.

closeness

Simulated Gaussian closeness outcome.

provided_support

Simulated provided-support score.

Details

The Gaussian outcome includes independent member-specific stationary AR(1) residual components. Their marginal standard deviation is 0.60 for both genders, with lag-1 correlations within these components of 0.55 for female members and 0.50 for male members. Separate shared and difference residuals induce same-occasion partner covariance.

Prepare with dyad = coupleID, member = personID, role = gender, and time = diaryday. These data contain three compositions. Use keep_compositions in prepare_dyad_data() when an analysis should retain only selected compositions; omit it when all supplied compositions should remain.

Source

Simulated for dyadMLM; design adapted from https://github.com/Pascal-Kueng/05DyadicDataAnalysis. See Küng, P. M. (2026). Distinguishable and Exchangeable Dyads: Bayesian Multilevel Modelling. Zenodo. doi:10.5281/zenodo.17400655.

Description

A simulated long-format dataset containing the same dyads, members, and dyad compositions as dyads_cross. Each member has one row. stress is the member's average across the 14 observations in dyads_nbinom_ild. The count outcome is a separate negative-binomial draw using the NB2 variance function; it is not an average or sum of the daily counts.

Usage

dyads_nbinom_cross

Format

A data frame with 720 rows and 5 variables:

personID

Unique person identifier.

coupleID

Dyad identifier.

gender

Gender role, with levels female and male.

conflict_count

Simulated number of conflictual interactions reported by a member.

stress

Mean simulated stress score across 14 days.

Details

Prepare with dyad = coupleID, member = personID, and role = gender. These data contain three compositions. Use keep_compositions in prepare_dyad_data() when an analysis should retain only selected compositions; omit it when all supplied compositions should remain.

Source

Simulated for dyadMLM; design adapted from https://github.com/Pascal-Kueng/05DyadicDataAnalysis. See Küng, P. M. (2026). Distinguishable and Exchangeable Dyads: Bayesian Multilevel Modelling. Zenodo. doi:10.5281/zenodo.17400655.

Description

A simulated long-format dataset containing the same dyads, members, diary days, and dyad compositions as dyads_ild. The outcome is a count generated from a negative-binomial distribution using the NB2 variance function.

Usage

dyads_nbinom_ild

Format

A data frame with 10,080 rows and 6 variables:

personID

Unique person identifier.

coupleID

Dyad identifier.

diaryday

Measurement day, from 0 through 13.

gender

Gender role, with levels female and male.

conflict_count

Simulated number of conflictual interactions reported by a member on that day.

stress

Simulated stress score.

Details

Prepare with dyad = coupleID, member = personID, role = gender, and time = diaryday. These data contain three compositions. Use keep_compositions in prepare_dyad_data() when an analysis should retain only selected compositions; omit it when all supplied compositions should remain.

Source

Simulated for dyadMLM; design adapted from https://github.com/Pascal-Kueng/05DyadicDataAnalysis. See Küng, P. M. (2026). Distinguishable and Exchangeable Dyads: Bayesian Multilevel Modelling. Zenodo. doi:10.5281/zenodo.17400655.

Description

Validates dyadic data, records the structural variables, and adds metadata and model-ready columns for dyadic multilevel model parameterizations.

Usage

prepare_dyad_data(
  data,
  dyad,
  member,
  role = NULL,
  time = NULL,
  predictors = NULL,
  lag1_predictors = NULL,
  model_types = "apim",
  dsm_role_order = NULL,
  temporal_decomposition = c("auto", "2l", "none"),
  set_exchangeable_compositions = NULL,
  keep_compositions = NULL,
  pool_compositions = NULL,
  incomplete_dyads = c("error", "drop"),
  missing_role = c("error", "drop"),
  seed = NULL,
  short_colnames = TRUE,
  include_arbitrary_member_contrast = FALSE,
  add_apim_gmc_predictors = FALSE
)

Arguments

Details

Data must be in long format. Cross-sectional dyadic data may contain at most one row per member within dyad. Intensive longitudinal dyadic data may contain at most one row per member and observed measurement occasion within dyad. Measured variables may contain missing values. Numeric dyad, member, time, and selected predictor columns must not contain infinite values. Structural completeness is assessed across all rows. incomplete_dyads controls dyads with fewer than two members; dyads with more than two members always cause an error. When role is supplied, stable member roles are resolved across repeated rows before missing_role is applied.

Dyad composition labels are canonical: role labels are sorted alphabetically before being combined, so labels do not depend on row or member order.

Value

The original data as a tibble with class dyadMLM_data, .composition and .composition_role factor columns, ⁠.is_*⁠ numeric indicator columns, and numeric ⁠.member_contrast_*⁠ columns coded -1 and 1 for the two members of each matching composition and 0 otherwise. With one final composition, their default names omit the composition label. The dyadMLM attribute contains structural metadata, dyad_compositions, and predictor metadata such as temporal_decompositions, lag1_predictors, apim_predictors, and dim_predictors, as well as dsm_predictors and dsm_role_order when applicable. The generated_columns table records each package-generated column retained in the returned data. Modifying prepared data does not automatically regenerate package-generated columns.

Examples

data <- data.frame(
  dyad_id = c(1, 1, 2, 2, 3, 3),
  person_id = c(1, 2, 3, 4, 5, 6),
  role = c("female", "male", "female", "female", "male", "male"),
  x = c(4, 7, 5, 6, 3, 8)
)

prepared <- prepare_dyad_data(
  data,
  dyad = dyad_id,
  member = person_id,
  role = role,
  predictors = x,
  model_types = "apim"
)

print(prepared)

pooled <- prepare_dyad_data(
  data,
  dyad = dyad_id,
  member = person_id,
  role = role,
  predictors = x,
  model_types = "apim",
  set_exchangeable_compositions = "female-male",
  pool_compositions = list(
    romantic_couples = c("female-female", "male-male", "female-male")
  )
)

print(pooled)

ild_data <- data.frame(
  dyad_id = rep(c(1, 2), each = 4),
  person_id = rep(c(1, 2), times = 4),
  time = rep(c(1, 1, 2, 2), times = 2),
  x = c(4, 7, 5, 8, 3, 6, 4, 7)
)

ild_prepared <- prepare_dyad_data(
  ild_data,
  dyad = dyad_id,
  member = person_id,
  time = time,
  predictors = x,
  lag1_predictors = x,
  model_types = "apim",
  seed = 123
)

print(ild_prepared)

Description

Back-transforms covariance matrices from paired shared and member-difference random-effect blocks to the covariance structure of two exchangeable members. The result is on the fitted random effects' linear-predictor scale. In non-Gaussian models, it therefore describes a Gaussian latent covariance, not response-scale residual covariance. For the model specification, derivation, and interpretation, see the exchangeable APIM vignette.

Usage

recover_exchangeable_covariance(
  model,
  block_pairings = NULL,
  posterior = c("mean", "median", "draws")
)

Arguments

Details

Automatic matching recognizes the compact .member_contrast_arbitrary and composition-qualified ⁠.member_contrast_*_arbitrary⁠ names. It first looks for the corresponding ⁠.is_*⁠ shared block. It requires the two blocks to use the same grouping factor and the same underlying terms. Most models fitted with dyadMLM-generated columns therefore need only:

result <- dyadMLM::recover_exchangeable_covariance(model)
print(result)

Supply block_pairings when automatic matching is ambiguous or when a model uses custom indicators or deliberately omitted blocks or terms. To specify one pair with a custom difference indicator:

result <- dyadMLM::recover_exchangeable_covariance(
  model,
  block_pairings = list(
    shared_block = "(1 + time | coupleID)",
    difference_block = "(0 + my_diff + I(my_diff * time) | coupleID)",
    difference_indicator = "my_diff"
  )
)

For multiple covariance levels, wrap the pairings in an outer list. For example, in a Gaussian glmmTMB model fitted with dispformula = ~ 0, this call recovers both a stable dyad-level covariance with an omitted difference time slope and the same-occasion partner residual covariance:

result <- dyadMLM::recover_exchangeable_covariance(
  model,
  block_pairings = list(
    dyad = list(
      shared_block = "(1 + diaryday | coupleID)",
      difference_block = "(0 + .member_contrast_arbitrary | coupleID)",
      difference_indicator =
        ".member_contrast_arbitrary"
    ),
    same_occasion = list(
      shared_block = "(1 | coupleID:diaryday)",
      difference_block = "(0 + .member_contrast_arbitrary | coupleID:diaryday)",
      difference_indicator =
        ".member_contrast_arbitrary"
    )
  )
)

At the dyad level, the fitted model includes a shared time slope but no difference time slope. Thus, the two members' time random effects are identical at this level, with correlation +1 whenever the shared slope variance is non-zero; at zero variance, the correlation is undefined. Covariances involving the diary-day slope are therefore supplied entirely by the shared block.

The random-effect terms may be copied exactly from the model formula. Equivalent backend syntax is also recognized, such as (1 + time | group) and us(1 + time | group), or (0 + x || group) and diag(0 + x | group). A custom difference-indicator name is supplied literally. Difference slopes may be written in either interaction order or as a simple product inside I().

When a difference block is supplied and the fitted model frame retains the indicator columns, difference_indicator must assign -1 and +1 to the two arbitrary member positions consistently within each dyad. For composition-specific blocks, it must be zero where shared_indicator is zero.

For custom difference indicators supplied through block_pairings, the function checks whether both positions occur within each supported fitted grouping unit. It rejects coding when no group contains both positions, and warns when only some groups are one-sided, which can result from fitted-row filtering. It also warns when stable member assignments cannot be verified across repeated rows without a member identifier.

Value

An exchangeable_covariance object: a named list with one element per matched block pairing. Each element contains shared_term and difference_term (each a fitted term string or NULL), plus the member-level variance-covariance matrix in varcov and its standard-deviation/correlation representation in sdcor, with standard deviations on the diagonal and correlations off the diagonal. Covariance dimension labels member_1 and member_2 denote arbitrary exchangeable positions, not substantive roles or encodings. Caller-supplied outer names are preserved; unnamed pairings receive stable names pair_1, pair_2, and so on in resolved order. varcov and sdcor are matrices, except when posterior = "draws" for a brms model, where they are draw by coefficient by coefficient arrays. brms quantities are transformed draw by draw before summarizing.

What omitted blocks and terms mean

recover_exchangeable_covariance() only describes constraints that were already imposed when the model was fitted. It does not remove a block, set a variance to zero, or otherwise constrain the supplied model. Describe only the structure that was actually fitted.

If a term occurs in only one selected block, the function represents the missing coordinate as a structural zero:

• A term present only in the shared block has no difference component, so the two members have identical random effects for that term.

• A term present only in the difference block has no shared component, so the two members have equal-magnitude, opposite-sign random effects for that term.

Setting difference_block = NULL or shared_block = NULL applies the corresponding rule to the entire omitted block. This is valid only when that block is truly absent from the fitted model. Do not use NULL merely to ignore an existing block; the resulting back-transformation would be incorrect.

See the constrained-block example in the exchangeable APIM vignette.

Backend note

For glmmTMB, random effects in non-conditional components are ignored with a warning. See the vignette's extension to exchangeable random slopes section for the manual calculation.

In brms, cross-sectional and same-occasion partner dependence can be represented directly with unstr(time = member_position, gr = residual_group). With Gaussian outcomes, sigma ~ 1 supplies the common residual scale. For families such as Poisson and negative binomial, unstr() adds a latent residual scale and correlation on the linear-predictor scale. Here, member_position identifies the same two arbitrary positions within every group, and residual_group identifies dyads in cross-sectional data or dyad-occasions in longitudinal data. This direct specification applies when one covariance structure is sufficient. Separate composition-specific unstr() structures for mixed dyad types are not currently supported in a standard single-response brms model. For Gaussian mixed-dyad residual covariance, use glmmTMB. Shared/difference blocks remain relevant for higher-level random effects and can represent latent link-scale covariance in non-Gaussian models.

See Also

The exchangeable APIM vignette for the model specification, covariance derivation, and constrained-block example. Run vignette("apim", package = "dyadMLM") to open the installed version.

Examples

if (requireNamespace("glmmTMB", quietly = TRUE)) {
  example_data <- prepare_dyad_data(
    dyads_cross,
    dyad = coupleID,
    member = personID,
    role = gender,
    model_types = "none",
    # All three observed compositions in `dyads_cross` are detected and retained
    # by default. This example focuses on `female-female` dyads, so we restrict
    # the analysis here.
    keep_compositions = "female-female",
    seed = 123
  )

  model <- glmmTMB::glmmTMB(
    closeness ~ 1 +
      us(1 | coupleID) +
      us(0 + .member_contrast_arbitrary | coupleID),
    dispformula = ~ 0,
    data = example_data
  )

  recover_exchangeable_covariance(model)
}

Description

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

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.