Package {dyadMLM}

Title:	Tools for Dyadic Multilevel Models
Version:	0.3.0
Date:	2026-10-27
Description:	Provides tools for dyadic multilevel modeling with linear and generalized linear mixed-effects models. It validates and prepares long-format cross-sectional and intensive longitudinal data, including ecological momentary assessment designs, for distinguishable and exchangeable dyads. It also supports datasets containing multiple observed dyad compositions. It constructs composition-aware, model-ready variables for Actor-Partner Interdependence Models (APIMs), Dyad-Individual Models (DIMs), and Dyadic Score Models (DSMs). Prepared data can be used with model engines such as 'glmmTMB' and 'brms' for Gaussian and non-Gaussian outcomes, including counts, proportions, and skewed continuous responses. Post-estimation tools compare compatible fitted models and back-transform exchangeable sum-and-difference random-effect covariance structures into member-level quantities. Experimental predictive checks compare observed and simulated response variances and partner correlations in Gaussian and selected non-Gaussian models. Distributional diagnostics separately check response distributions and residual patterns, using simulations that retain partner and time dependence. The APIM and DSM specifications and their relationships follow Iida et al. (2018) <doi:10.1177/0265407517725407>; the multilevel sum-and-difference random-effects implementation for exchangeable dyads adapts del Rosario and West (2025) <doi:10.1177/25152459251351286>.
License:	GPL-3
Encoding:	UTF-8
RoxygenNote:	8.0.0
Suggests:	brms, glmmTMB (≥ 1.1.11), htmltools, knitr, lavaan, lifecycle, marginaleffects, qgam, rmarkdown, spelling, testthat (≥ 3.1.7)
Config/testthat/edition:	3
Imports:	cli, dplyr (≥ 1.1.1), pillar, rlang (≥ 0.4.0), tidyselect, tibble (≥ 3.0.0), withr (≥ 3.0.3)
VignetteBuilder:	knitr, rmarkdown
Depends:	R (≥ 4.1)
LazyData:	true
URL:	https://pascal-kueng.github.io/dyadMLM/, https://github.com/Pascal-Kueng/dyadMLM
BugReports:	https://github.com/Pascal-Kueng/dyadMLM/issues
Language:	en-US
NeedsCompilation:	no
Packaged:	2026-10-28 00:00:01 UTC; pascalkueng
Author:	Pascal Küng [aut, cre, cph]
Maintainer:	Pascal Küng <kueng.pascal@gmail.com>
Repository:	CRAN
Date/Publication:	2026-10-28 00:20:02 UTC

dyadMLM: Tools for Dyadic Multilevel Models

Description

dyadMLM validates and prepares cross-sectional and intensive longitudinal dyadic data for multilevel modeling. It creates model-ready variables for APIMs, DIMs, and DSMs and provides supporting post-estimation tools.

Main functions

• prepare_dyad_data() validates long-format dyadic data and creates model-ready variables.

• compare_nested_models() compares compatible nested glmmTMB models fitted to equivalent data.

• recover_exchangeable_covariance() back-transforms exchangeable shared/difference random-effect covariance structures into member-level quantities.

• simulate_dyad_responses() and check_partner_dependence() provide experimental predictive comparisons of response variances and partner correlations.

• check_dyad_residuals() and check_dyad_outcomes() compare residual and outcome distributions using the same simulations.

Example data

See dyads_cross and dyads_ild for Gaussian examples, and dyads_nbinom_cross and dyads_nbinom_ild for negative-binomial examples.

Getting started

See vignette("getting-started", package = "dyadMLM") for an overview. The APIM, DIM, and DSM vignettes provide model-specific examples.

Author(s)

Maintainer: Pascal Küng kueng.pascal@gmail.com (ORCID) [copyright holder]

Authors:

• Pascal Küng kueng.pascal@gmail.com (ORCID) [copyright holder]

See Also

Useful links:

• https://pascal-kueng.github.io/dyadMLM/

• https://github.com/Pascal-Kueng/dyadMLM

• Report bugs at https://github.com/Pascal-Kueng/dyadMLM/issues

Add actor and partner predictor columns

Description

Adds APIM-style actor and partner columns for the predictors recorded in an dyadMLM_data object. For uncentered predictors, this will create actor and partner versions of the raw predictor. For centered intensive longitudinal predictors, this will create actor and partner versions of the raw predictor and each recorded predictor component, such as the within-person and between-person components created by center_predictors(). Selected lag predictors also create actor and partner columns for eligible lagged components.

Usage

add_actor_partner_columns(data)

Arguments

Details

The function will use the predictor decomposition metadata stored in attr(data, "dyadMLM")$temporal_decompositions, so downstream code does not need to infer generated predictor columns from their names. It stores the constructed APIM columns in attr(data, "dyadMLM")$apim_predictors.

Value

A dyadMLM_data object with actor and partner predictor columns added and APIM predictor metadata recorded.

Add dyad-individual predictor columns

Description

Adds Dyad-Individual Model (DIM) style dyad-mean and within-dyad member-deviation columns for the predictors recorded in a dyadMLM_data object. DIM requires one final exchangeable dyad composition, which may result from pooling exchangeable compositions. Distinguishable dyads are not currently supported. For intensive longitudinal predictors decomposed by center_predictors(), raw predictors and within-person components are decomposed within each dyad-time occasion, while between-person components are decomposed once within each dyad. For raw predictors, the dyad-mean column is centered around the grand mean of dyad means, or dyad-occasion means in longitudinal data, while the within-dyad member-deviation column is the person's deviation from the uncentered dyad mean. Selected lag predictors additionally create lag-1 raw and within-person dyad-mean and within-dyad member-deviation columns.

Usage

add_dyad_individual_columns(data)

Arguments

Details

The function reads attr(data, "dyadMLM")$temporal_decompositions and stores the constructed DIM columns in attr(data, "dyadMLM")$dim_predictors.

Value

A dyadMLM_data object with dyad-mean and within-dyad member-deviation predictor columns added and DIM predictor metadata recorded.

Add dyadic-score model predictor columns and contrast

Description

Adds Dyadic Score Model (DSM) dyad-mean and signed dyad-difference columns for the predictors recorded in a dyadMLM_data object, together with a DSM role contrast coded +0.5 and -0.5. DSM differences follow the role order recorded in attr(data, "dyadMLM")$dsm_role_order. The supported DSM structure contains one distinguishable dyad composition; exchangeable dyads and multiple compositions are not supported.

Usage

add_dyadic_score_columns(data)

Arguments

Details

For ILD predictors, raw and within-person scores are computed within dyad-time and between-person scores within dyad. Raw dyad means are grand-mean centered. Both partners' predictor values are required for each score pair. Selected lag predictors additionally create lag-1 raw and within-person dyad-mean and signed-difference columns.

Constructed predictor columns are recorded in attr(data, "dyadMLM")$dsm_predictors, and the contrast column name is recorded in attr(data, "dyadMLM")$dsm_role_contrast_column.

Value

A dyadMLM_data object with dyad-mean and signed dyad-difference predictor columns and a DSM role contrast added.

Add lagged temporal predictor columns

Description

Adds eligible lag-1 raw, within-person, and APIM GMC source columns. Values are matched at exactly time - 1; gaps are not bridged. Stable between-person components are not lagged.

Usage

add_temporal_lag_columns(data)

Arguments

Value

A dyadMLM_data object with lagged temporal predictor columns and updated predictor metadata.

Assign arbitrary member roles within dyads

Description

Creates a member-level lookup with one arbitrary role per observed member. The assignment is stable across longitudinal rows because it is made once per group x member.

Usage

assign_arbitrary_member_roles(data, group_name, member_name, seed = NULL)

Arguments

Value

A data frame with group_name, member_name, and .dy_arbitrary_role.

Create a canonical dyad composition label

Description

Sorts role labels before pasting them so composition labels do not depend on row order.

Usage

canonical_composition(roles, sep = dyad_composition_sep)

composition_role_label(composition, role, sep = dyad_composition_role_sep)

Arguments

Value

A single composition label.

Center predictor variables for dyadic models

Description

Adds centered predictor columns to a dyadMLM_data object. It currently supports two-level temporal centering for intensive longitudinal predictors: a within-person component and a between-person component. The original predictor remains available as a raw component for model-specific column construction. It can also add GMC source components for numeric APIM predictors. For two-level temporal centering, the between-person component is centered around the grand mean of person means, not the grand mean of all observed rows. This gives each person equal weight even when people have different numbers of observed measurement occasions.

Usage

center_predictors(data, add_apim_gmc_predictors = FALSE)

Arguments

Details

The function uses the structural metadata stored by prepare_dyad_data().

Value

A dyadMLM_data object with centered predictor columns added and updated predictor metadata.

Check outcome distributions against model simulations

Description

A model should reproduce the distribution of the observed outcomes. This check compares the observed outcome distribution and summaries (response SD after subtracting predictions, largest absolute deviation from prediction, number of zeros) with complete datasets from simulate_dyad_responses(). This helps identify mismatches in the model's assumptions, such as the distribution family, dispersion, or zero inflation. For PIT residual patterns, also run check_dyad_residuals() with the same simulations. For guidance on how to read the plots, see plot.dyadMLM_outcome_check().

Usage

check_dyad_outcomes(
  simulations,
  dyad = NULL,
  role = NULL,
  member = NULL,
  plot = TRUE,
  ask = NULL,
  panels = TRUE,
  data = NULL
)

Arguments

Value

Invisibly returns a dyadMLM_outcome_check list with role-specific summary matrices and outcome distributions in compositions. Matrices have one column per dataset (observed first, then all simulated). Optionally, suppress plotting with plot = FALSE and draw the saved result later with plot(result), which also takes ask and panels.

Scope

Complete simulated datasets keep the fitted partner and time dependence and any modeled differences in variability between roles. Unlike DHARMa's default (since version 0.5.0), simulations redraw all random effects. Dyadic models often fix the residual variance near zero (dispformula = ~ 0) and model the partners' residual variances and covariance as random effects. If held fixed, these random effects would reproduce the data almost exactly. Parameters stay fixed at their fitted values, so parameter uncertainty is not included. These are descriptive checks, not significance tests, and differ from DHARMa's plots and tests (see the distribution-checks article). They do not directly assess serial correlation.

Predicted outcomes are mean predictions with random effects set to zero (see simulate_dyad_responses()). The checks use the same fitted rows and families as simulate_dyad_responses() and weight all observations equally. Check partner correlations with check_partner_dependence().

See Also

check_dyad_residuals(), check_partner_dependence()

Examples

example_data <- prepare_dyad_data(
  dyads_cross,
  dyad = coupleID,
  member = personID,
  model_types = "none",
  seed = 123
)
# Exchangeable dyads: the two random-effect blocks represent the partners'
# residual variances and covariance.
model <- glmmTMB::glmmTMB(
  closeness ~ gender + provided_support +
    us(1 | coupleID) + us(0 + .member_contrast_arbitrary | coupleID),
  dispformula = ~ 0,
  data = example_data
)
simulations <- simulate_dyad_responses(
  model,
  nsim = 100, # Low for a quick example. Use at least 1,000.
  seed = 123
)
check_dyad_outcomes(
  simulations,
  dyad = coupleID,
  role = gender,
  ask = FALSE
)

Check residual distributions and patterns

Description

A model should reproduce how responses are distributed around its predictions. This check finds where each observed outcome falls among its simulated values (PIT residuals), using complete datasets from simulate_dyad_responses(). This helps identify mismatches in the model's assumptions, such as the distribution family, dispersion, or predictor effects. For checks on the outcome scale, also run check_dyad_outcomes() with the same simulations. For guidance on how to read the plots, see plot.dyadMLM_residual_check().

Usage

check_dyad_residuals(
  simulations,
  dyad = NULL,
  role = NULL,
  member = NULL,
  predictors = NULL,
  seed = 123,
  plot = TRUE,
  ask = NULL,
  panels = TRUE,
  data = NULL,
  quartile_lines = c("bins", "qgam")
)

Arguments

Value

Invisibly returns a dyadMLM_residual_check list with the pit matrix (after any centring) and calculated summaries in compositions. pit has one row per fitted observation and one column per dataset (observed first, then all simulated). Optionally, suppress plotting with plot = FALSE and draw the saved result later with plot(result), which also takes ask and panels.

Scope

Complete simulated datasets keep the fitted partner and time dependence and any modeled differences in variability between roles. Unlike DHARMa's default (since version 0.5.0), simulations redraw all random effects. Dyadic models often fix the residual variance near zero (dispformula = ~ 0) and model the partners' residual variances and covariance as random effects. If held fixed, these random effects would reproduce the data almost exactly. Parameters stay fixed at their fitted values, so parameter uncertainty is not included. These are descriptive checks, not significance tests, and differ from DHARMa's plots and tests (see the distribution-checks article). They do not directly assess serial correlation.

Predicted outcomes are mean predictions with random effects set to zero (see simulate_dyad_responses()). The checks use the same fitted rows and families as simulate_dyad_responses() and weight all observations equally. Check partner correlations with check_partner_dependence().

Method and credit

The internal calculation follows the simulation-based PIT approach used in Florian Hartig's DHARMa package and the randomized quantile residual method of Dunn and Smyth (1996). We implement the definition directly and keep the uniform 0–1 scale.

At each fitted row, rank the observed and all simulated outcomes together, breaking ties randomly. For rank r among m datasets, use (r - U) / m, where U is uniform between 0 and 1.

When the model freely estimates an overall location (an intercept or separate role means) and is not ordinal, PIT is centred automatically: each dataset's PIT is transformed to normal scores, centred on its median across all fitted rows, and transformed back. Panels then show residuals relative to this overall location, so all panels can show less of an overall shift. Centring also partly offsets how fitting makes departures look smaller.

Global envelopes use extreme-rank-length ordering (Myllymäki et al., 2017).

References

Dunn, P. K., & Smyth, G. K. (1996). Randomized quantile residuals. Journal of Computational and Graphical Statistics, 5(3), 236–244. doi:10.2307/1390802

Fasiolo, M., Wood, S. N., Zaffran, M., Nedellec, R., & Goude, Y. (2021). qgam: Bayesian nonparametric quantile regression modeling in R. Journal of Statistical Software, 100(9), 1–31. doi:10.18637/jss.v100.i09

Hartig, F. (2026). DHARMa: Residual diagnostics for hierarchical (multi-level/mixed) regression models (Version 0.5.0) [R package]. doi:10.32614/CRAN.package.DHARMa

Myllymäki, M., Mrkvička, T., Grabarnik, P., Seijo, H., & Hahn, U. (2017). Global envelope tests for spatial processes. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 79(2), 381–404. doi:10.1111/rssb.12172

See Also

check_dyad_outcomes(), check_partner_dependence()

Examples

example_data <- prepare_dyad_data(
  dyads_cross,
  dyad = coupleID,
  member = personID,
  model_types = "none",
  seed = 123
)
# Exchangeable dyads: the two random-effect blocks represent the partners'
# residual variances and covariance.
model <- glmmTMB::glmmTMB(
  closeness ~ gender + provided_support +
    us(1 | coupleID) + us(0 + .member_contrast_arbitrary | coupleID),
  dispformula = ~ 0,
  data = example_data
)
simulations <- simulate_dyad_responses(
  model,
  nsim = 200, # Low for a quick example. Use at least 1,000.
  seed = 123
)
check_dyad_residuals(
  simulations,
  dyad = coupleID,
  role = gender,
  predictors = "provided_support",
  ask = FALSE
)
# Save a pooled check without drawing, then plot each figure separately.
result <- check_dyad_residuals(
  simulations,
  role = NULL,
  predictors = "provided_support",
  plot = FALSE
)
plot(result, panels = FALSE, ask = FALSE)

Check whether a fitted model reproduces partner dependence

Description

A dyadic model should reproduce how much responses vary and how strongly partners' responses are related. This check compares the observed response variances and partner correlations with those of responses simulated from the fitted model. This helps identify mismatches in the model's assumptions. This check is especially useful when a less restricted model does not converge and so cannot be used for model comparison. For repeated occasions, supply member and time.

Usage

check_partner_dependence(
  simulations,
  dyad,
  role = NULL,
  member = NULL,
  time = 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. Rows with missing dyad IDs or roles and incomplete dyads are omitted with a warning. The object includes omission counts and settings, and can be saved and plotted later. With repeated occasions, compositions has a level column ("between" or "within"), and each label ends with its level. Within-level rows count the dyads that share at least one occasion.

Reading the plots

Blue 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 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 the same information using dyad averages and half-differences:

• 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.

An observed value far from most simulated values may indicate that the model does not reproduce that feature of the data well. Even if simulated and observed values agree, this does not show 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 separately can reveal differences hidden by pooling. With few dyads in a composition, the simulated ranges are wide, so its panels lack sensitivity. Flags (observed values outside the middle 95%, marked * when printed) can occur by chance, especially when checking several summaries. They are a starting point for further investigation, not necessarily a reason to reject the model. The covariance-pooling study illustrates detection and false alarms when checking each composition.

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

Repeated occasions

With repeated occasions, partners can be related in two ways: in their mean responses, and in being above or below their own means on the same occasions. A model should reproduce both. With member and time, every dataset (observed and simulated) is split into two levels after any model-centring. Each level is checked like cross-sectional data:

• Between: each member's mean over all their fitted occasions, paired within the dyad. Each dyad counts once.

• Within: each response minus its member's mean, paired with the partner's deviation on the same occasion. Each shared occasion counts once.

An occasion without the partner still counts in the member's mean.

The within check is related to repeated-measures correlation (Bakdash & Marusich, 2017).

The longitudinal study examines detection and false alarms. With 50 or fewer dyads, the check often missed omitted partner covariance.

Known limitations and caveats:

• A misfit at one level can also impact the other level's summaries, especially with few occasions (due to member means including occasion-level variance) or nonlinear links. Thus, both levels should be interpreted together.

• Check dependence across occasions (e.g., autocorrelation) separately. If the model omits temporal dependence, its covariance estimates can be biased, and these checks may not reveal it.

• These checks can be misleading when the model includes lagged outcomes as predictors.

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

Bakdash, J. Z., & Marusich, L. R. (2017). Repeated measures correlation. Frontiers in Psychology, 8, 456. doi:10.3389/fpsyg.2017.00456

See Also

check_dyad_residuals(), check_dyad_outcomes()

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 whether pooling fits the data. Each composition is checked
# separately, with role-specific summaries for distinguishable compositions.
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
)

plot(check, ask = FALSE, panels = TRUE)
print(check)

# For repeated occasions, also supply `member` and `time`.
example_diary_data <- prepare_dyad_data(
  dyads_ild[dyads_ild$coupleID <= 80, ], # Couples 1 to 80 are female-male.
  dyad = coupleID,
  member = personID,
  role = gender,
  time = diaryday,
  model_types = "none",
  seed = 123
)

diary_model <- glmmTMB::glmmTMB(
  closeness ~ 0 + gender + gender:diaryday +

    # In this example, we omit the correlation between partners' random
    # intercepts on level 2 with a diagonal covariance structure. This is
    # sometimes done when a model does not otherwise converge cleanly.
    diag(0 + gender | coupleID) +

    # Keep the full structure for the level 1 residuals.
    us(0 + gender | coupleID:diaryday),
  dispformula = ~ 0,
  data = example_diary_data
)

# The between level flags the omitted partner correlation.
check_partner_dependence(
  simulate_dyad_responses(diary_model, nsim = 100, seed = 123),
  dyad = coupleID,
  role = gender,
  member = personID,
  time = diaryday,
  # Supply the fitting data because personID is not in the model formula.
  data = example_diary_data
)

Compare nested glmmTMB models fitted to equivalent data

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.

Because the test requires maximum-likelihood fits, models fitted with priors in glmmTMB::glmmTMB() are rejected. Their estimates are maximum a posteriori, not maximum likelihood.

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

Collect dyadMLM-generated columns

Description

Returns a normalized, one-row-per-column view of the generated columns recorded in a dyadMLM attribute. Shared semantics and display text are supplied by generated_column_spec_lookup().

Usage

dyad_generated_columns(meta)

Arguments

Value

A tibble with one row per generated column. The lag column is 0 for contemporaneous columns and 1 for lag-1 columns.

Example Gaussian cross-sectional dyadic data

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.

Example Gaussian intensive longitudinal dyadic data

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.

Example negative-binomial cross-sectional dyadic data

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.

Example negative-binomial intensive longitudinal dyadic data

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.

Extract exchangeable random-effect blocks from a fitted model

Description

Normalizes the random-effect coefficients and fitted covariance parameters needed by recover_exchangeable_covariance() while keeping backend-specific work in two small adapters.

Usage

extract_exchangeable_residual_blocks(model)

Arguments

Value

A list containing the model backend, one normalized record per random-effect block, and one grouping-factor ID per fitted row. Every block record contains group, coefficients, correlated, term, and an estimate/draw-by-coefficients covariance array.

Format counted group identifiers for validation messages

Description

Converts a vector of dyad or group identifiers into text that includes both the number of groups and a compact list of their identifiers.

Usage

format_group_count(groups, singular = "dyad", plural = "dyads", max = 10)

Arguments

Value

A single character string.

Format group identifiers for validation messages

Description

Converts a vector of dyad or group identifiers into a compact comma-separated string for use in validation errors and warnings.

Usage

format_group_list(groups, max = 10)

Arguments

Value

A single character string.

Infer dyad compositions

Description

Builds a dyad-level summary of role compositions from a validated dyadMLM_data object.

Usage

infer_dyad_compositions(
  data,
  seed = NULL,
  keep_compositions = NULL,
  set_exchangeable_compositions = NULL,
  pool_compositions = NULL,
  short_colnames = TRUE,
  include_arbitrary_member_contrast = FALSE
)

Arguments

Value

A dyadMLM_data object with added .composition and .composition_role factor columns, ⁠.is_*⁠ numeric indicator columns, composition-specific numeric ⁠.member_contrast_*⁠ columns coded -1 and 1 for the two members of each matching composition and 0 otherwise, and dyad composition metadata.

Create safe suffixes for generated dyadMLM columns

Description

Create safe suffixes for generated dyadMLM columns

Usage

make_dyad_suffixes(
  labels,
  label_type = "labels",
  rename_hint = "role or composition labels"
)

Arguments

Value

A named character vector. Names are the original labels; values are sanitized column-name suffixes.

Read outcome check plots

Description

How to read the pages drawn by check_dyad_outcomes() and how to follow up flags. plot() redraws a check saved with plot = FALSE.

Usage

## S3 method for class 'dyadMLM_outcome_check'
plot(x, ask = NULL, panels = TRUE, ...)

Arguments

Value

x, invisibly.

Outcome check pages

The first row compares outcome distributions. Ordinal responses use bars for all defined categories, with separate 95% ranges per category, so one red point outside can occur by chance. Other outcomes show the proportion at or below each value (ECDF) for the observations and up to 30 simulated datasets. These simulated curves are not an envelope, so short stretches of the observed curve outside them can occur by chance.

The remaining rows compare the response SD (standard deviation of outcome minus prediction), the largest absolute deviation from prediction, and, when the composition has simulated zeros (or observed zeros, for count and Tweedie families), the number of zeros. The SD includes random effects, so it does not isolate the model's residual variance or dispersion parameter. For complete cross-sectional dyads with distinct roles, these SDs match the role SDs of check_partner_dependence() with its default response = "model-centred". Ordinal summaries depend on the category scores ⁠1, 2, ..., K⁠. Use a tall plotting window to keep all rows readable.

Views and compositions

Start with everyone pooled (role = NULL), then also inspect meaningful roles, such as gender, even when the model treats partners as exchangeable: role-specific problems show only there. Use predictors of check_dyad_residuals() to examine other measured characteristics. For separate plots and summaries by a person-level categorical characteristic, use that variable as role.

Each composition gets one pooled column for same-role partners, or separate columns for distinct roles. Supplying roles does not change the model's assumptions. Observations whose composition or role is unknown, for example because the partner's row is missing and data was not supplied, are omitted with a warning.

Reading flags

A flag is a red curve or point outside its blue envelope or range, or a red line outside its dashed limits. Each panel alone rarely flags a correct model, but flags add up across panels, roles, and predictor pages, so look for a consistent pattern, not any single flag. The simulation study reports how often correct and mismatched models were flagged.

With about 40 dyads, most mismatches studied were hard to detect, so the absence of flags does not show a good fit. Panels suggest possible problems, not their cause. Neither a flag nor its absence shows whether estimates or standard errors are wrong. After a flag, fit the model the pattern suggests and compare the estimates and standard errors you report. Mismatches showed mainly in these panels of both checks:

• Poisson fitted to overdispersed counts: mean PIT distance, number of zeros, QQ.

• Extra zeros: QQ and PIT quartiles (the number of zeros less often).

• Omitted role means: QQ and PIT quartiles across predicted outcomes, by role.

• Different role SDs: response SD and mean PIT distance, by role.

• Spread rising with a predictor: largest absolute deviation, mean PIT distance, QQ, and PIT quartiles across that predictor.

• Omitted partner effect: PIT quartiles across the partner's predictor (supply it in predictors of check_dyad_residuals()).

See Also

check_dyad_outcomes(), plot.dyadMLM_residual_check()

Plot a saved partner-dependence check

Description

Draws the comparison plots without repeating simulations.

Usage

## S3 method for class 'dyadMLM_partner_check'
plot(x, ask = NULL, panels = TRUE, ...)

Arguments

Details

See check_partner_dependence() for interpretation, panel layouts, and technical details.

Value

Invisibly, x.

Examples

example_data <- prepare_dyad_data(
  dyads_cross[dyads_cross$coupleID <= 40, ],
  dyad = coupleID,
  member = personID,
  model_types = "none",
  seed = 123
)

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

simulations <- simulate_dyad_responses(
  model,
  nsim = 50,
  seed = 123
)

check <- check_partner_dependence(
  simulations,
  dyad = coupleID,
  role = gender,
  data = example_data,
  plot = FALSE
)

plot(check, ask = FALSE)

Read residual check plots

Description

How to read the pages drawn by check_dyad_residuals() and how to follow up flags. plot() redraws a check saved with plot = FALSE.

Usage

## S3 method for class 'dyadMLM_residual_check'
plot(x, ask = NULL, panels = TRUE, ...)

Arguments

Value

x, invisibly.

Residual check pages

For each composition, the first page shows uniform QQ and PIT quartiles against predicted outcomes, like DHARMa's main plot. The second shows mean PIT distance. Each supplied predictor gets a page of PIT quartiles.

PIT (probability integral transform) residuals rank each outcome from 0 (low) to 1 (high) relative to its own simulated values. A suitable model should produce ranks roughly evenly spread between 0 and 1, allowing for the variation shown in blue. PIT distance is 2 * abs(PIT - 0.5). Blue envelopes are global: they are built so that the whole red curve of a correct model stays inside in a nominal 95% of datasets. This holds within each panel, not across roles, predictors, or pages. Because the model was fitted to these same data, departures often look smaller than they are, as in the settings studied by Robins et al. (2000) and in our simulations.

Predicted outcomes and numeric predictors with more than eight distinct values are grouped into up to eight bins per role, with at least about 20 observations per bin. As in DHARMa, these axes show ranks. Tick labels give the original values at the minimum, quartiles and maximum. Grey points show the observed PIT residuals. With four or more bins, PIT quartiles are smoothed across neighboring bins, using distances on the original scale. quartile_lines = "qgam" in check_dyad_residuals() replaces these smoothed quartiles with qgam lines. Because simulated datasets reuse the qgam tuning of the observed data, the 95% coverage of these envelopes is approximate. Categories, and predicted outcomes or numeric predictors with up to eight distinct values, get a box per value, as in DHARMa: a red box spans the observed 25th to 75th percentile, with the median in bold. Red points mark each observed quartile on its own blue range (25th left, median centre, 75th right). Flags come from these points, not from the box lines. All PIT quartile panels show full-model residuals, not partial effects. Bins and separate predictors can hide local patterns and interactions.

Views and compositions

Start with everyone pooled (role = NULL), then also inspect meaningful roles, such as gender, even when the model treats partners as exchangeable: role-specific problems show only there. Use predictors of check_dyad_residuals() to examine other measured characteristics. For separate plots and summaries by a person-level categorical characteristic, use that variable as role.

Each composition gets one pooled column for same-role partners, or separate columns for distinct roles. Supplying roles does not change the model's assumptions. Observations whose composition or role is unknown, for example because the partner's row is missing and data was not supplied, are omitted with a warning.

Reading flags

A flag is a red curve or point outside its blue envelope or range, or a red line outside its dashed limits. Each panel alone rarely flags a correct model, but flags add up across panels, roles, and predictor pages, so look for a consistent pattern, not any single flag. The simulation study reports how often correct and mismatched models were flagged.

With about 40 dyads, most mismatches studied were hard to detect, so the absence of flags does not show a good fit. Panels suggest possible problems, not their cause. Neither a flag nor its absence shows whether estimates or standard errors are wrong. After a flag, fit the model the pattern suggests and compare the estimates and standard errors you report. Mismatches showed mainly in these panels of both checks:

• Poisson fitted to overdispersed counts: mean PIT distance, number of zeros, QQ.

• Extra zeros: QQ and PIT quartiles (the number of zeros less often).

• Omitted role means: QQ and PIT quartiles across predicted outcomes, by role.

• Different role SDs: response SD and mean PIT distance, by role.

• Spread rising with a predictor: largest absolute deviation, mean PIT distance, QQ, and PIT quartiles across that predictor.

• Omitted partner effect: PIT quartiles across the partner's predictor (supply it in predictors of check_dyad_residuals()).

References

Robins, J. M., van der Vaart, A., & Ventura, V. (2000). Asymptotic distribution of P values in composite null models. Journal of the American Statistical Association, 95(452), 1143–1156. doi:10.1080/01621459.2000.10474310

See Also

check_dyad_residuals(), plot.dyadMLM_outcome_check()

Prepare dyadic data for multilevel models

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)

Print a summary of the partner-dependence predictive check object

Description

Shows each observed statistic with the middle 95% of its simulated values and marks observed values outside that range. Use plot.dyadMLM_partner_check() to view the comparisons.

Usage

## S3 method for class 'dyadMLM_partner_check'
print(x, digits = 3L, ...)

Arguments

Value

x, invisibly.

Print simulated response datasets

Description

Prints the number and type of simulated datasets.

Usage

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

Arguments

Value

x, invisible.

Print recovered exchangeable member-level covariance

Description

Print recovered exchangeable member-level covariance

Usage

## S3 method for class 'exchangeable_covariance'
print(x, representation = c("both", "varcov", "sdcor"), digits = 3L, ...)

Arguments

Value

x, invisibly.

Recover member-level residual covariance from exchangeable random-effect blocks

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.

The function checks whether both positions occur within each supported fitted grouping unit. It rejects coding when no group contains both positions. For custom difference indicators supplied through block_pairings, it also warns when only some groups are one-sided, which can result from fitted-row filtering, and 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)
}

Resolve user-supplied composition references

Description

Resolve user-supplied composition references

Usage

resolve_composition_references(references, observed_compositions, arg_name)

Arguments

Value

Canonical observed composition labels.

Resolve the DSM role order

Description

Checks the relationship between model_types, role, and dsm_role_order, validates the requested role order, and returns its stored representation.

Usage

resolve_dsm_role_order(dsm_role_order, model_types, has_role)

Arguments

Value

The validated, unnamed role order, or NULL.

Restore rows observed in the supplied data

Description

Removes rows added by temporarily_complete_dyad_occasions() and restores the original row order. Rows removed by an explicit preparation option, such as composition filtering, remain removed.

Usage

restore_observed_dyad_rows(data)

Arguments

Value

A dyadMLM_data object containing only originally observed rows.

Simulate response datasets for predictive checks

Description

If a model fits, data simulated from it should resemble the observed data. This function generates new response datasets for the same observations and predictors. Reuse them with check_partner_dependence(), check_dyad_residuals(), or check_dyad_outcomes() 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(), check_dyad_residuals(), and check_dyad_outcomes().

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() (glmmTMB 1.1.15.2 or later)

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 checks compare observed and simulated scores. The scores do not measure distances on an underlying continuous scale. Because each dyad contains only two people, glmmTMB tends to underestimate the dyad variance in ordinal models. Simulations then show less partner dependence than the data, so the partner check can flag a correct ordinal model, especially with many dyads.

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

Technical details

Each simulation draws new random effects, if the model has any, 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 (Gelman et al., 1996). If all grouping factors are dyads or nested within dyads (such as dyad-days), 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).

Partner-dependence checks (by default) and the response SD and largest absolute deviation of check_dyad_outcomes() 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.

References

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

Summarize prepared dyadic data

Description

Prints the dyadic structure followed by standard summaries for all columns.

Usage

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

Arguments

Value

Invisibly, the standard summary of all columns.

Temporarily complete observed dyad occasions

Description

Adds the missing member row when only one member is present at an observed dyad occasion. Structural identifiers and a resolved stable role are kept; all measured variables are missing on the added row. The added rows are used only while model-ready columns are constructed.

Usage

temporarily_complete_dyad_occasions(data)

Arguments

Value

A dyadMLM_data object containing temporary missing-member rows and an internal column recording the original row order.

Validate DIM compatibility

Description

Checks whether a dyadMLM_data object can be used for the currently supported undirected DIM construction. DIM requires one final exchangeable dyad composition, which may result from pooling exchangeable compositions. Distinguishable dyads are not currently supported.

Usage

validate_dim_compatibility(data)

Arguments

Value

Invisibly returns data when compatible.

Validate DSM compatibility

Description

Checks whether prepared data contain the single distinguishable dyad composition required by the DSM and whether its observed roles match the declared directional role order.

Usage

validate_dsm_compatibility(data)

Arguments

Value

Invisibly returns data when compatible.

Validate dyadic input data

Description

Checks whether data has a valid long-format dyadic structure and returns it as a tibble with an additional dyadMLM_data class. Cross-sectional data may contain at most one row per member within each dyad. Intensive longitudinal data may contain at most one row per member and measurement occasion within each dyad.

Usage

validate_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"),
  incomplete_dyads = c("error", "drop"),
  missing_role = c("error", "drop")
)

Arguments

Value

A tibble with class dyadMLM_data and metadata about the dyad, member, optional role, and optional time columns.