Nothing
Actor-Partner Interdependence Model (APIM)
In dyadMLM: Tools for Dyadic Multilevel Models
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
diagram_helper <- if (file.exists("diagram-helpers.Rinc")) {
"diagram-helpers.Rinc"
} else {
file.path("vignettes", "diagram-helpers.Rinc")
}
sys.source(diagram_helper, envir = knitr::knit_global())
diagram_device <- .vignette_diagram_device()
knitr::opts_chunk$set(
dev = diagram_device,
fig.ext = diagram_device,
dev.args = list(bg = "transparent")
)
library(dyadMLM)
has_glmmTMB <- requireNamespace("glmmTMB", quietly = TRUE)
apim_distinguishable_fitted_alt <-
"Fitted distinguishable APIM diagram unavailable."
apim_exchangeable_fitted_alt <-
"Fitted exchangeable APIM diagram unavailable."
This vignette focuses on Gaussian cross-sectional and intensive longitudinal
Actor-Partner Interdependence Models for distinguishable and exchangeable
dyads. The intensive longitudinal examples cover concurrent associations,
member-specific residual persistence, and a brief lagged-outcome extension.
For the broader package workflow and an overview of the available
model-specific vignettes, including the
Dyad-Individual Model and
Dyadic Score Model, see the
online package overview.
A vignette with non-Gaussian APIM examples is planned.
Cross-sectional APIMs
These examples use add_apim_gmc_predictors = TRUE, which retains raw APIM
columns and adds variants centered over all retained non-missing source values.
Using raw columns instead changes only the intercept reference. Do not include
both variants in a model with an intercept.
The distinguishable APIM
A conceptual example for distinguishable female-male dyads:
```r$, actor effects $a$, and partner effects $p$ can differ by the role of the outcome member (F and M), and the two outcome residuals covary within dyads.", fig.alt="Path diagram for a distinguishable cross-sectional APIM. Grand-mean-centered female and male predictors each have an actor path to their own outcome and a partner path to the other member's outcome. Female and male outcomes have separate intercepts, and their residuals covary."}
draw_apim_diagram("distinguishable", predictors_centered = TRUE)
For univariate MLM software like `glmmTMB`, this model is fitted in long format
with one outcome row per member, which can be visualized as:
```r
draw_apim_member_diagram("distinguishable", predictors_centered = TRUE)
Residual random-effects structure {#distinguishable-residual-structure}
For a distinguishable female-male dyad, the two members can have different
residual variances. In the notation below, the member role is written first and
$i$ indexes dyads. The within-dyad residual covariance block (shared across
dyads) is:
$$
\operatorname{Cov}
\begin{pmatrix}
\epsilon_{Fi} \
\epsilon_{Mi}
\end{pmatrix}
= \Sigma_{\epsilon}
= \begin{bmatrix}
\sigma_{\epsilon_F}^{2}
& \rho_{\epsilon_F\epsilon_M}\sigma_{\epsilon_F}\sigma_{\epsilon_M} \
\rho_{\epsilon_F\epsilon_M}\sigma_{\epsilon_F}\sigma_{\epsilon_M}
& \sigma_{\epsilon_M}^{2}
\end{bmatrix}
$$
And the full residual covariance matrix for all dyads (first three shown) is then block-diagonal:
$$
\Sigma_{\mathrm{model}}
= \begin{bmatrix}
\Sigma_{\epsilon} & 0 & 0 & \cdots \
0 & \Sigma_{\epsilon} & 0 & \cdots \
0 & 0 & \Sigma_{\epsilon} & \cdots \
\vdots & \vdots & \vdots & \ddots
\end{bmatrix}
$$
This structure is estimated with an unstructured random-effects block such as
us(0 + .is_female + .is_male | coupleID)
and dispformula = ~ 0.
Fitting the distinguishable APIM with glmmTMB
We first prepare the example data with dyadMLM::prepare_dyad_data():
apim_distinguishable_data <- dyadMLM::prepare_dyad_data(
data = dyads_cross,
dyad = coupleID,
member = personID,
role = gender,
predictors = provided_support,
model_types = "apim",
add_apim_gmc_predictors = TRUE,
# 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
)
print(apim_distinguishable_data, n=4)
The optional member contrast is used later for the restricted model. It does
not change the composition's distinguishable metadata or role indicators.
apim_distinguishable_model <- glmmTMB::glmmTMB(
closeness ~
# Gender-specific intercepts
0 +
.is_female +
.is_male +
# Gender-specific actor effects
.is_female:.provided_support_gmc_actor +
.is_male:.provided_support_gmc_actor +
# Gender-specific partner effects
.is_female:.provided_support_gmc_partner +
.is_male:.provided_support_gmc_partner +
# Dyad-level unstructured random effects represent the two partner
# residual variances and their covariance when dispformula = ~ 0.
# This is glmmTMB-specific syntax! `brms` uses different syntax.
us(0 +
.is_female +
.is_male
| coupleID)
, dispformula = ~ 0
, family = gaussian()
, data = apim_distinguishable_data
)
summary(apim_distinguishable_model)
apim_distinguishable_diagram_values <- .extract_apim_diagram_values(
apim_distinguishable_model,
type = "distinguishable"
)
apim_distinguishable_diagram_estimates <-
apim_distinguishable_diagram_values$estimates
apim_distinguishable_diagram_residuals <-
apim_distinguishable_diagram_values$residuals
apim_distinguishable_fitted_alt <- sprintf(
paste(
"Fitted distinguishable APIM. Female and male intercepts %.2f and %.2f;",
"actor effects %.2f and %.2f; partner effects %.2f and %.2f;",
"residual SDs %.2f and %.2f, with correlation %.2f."
),
apim_distinguishable_diagram_estimates[["intercept_female"]],
apim_distinguishable_diagram_estimates[["intercept_male"]],
apim_distinguishable_diagram_estimates[["actor_female"]],
apim_distinguishable_diagram_estimates[["actor_male"]],
apim_distinguishable_diagram_estimates[["partner_female"]],
apim_distinguishable_diagram_estimates[["partner_male"]],
apim_distinguishable_diagram_residuals[["sd_female"]],
apim_distinguishable_diagram_residuals[["sd_male"]],
apim_distinguishable_diagram_residuals[["correlation"]]
)
With the common centering, the two intercepts are the expected female and male
closeness scores when both partners' provided support equals the pooled sample
mean. The estimated coefficients map as follows:
draw_apim_diagram(
"distinguishable",
model = apim_distinguishable_model,
labels = c(predictor = "Provided support", outcome = "Closeness"),
predictors_centered = TRUE
)
The exchangeable APIM
Conceptually, the exchangeable APIM constrains several of the effects to be equal:
draw_apim_diagram("exchangeable", predictors_centered = TRUE)
Because the member labels are arbitrary, swapping members 1 and 2 does not
change the model.
To estimate this model in a univariate MLM framework, we can draw a conceptual
diagram as such:
draw_apim_member_diagram("exchangeable", predictors_centered = TRUE)
Modeling the residual random-effects structure {#exchangeable-residual-structure}
For a simple random-intercept structure, equal member variances and their
covariance can be specified directly. The specification becomes more difficult
when the same exchangeability constraints must cover several random slopes. A
single homogeneous structure would impose one variance across intercepts and
slopes, whereas separate structures would omit their correlations.
The shared/difference representation works for residuals and other random-effect
terms, including random slopes. Following
@delrosarioPracticalGuideSpecifying2025, dyadMLM::prepare_dyad_data()
generates an arbitrary member-difference column, named .member_contrast_*.
This contrast is +1 for one member and -1 for the other. The exchangeable
residual structure is represented by two separate random-effects terms: a
shared dyad random intercept and a random coefficient for this difference
column. Additional random slopes can be included in both blocks without
changing this logic.
We fit a simple exchangeable APIM and then use
dyadMLM::recover_exchangeable_covariance() to back-transform its
shared/difference representation into a member-level residual covariance matrix.
Fitting the restricted exchangeable APIM with glmmTMB
The restricted model omits gender-specific fixed effects and uses the
shared/difference residual representation. We reuse
apim_distinguishable_data, so both models use exactly the same prepared rows
and centering. set_exchangeable_compositions would instead reclassify the
composition when it should be treated as exchangeable throughout an analysis.
We then use the columns to fit the model as follows:
apim_exchangeable_model <- glmmTMB::glmmTMB(
closeness ~
# Pooled single intercept
1 +
# Pooled single actor and partner effects
.provided_support_gmc_actor +
.provided_support_gmc_partner +
# Residual variance covariance matrix via the shared/difference
# specification in two uncorrelated blocks
us(1 | coupleID) +
us(0 + .member_contrast_arbitrary | coupleID),
dispformula = ~ 0,
family = gaussian(),
data = apim_distinguishable_data
)
summary(apim_exchangeable_model)
We use dyadMLM::recover_exchangeable_covariance() to recover the interpretable
variance-covariance matrix:
backtransformed <- dyadMLM::recover_exchangeable_covariance(apim_exchangeable_model)
# residual variance-covariance and SD-correlation representations
print(backtransformed)
The back-transformation follows directly from the shared and member-difference
random effects. If $u_j$ is the shared effect for dyad $j$ and
$\widetilde{u}_j$ its member-difference effect, the two member effects are
$$
u_{1j} = u_j + \widetilde{u}j,
\qquad
u{2j} = u_j - \widetilde{u}_j.
$$
Because the two fitted blocks are independent,
$$
\operatorname{Var}(u_{1j}) = \operatorname{Var}(u_{2j})
= \operatorname{Var}(u_j) + \operatorname{Var}(\widetilde{u}j),
\qquad
\operatorname{Cov}(u{1j}, u_{2j})
= \operatorname{Var}(u_j) - \operatorname{Var}(\widetilde{u}_j).
$$
apim_exchangeable_diagram_values <- .extract_apim_diagram_values(
apim_exchangeable_model,
type = "exchangeable"
)
apim_exchangeable_diagram_estimates <-
apim_exchangeable_diagram_values$estimates
apim_exchangeable_diagram_residuals <-
apim_exchangeable_diagram_values$residuals
apim_exchangeable_fitted_alt <- sprintf(
paste(
"Fitted exchangeable APIM. Intercept %.2f, actor effect %.2f, partner",
"effect %.2f, common residual SD %.2f, and residual correlation %.2f."
),
apim_exchangeable_diagram_estimates[["intercept"]],
apim_exchangeable_diagram_estimates[["actor"]],
apim_exchangeable_diagram_estimates[["partner"]],
apim_exchangeable_diagram_residuals[["sd"]],
apim_exchangeable_diagram_residuals[["correlation"]]
)
With the common centering, the shared intercept is the expected closeness of
either member when both partners' provided support equals the pooled sample
mean. The output can now be mapped as follows:
draw_apim_diagram(
"exchangeable",
model = apim_exchangeable_model,
labels = c(predictor = "Provided support", outcome = "Closeness"),
predictors_centered = TRUE
)
Testing distinguishability {#testing-distinguishability}
Distinguishability can be evaluated by comparing a full model in which the two
roles may differ with a restricted exchangeable model. This comparison tests
the imposed equality constraints jointly. Here, they concern the fixed
intercepts, actor effects, partner effects, and residual variances.
Both models use the same prepared data object. The full model uses its retained
role indicators. The restricted model instead uses the opt-in arbitrary member
contrast for the exchangeability-constrained residual structure.
dyadMLM::compare_nested_models() verifies that both models use equivalent original
observations before performing the likelihood-ratio test:
dyadMLM::compare_nested_models(
apim_exchangeable_model,
apim_distinguishable_model
)
The test provides evidence against all restrictions jointly, but it does not
show which parameter differs. The helper also cannot determine whether the
models are mathematically nested. That remains a modeling requirement. The
usual chi-squared reference distribution may be unreliable when a tested
variance parameter lies on its boundary.
Checking partner dependence {#checking-partner-dependence}
The model comparison tests the restrictions jointly. Predictive checks show
whether the restricted model reproduces each role's variance and the partner
correlation. We supply the original roles to examine their variances separately.
::: {.no-invert-dark-mode}
apim_simulations <- dyadMLM::simulate_dyad_responses(
apim_exchangeable_model, seed = 123
)
dyadMLM::check_partner_dependence(
apim_simulations,
dyad = coupleID,
role = gender,
# Supply the fitting data because gender is not in the model formula.
data = apim_distinguishable_data
)
:::
The top row shows role-specific SDs and partner correlation. The bottom
row shows the same information using dyad averages and half the differences
between partners. Histograms show simulated values and red lines mark observed
values. A red line far from most simulations suggests a mismatch.
Here, the restricted model produces a lower partner correlation than observed.
It also overstates how much the partner differences vary across dyads. The
plots do not identify which restriction causes the mismatch.
By default, the check subtracts the model's fixed-effect predictions from
observed and simulated responses. Good agreement is often expected for variances
and correlations estimated from the same data. See the
function help
for details and interpretation. Simulation studies show how often the check
detects omitted partner dependence
and incorrect covariance pooling,
and how often it raises false alarms. For intensive longitudinal data, also supply
member and time. A separate study covers this case.
Checking residuals and outcomes {#checking-residuals-and-outcomes}
We reuse apim_simulations with
check_dyad_residuals()
and
check_dyad_outcomes().
The residual check uses PIT residuals, which rank each outcome from 0 (low)
to 1 (high) among its simulated values. The outcome check compares outcome
distributions and summaries.
We check by gender (role = gender), even though the restricted model treats
partners as exchangeable. We do this because role differences the model misses
can cancel out when all observations are pooled: if one role's outcomes are
higher than simulated and the other's lower, a pooled check can look fine. Here,
the pooled QQ plot (not shown) stays inside its envelope. Any other column that
gives each member a fixed role, such as older and younger partner, can be
checked the same way.
With predictors, we add a page of PIT quartiles across the partner's provided
support. If the model captures the partner effect, the quartiles stay near their
dashed lines at 0.25, 0.5, and 0.75 for all levels of partner support. With a
missing or misspecified partner effect, they move away from these lines as
partner support changes.
::: {.no-invert-dark-mode}
dyadMLM::check_dyad_residuals(
apim_simulations,
dyad = coupleID,
role = gender,
predictors = ".provided_support_gmc_partner",
data = apim_distinguishable_data
)
:::
::: {.no-invert-dark-mode}
dyadMLM::check_dyad_outcomes(
apim_simulations,
dyad = coupleID,
role = gender,
data = apim_distinguishable_data
)
:::
Look for flags (red curves outside their blue envelopes, or red lines outside
their dashed limits) that form a pattern across plots. Here, women's observed
closeness is higher than simulated and men's closeness is lower, in the QQ
plots, the PIT quartiles, and the outcome ECDFs. For men, the gap grows at high
predicted outcomes, as well as at high partner support. Mean PIT distance,
response SDs, and largest absolute deviations stay within their limits. These
panels mainly react to spread, so the model seems to reproduce how much
closeness varies around its predictions. The pattern suggests that the model
omits relevant gender differences, at least in mean closeness.
Adding gender-specific intercepts (0 + .is_female + .is_male instead of 1)
makes the actor estimate smaller and the partner estimate larger. Women provide
slightly more support than men here, so the omitted gender difference in mean
closeness had shifted both effects. Rechecking this model (not shown) flags only
men's PIT quartiles across predicted outcomes. The distinguishable APIM, which
also lets actor and partner effects and residual variances differ by gender,
shows no flags. The
simulation study
shows how often these panels flag correct and mismatched models.
If you know DHARMa: the uniform QQ
plot and the PIT quartiles across predicted outcomes use the same principle as
its residual plots and are read the same way. Two things differ:
- Our simulations redraw all random effects. DHARMa's default (since version
0.5.0) holds them at their fitted values and redraws only the residual
variation. In APIMs fitted with
dispformula = ~ 0, as here, the random
effects absorb all residual variation, so these simulations would reproduce
the data almost exactly. Redrawing all random effects also keeps partner
dependence in the reference.
- The blue envelopes are global envelopes [@myllymakiGlobalEnvelopeTests2017]
built from whole simulated datasets, which keep the fitted partner and time
dependence. They are an alternative to DHARMa's tests, which assume
independent residuals.
See why we do not use DHARMa directly.
Intensive longitudinal APIMs
For longitudinal APIMs, time-varying predictors are decomposed into
within-person and between-person components before actor and partner variables
are constructed [@bolgerIntensiveLongitudinalMethods2013;
@gistelinckMultilevelAutoregressiveDyadic2020]. The default "auto" selects
"2l" when both time and predictors are supplied. The within-person (cwp)
component captures occasion-specific deviations from each member's observed
mean, whereas the between-person (cbp) component captures each member's
observed mean relative to the sample grand mean.
Note that observed person means used to construct the between-person (cbp)
predictors can be unreliable when each member contributes few occasions,
which can bias between-person estimates
[@gottfredsonStraightforwardApproachCoping2019].
Concurrent ILD Gaussian APIM for distinguishable dyads
The decomposition above can be combined with role-specific effects. We first
retain the female-male distinction when preparing the data:
ild_distinguishable_data <- dyadMLM::prepare_dyad_data(
dyads_ild,
dyad = coupleID,
member = personID,
role = gender,
time = diaryday,
predictors = provided_support,
model_types = "apim",
keep_compositions = "female-male"
) |>
dplyr::mutate(
# we grand-mean center in this example to help convergence and
# interpretation
diaryday_gmc = diaryday - mean(diaryday)
)
The following model allows the intercepts, time trends, and within- and
between-person actor and partner effects to differ between female and male
members. Centering diaryday makes the role-specific intercepts refer to the
average study day:
ild_distinguishable_model <- glmmTMB::glmmTMB(
closeness ~
0 +
# Role-specific intercepts
.is_female +
.is_male +
# Role-specific time trends
.is_female:diaryday_gmc +
.is_male:diaryday_gmc +
# Role-specific within-person actor effects
.is_female:.provided_support_cwp_actor +
.is_male:.provided_support_cwp_actor +
# Role-specific within-person partner effects
.is_female:.provided_support_cwp_partner +
.is_male:.provided_support_cwp_partner +
# Role-specific between-person actor effects
.is_female:.provided_support_cbp_actor +
.is_male:.provided_support_cbp_actor +
# Role-specific between-person partner effects
.is_female:.provided_support_cbp_partner +
.is_male:.provided_support_cbp_partner +
# Stable dyad-level covariance
us(0 + .is_female + .is_male | coupleID) +
# Same-occasion covariance
us(0 + .is_female + .is_male |
coupleID:diaryday),
dispformula = ~ 0,
family = gaussian(),
data = ild_distinguishable_data
)
Adding role-specific residual AR(1) {#distinguishable-residual-ar1}
The same-occasion block above allows the partners' residuals to covary within a
day, but it does not describe persistence across days. In glmmTMB, we can
model independent persistent AR(1) components per member.
For distinguishable dyads, we add one term for each role while retaining the
same-occasion covariance. We spell out the complete model so the concurrent
and persistent components are visible together.
ild_distinguishable_ar_model <- glmmTMB::glmmTMB(
closeness ~
0 +
# Role-specific intercepts
.is_female +
.is_male +
# Role-specific time trends
.is_female:diaryday_gmc +
.is_male:diaryday_gmc +
# Role-specific within-person actor effects
.is_female:.provided_support_cwp_actor +
.is_male:.provided_support_cwp_actor +
# Role-specific within-person partner effects
.is_female:.provided_support_cwp_partner +
.is_male:.provided_support_cwp_partner +
# Role-specific between-person actor effects
.is_female:.provided_support_cbp_actor +
.is_male:.provided_support_cbp_actor +
# Role-specific between-person partner effects
.is_female:.provided_support_cbp_partner +
.is_male:.provided_support_cbp_partner +
# Stable dyad-level covariance
us(0 + .is_female + .is_male | coupleID) +
# Same-occasion covariance
us(0 + .is_female + .is_male |
coupleID:diaryday) +
# Role-specific residual persistence
ar1(0 + .is_female:factor(diaryday) | coupleID) +
ar1(0 + .is_male:factor(diaryday) | coupleID),
dispformula = ~ 0,
family = gaussian(),
data = ild_distinguishable_data,
# Non-default settings help this example converge.
control = glmmTMB::glmmTMBControl(profile = TRUE)
)
glmmTMB::VarCorr(ild_distinguishable_ar_model)
The two AR standard deviations and correlations describe each role's persistent
component. They are not correlations of the complete residual and do not
represent partner influence. The dyad-day block continues to describe the
remaining same-day variances and partner covariance.
Random slopes
For example, role-specific within-person actor random slopes can be added by
replacing the stable dyad-level block above with:
us(
0 +
.is_female +
.is_male +
.is_female:.provided_support_cwp_actor +
.is_male:.provided_support_cwp_actor
| coupleID
)
This block estimates the covariance among the two role-specific random
intercepts and actor slopes.
Concurrent ILD Gaussian APIM for exchangeable dyads
In longitudinal Gaussian exchangeable APIMs, the sum-and-difference
parametrization from
@delrosarioPracticalGuideSpecifying2025 can be extended to the dyad-occasion
level to represent same-occasion residual dependence.
We first prepare within-person (cwp) and between-person (cbp) actor and
partner predictors. We retain the female-female dyads as one substantively
exchangeable composition:
ild_apim_data <- dyadMLM::prepare_dyad_data(
dyads_ild,
dyad = coupleID,
member = personID,
role = gender,
time = diaryday,
predictors = provided_support,
model_types = "apim",
keep_compositions = "female-female",
seed = 123
)
print(ild_apim_data, n = 4)
The example below estimates same-day associations between support and closeness
and includes diaryday to adjust for a linear time trend. Its stable and
same-occasion covariance blocks impose exchangeability through paired shared
and member-contrast terms.
ild_apim_model <- glmmTMB::glmmTMB(
closeness ~
1 +
diaryday +
# Within-person actor and partner effects
.provided_support_cwp_actor +
.provided_support_cwp_partner +
# Between-person actor and partner effects
.provided_support_cbp_actor +
.provided_support_cbp_partner +
# Stable exchangeable dyad-level covariance
us(1 | coupleID) +
us(0 + .member_contrast_arbitrary | coupleID) +
# Same-occasion exchangeable covariance
us(1 | coupleID:diaryday) +
us(0 + .member_contrast_arbitrary | coupleID:diaryday)
, dispformula = ~ 0
, family = gaussian()
, data = ild_apim_data
)
summary(ild_apim_model)
Adding pooled residual AR(1) {#exchangeable-residual-ar1}
Under exchangeability, each member still has a separate time series, but the AR
standard deviation and correlation are pooled. We add one term grouped by the
member nested within the dyad:
ild_exchangeable_ar_model <- glmmTMB::glmmTMB(
closeness ~
1 +
diaryday +
# Within-person actor and partner effects
.provided_support_cwp_actor +
.provided_support_cwp_partner +
# Between-person actor and partner effects
.provided_support_cbp_actor +
.provided_support_cbp_partner +
# Stable exchangeable dyad-level covariance
us(1 | coupleID) +
us(0 + .member_contrast_arbitrary | coupleID) +
# Same-occasion exchangeable covariance
us(1 | coupleID:diaryday) +
us(0 + .member_contrast_arbitrary | coupleID:diaryday) +
# Pooled residual persistence across member series
ar1(0 + factor(diaryday) | coupleID:personID),
dispformula = ~ 0,
family = gaussian(),
data = ild_apim_data,
# Non-default settings help this example converge.
control = glmmTMB::glmmTMBControl(profile = TRUE)
)
glmmTMB::VarCorr(ild_exchangeable_ar_model)
The pooled AR parameters describe the persistent within-member component. The
dyad-day blocks continue to describe same-day partner covariance.
For the fitted AR model, we can recover the member-level covariance matrices for
both the stable dyad effects and the same-occasion residual dependence. The two
matched block pairs are returned separately:
recovered_covariance <- dyadMLM::recover_exchangeable_covariance(
ild_exchangeable_ar_model
)
print(recovered_covariance)
The cwp terms estimate actor and partner associations for occasion-specific
deviations from each member's usual support. The cbp terms estimate actor and
partner associations involving members' usual support levels. These are
concurrent associations. They do not by themselves represent temporal
carryover.
Extension to exchangeable random slopes {#exchangeable-random-slope-back-transformation}
For a focused random-slope demonstration, we branch from the concurrent base
model and replace its stable intercept blocks with shared and member-contrast
intercept-and-slope blocks. This advanced fit does not include the AR term
above, so its covariance estimates illustrate the parameterization
rather than complete recovery of the simulation.
ild_apim_random_slope_model <- update(
ild_apim_model,
formula = . ~ . -
us(1 | coupleID) -
us(0 + .member_contrast_arbitrary | coupleID) +
us(1 + .provided_support_cwp_actor | coupleID) +
us(0 +
.member_contrast_arbitrary +
.member_contrast_arbitrary:
.provided_support_cwp_actor
| coupleID),
# Non-default settings help this example converge.
control = glmmTMB::glmmTMBControl(
profile = TRUE,
optimizer = stats::optim,
optArgs = list(method = "BFGS")
)
)
random_slope_covariance <- dyadMLM::recover_exchangeable_covariance(
ild_apim_random_slope_model
)
print(random_slope_covariance)
The same shared/difference back-transformation described above applies
separately to every random slope. It also applies to covariances among the
random intercept, actor slope, and partner slope.
Testing random-effect constraints {#fitted-constraints-and-omitted-blocks}
The random-slope model above estimates both a shared and a member-contrast actor
random slope. We can first test a smaller model that omits only the actor random
slope from the member-contrast block. The member-contrast random intercept and
both same-occasion blocks remain in the model:
ild_apim_no_contrast_slope <- update(
ild_apim_random_slope_model,
formula = . ~ . -
us(0 +
.member_contrast_arbitrary +
.member_contrast_arbitrary:
.provided_support_cwp_actor
| coupleID) +
us(0 + .member_contrast_arbitrary | coupleID)
)
dyadMLM::compare_nested_models(
ild_apim_no_contrast_slope,
ild_apim_random_slope_model
)
Without the member-contrast slope, the two members have identical actor random
slopes at the stable dyad level. We tell the back-transformation which fitted
blocks represent this constraint.
Since we omitted terms, automatic matching is no longer possible and we
need to tell dyadMLM::recover_exchangeable_covariance() what blocks
belong together.
no_contrast_slope_covariance <- dyadMLM::recover_exchangeable_covariance(
ild_apim_no_contrast_slope,
block_pairings = list(
dyad = list(
shared_block =
"us(1 + .provided_support_cwp_actor | coupleID)",
difference_block =
"us(0 + .member_contrast_arbitrary | coupleID)",
difference_indicator =
".member_contrast_arbitrary"
)
)
)
print(no_contrast_slope_covariance, representation = "sdcor")
We can impose the stronger constraint by omitting the full member-contrast
block at the stable dyad level. The same-occasion member-contrast block again
remains in the model:
ild_apim_no_contrast_block <- update(
ild_apim_random_slope_model,
formula = . ~ . -
us(0 +
.member_contrast_arbitrary +
.member_contrast_arbitrary:
.provided_support_cwp_actor
| coupleID)
)
dyadMLM::compare_nested_models(
ild_apim_no_contrast_block,
ild_apim_random_slope_model
)
Here, both members have identical random intercepts and actor random slopes at
the stable dyad level. Because the full member-contrast block is absent, we
specify it as NULL:
no_contrast_block_covariance <- dyadMLM::recover_exchangeable_covariance(
ild_apim_no_contrast_block,
block_pairings = list(
dyad = list(
shared_block =
"us(1 + .provided_support_cwp_actor | coupleID)",
difference_block = NULL,
difference_indicator =
".member_contrast_arbitrary"
)
)
)
print(no_contrast_block_covariance, representation = "sdcor")
These are constraints on the stable dyad-level random effects, not on the
same-occasion residual structure. Because variance constraints lie on the
boundary of the parameter space, the usual chi-squared reference distribution
for the likelihood-ratio tests should be interpreted cautiously.
When carryover is the research question {#dynamic-models}
Residual AR(1) and lagged-outcome models answer different questions. The AR(1)
terms above retain the concurrent APIM mean model while describing residual
persistence. When temporal carryover is itself the research question, the
member's own and the partner's lagged outcomes can enter the mean model
[@gistelinckMultilevelAutoregressiveDyadic2020;
@delrosarioPracticalGuideSpecifying2025]. A partner-lag coefficient is an
association. Interpreting it as influence requires stronger assumptions.
Residual dependence should still be assessed after adding outcome lags.
By adding the outcome to predictors and selecting it with lag1_predictors,
dyadMLM::prepare_dyad_data() returns lag-1 raw and within-person scores alongside
the contemporaneous scores. Between-person scores are not lagged because they
describe stable differences between members.
Here, lag1_predictors creates the lagged actor and partner outcome columns:
ild_apim_data_dynamic <- dyadMLM::prepare_dyad_data(
dyads_ild,
dyad = coupleID,
member = personID,
role = gender,
time = diaryday,
predictors = closeness,
lag1_predictors = closeness,
model_types = "apim",
keep_compositions = "female-female",
seed = 123
)
This returns lag-1 raw and within-person closeness scores. Lags match exactly
diaryday - 1, so omitted diary days are not bridged.
These data were generated with residual AR(1), not a lagged-outcome mean model.
The following unevaluated model demonstrates the specification only; its
coefficients would not recover the data-generating process.
lagged_outcome_example <- glmmTMB::glmmTMB(
closeness ~ 1 +
# Own-outcome lag (stability association)
.closeness_actor_lag1 +
# Partner-outcome lag association
.closeness_partner_lag1 +
# Linear time trend
diaryday +
# Stable exchangeable dyad-level covariance
us(1 | coupleID) +
us(0 + .member_contrast_arbitrary | coupleID) +
# Same-day exchangeable dyad-level covariance
us(1 | coupleID:diaryday) +
us(0 + .member_contrast_arbitrary | coupleID:diaryday)
, dispformula = ~ 0
, family = gaussian()
, data = ild_apim_data_dynamic
# Non-default settings help this example converge.
, control = glmmTMB::glmmTMBControl(
profile = TRUE,
optimizer = stats::optim,
optArgs = list(method = "BFGS")
)
)
The model uses raw outcome lags. Within-person-centered alternatives end in
_cwp_actor_lag1 and _cwp_partner_lag1, but person-mean centering can bias
carryover downward in short panels [@hamakerCenterNotCenter2015;
@nickellBiasesDynamicModels1981]. Raw lags avoid that centering bias but can
remain related to stable member levels
[@gistelinckMultilevelAutoregressiveSmall2021]. For short panels, an LD-APIM can
model the initial outcomes jointly with those stable levels
[@gistelinckMultilevelAutoregressiveDyadic2020].
Independent member AR(1) terms do not represent cross-partner residual
carryover. A full residual VAR or RDSEM is a different model and is not directly
available through the glmmTMB interface used here
[@asparouhovComparisonModelsAnalysis2020; @mcneishPrimerTwoLevel2020].
From here, choose the model-specific vignette that matches the research
question:
- Dyad-Individual Model vignette for the exchangeable DIM
parameterization; or
- Dyadic Score Model vignette for the distinguishable DSM
parameterization.
Or return to the
online package overview.
References
::: {#refs}
:::
Try the dyadMLM package in your browser
Any scripts or data that you put into this service are public.
dyadMLM documentation built on Oct. 28, 2026, 1:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) diagram_helper <- if (file.exists("diagram-helpers.Rinc")) { "diagram-helpers.Rinc" } else { file.path("vignettes", "diagram-helpers.Rinc") } sys.source(diagram_helper, envir = knitr::knit_global()) diagram_device <- .vignette_diagram_device() knitr::opts_chunk$set( dev = diagram_device, fig.ext = diagram_device, dev.args = list(bg = "transparent") )
library(dyadMLM)
has_glmmTMB <- requireNamespace("glmmTMB", quietly = TRUE) apim_distinguishable_fitted_alt <- "Fitted distinguishable APIM diagram unavailable." apim_exchangeable_fitted_alt <- "Fitted exchangeable APIM diagram unavailable."
This vignette focuses on Gaussian cross-sectional and intensive longitudinal Actor-Partner Interdependence Models for distinguishable and exchangeable dyads. The intensive longitudinal examples cover concurrent associations, member-specific residual persistence, and a brief lagged-outcome extension.
For the broader package workflow and an overview of the available model-specific vignettes, including the Dyad-Individual Model and Dyadic Score Model, see the online package overview.
A vignette with non-Gaussian APIM examples is planned.
Cross-sectional APIMs
These examples use add_apim_gmc_predictors = TRUE, which retains raw APIM
columns and adds variants centered over all retained non-missing source values.
Using raw columns instead changes only the intercept reference. Do not include
both variants in a model with an intercept.
The distinguishable APIM
A conceptual example for distinguishable female-male dyads:
```r$, actor effects $a$, and partner effects $p$ can differ by the role of the outcome member (F and M), and the two outcome residuals covary within dyads.", fig.alt="Path diagram for a distinguishable cross-sectional APIM. Grand-mean-centered female and male predictors each have an actor path to their own outcome and a partner path to the other member's outcome. Female and male outcomes have separate intercepts, and their residuals covary."} draw_apim_diagram("distinguishable", predictors_centered = TRUE)
For univariate MLM software like `glmmTMB`, this model is fitted in long format
with one outcome row per member, which can be visualized as:
```r
draw_apim_member_diagram("distinguishable", predictors_centered = TRUE)
Residual random-effects structure {#distinguishable-residual-structure}
For a distinguishable female-male dyad, the two members can have different residual variances. In the notation below, the member role is written first and $i$ indexes dyads. The within-dyad residual covariance block (shared across dyads) is:
$$ \operatorname{Cov} \begin{pmatrix} \epsilon_{Fi} \ \epsilon_{Mi} \end{pmatrix} = \Sigma_{\epsilon} = \begin{bmatrix} \sigma_{\epsilon_F}^{2} & \rho_{\epsilon_F\epsilon_M}\sigma_{\epsilon_F}\sigma_{\epsilon_M} \ \rho_{\epsilon_F\epsilon_M}\sigma_{\epsilon_F}\sigma_{\epsilon_M} & \sigma_{\epsilon_M}^{2} \end{bmatrix} $$
And the full residual covariance matrix for all dyads (first three shown) is then block-diagonal:
$$ \Sigma_{\mathrm{model}} = \begin{bmatrix} \Sigma_{\epsilon} & 0 & 0 & \cdots \ 0 & \Sigma_{\epsilon} & 0 & \cdots \ 0 & 0 & \Sigma_{\epsilon} & \cdots \ \vdots & \vdots & \vdots & \ddots \end{bmatrix} $$
This structure is estimated with an unstructured random-effects block such as
us(0 + .is_female + .is_male | coupleID)
and dispformula = ~ 0.
Fitting the distinguishable APIM with glmmTMB
We first prepare the example data with dyadMLM::prepare_dyad_data():
apim_distinguishable_data <- dyadMLM::prepare_dyad_data( data = dyads_cross, dyad = coupleID, member = personID, role = gender, predictors = provided_support, model_types = "apim", add_apim_gmc_predictors = TRUE, # 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 ) print(apim_distinguishable_data, n=4)
The optional member contrast is used later for the restricted model. It does not change the composition's distinguishable metadata or role indicators.
apim_distinguishable_model <- glmmTMB::glmmTMB( closeness ~ # Gender-specific intercepts 0 + .is_female + .is_male + # Gender-specific actor effects .is_female:.provided_support_gmc_actor + .is_male:.provided_support_gmc_actor + # Gender-specific partner effects .is_female:.provided_support_gmc_partner + .is_male:.provided_support_gmc_partner + # Dyad-level unstructured random effects represent the two partner # residual variances and their covariance when dispformula = ~ 0. # This is glmmTMB-specific syntax! `brms` uses different syntax. us(0 + .is_female + .is_male | coupleID) , dispformula = ~ 0 , family = gaussian() , data = apim_distinguishable_data ) summary(apim_distinguishable_model)
apim_distinguishable_diagram_values <- .extract_apim_diagram_values( apim_distinguishable_model, type = "distinguishable" ) apim_distinguishable_diagram_estimates <- apim_distinguishable_diagram_values$estimates apim_distinguishable_diagram_residuals <- apim_distinguishable_diagram_values$residuals apim_distinguishable_fitted_alt <- sprintf( paste( "Fitted distinguishable APIM. Female and male intercepts %.2f and %.2f;", "actor effects %.2f and %.2f; partner effects %.2f and %.2f;", "residual SDs %.2f and %.2f, with correlation %.2f." ), apim_distinguishable_diagram_estimates[["intercept_female"]], apim_distinguishable_diagram_estimates[["intercept_male"]], apim_distinguishable_diagram_estimates[["actor_female"]], apim_distinguishable_diagram_estimates[["actor_male"]], apim_distinguishable_diagram_estimates[["partner_female"]], apim_distinguishable_diagram_estimates[["partner_male"]], apim_distinguishable_diagram_residuals[["sd_female"]], apim_distinguishable_diagram_residuals[["sd_male"]], apim_distinguishable_diagram_residuals[["correlation"]] )
With the common centering, the two intercepts are the expected female and male closeness scores when both partners' provided support equals the pooled sample mean. The estimated coefficients map as follows:
draw_apim_diagram( "distinguishable", model = apim_distinguishable_model, labels = c(predictor = "Provided support", outcome = "Closeness"), predictors_centered = TRUE )
The exchangeable APIM
Conceptually, the exchangeable APIM constrains several of the effects to be equal:
draw_apim_diagram("exchangeable", predictors_centered = TRUE)
Because the member labels are arbitrary, swapping members 1 and 2 does not change the model.
To estimate this model in a univariate MLM framework, we can draw a conceptual diagram as such:
draw_apim_member_diagram("exchangeable", predictors_centered = TRUE)
Modeling the residual random-effects structure {#exchangeable-residual-structure}
For a simple random-intercept structure, equal member variances and their covariance can be specified directly. The specification becomes more difficult when the same exchangeability constraints must cover several random slopes. A single homogeneous structure would impose one variance across intercepts and slopes, whereas separate structures would omit their correlations.
The shared/difference representation works for residuals and other random-effect
terms, including random slopes. Following
@delrosarioPracticalGuideSpecifying2025, dyadMLM::prepare_dyad_data()
generates an arbitrary member-difference column, named .member_contrast_*.
This contrast is +1 for one member and -1 for the other. The exchangeable
residual structure is represented by two separate random-effects terms: a
shared dyad random intercept and a random coefficient for this difference
column. Additional random slopes can be included in both blocks without
changing this logic.
We fit a simple exchangeable APIM and then use
dyadMLM::recover_exchangeable_covariance() to back-transform its
shared/difference representation into a member-level residual covariance matrix.
Fitting the restricted exchangeable APIM with glmmTMB
The restricted model omits gender-specific fixed effects and uses the
shared/difference residual representation. We reuse
apim_distinguishable_data, so both models use exactly the same prepared rows
and centering. set_exchangeable_compositions would instead reclassify the
composition when it should be treated as exchangeable throughout an analysis.
We then use the columns to fit the model as follows:
apim_exchangeable_model <- glmmTMB::glmmTMB( closeness ~ # Pooled single intercept 1 + # Pooled single actor and partner effects .provided_support_gmc_actor + .provided_support_gmc_partner + # Residual variance covariance matrix via the shared/difference # specification in two uncorrelated blocks us(1 | coupleID) + us(0 + .member_contrast_arbitrary | coupleID), dispformula = ~ 0, family = gaussian(), data = apim_distinguishable_data ) summary(apim_exchangeable_model)
We use dyadMLM::recover_exchangeable_covariance() to recover the interpretable
variance-covariance matrix:
backtransformed <- dyadMLM::recover_exchangeable_covariance(apim_exchangeable_model) # residual variance-covariance and SD-correlation representations print(backtransformed)
The back-transformation follows directly from the shared and member-difference random effects. If $u_j$ is the shared effect for dyad $j$ and $\widetilde{u}_j$ its member-difference effect, the two member effects are
$$ u_{1j} = u_j + \widetilde{u}j, \qquad u{2j} = u_j - \widetilde{u}_j. $$
Because the two fitted blocks are independent,
$$ \operatorname{Var}(u_{1j}) = \operatorname{Var}(u_{2j}) = \operatorname{Var}(u_j) + \operatorname{Var}(\widetilde{u}j), \qquad \operatorname{Cov}(u{1j}, u_{2j}) = \operatorname{Var}(u_j) - \operatorname{Var}(\widetilde{u}_j). $$
apim_exchangeable_diagram_values <- .extract_apim_diagram_values( apim_exchangeable_model, type = "exchangeable" ) apim_exchangeable_diagram_estimates <- apim_exchangeable_diagram_values$estimates apim_exchangeable_diagram_residuals <- apim_exchangeable_diagram_values$residuals apim_exchangeable_fitted_alt <- sprintf( paste( "Fitted exchangeable APIM. Intercept %.2f, actor effect %.2f, partner", "effect %.2f, common residual SD %.2f, and residual correlation %.2f." ), apim_exchangeable_diagram_estimates[["intercept"]], apim_exchangeable_diagram_estimates[["actor"]], apim_exchangeable_diagram_estimates[["partner"]], apim_exchangeable_diagram_residuals[["sd"]], apim_exchangeable_diagram_residuals[["correlation"]] )
With the common centering, the shared intercept is the expected closeness of either member when both partners' provided support equals the pooled sample mean. The output can now be mapped as follows:
draw_apim_diagram( "exchangeable", model = apim_exchangeable_model, labels = c(predictor = "Provided support", outcome = "Closeness"), predictors_centered = TRUE )
Testing distinguishability {#testing-distinguishability}
Distinguishability can be evaluated by comparing a full model in which the two roles may differ with a restricted exchangeable model. This comparison tests the imposed equality constraints jointly. Here, they concern the fixed intercepts, actor effects, partner effects, and residual variances.
Both models use the same prepared data object. The full model uses its retained role indicators. The restricted model instead uses the opt-in arbitrary member contrast for the exchangeability-constrained residual structure.
dyadMLM::compare_nested_models() verifies that both models use equivalent original
observations before performing the likelihood-ratio test:
dyadMLM::compare_nested_models( apim_exchangeable_model, apim_distinguishable_model )
The test provides evidence against all restrictions jointly, but it does not show which parameter differs. The helper also cannot determine whether the models are mathematically nested. That remains a modeling requirement. The usual chi-squared reference distribution may be unreliable when a tested variance parameter lies on its boundary.
Checking partner dependence {#checking-partner-dependence}
The model comparison tests the restrictions jointly. Predictive checks show whether the restricted model reproduces each role's variance and the partner correlation. We supply the original roles to examine their variances separately.
::: {.no-invert-dark-mode}
apim_simulations <- dyadMLM::simulate_dyad_responses( apim_exchangeable_model, seed = 123 ) dyadMLM::check_partner_dependence( apim_simulations, dyad = coupleID, role = gender, # Supply the fitting data because gender is not in the model formula. data = apim_distinguishable_data )
:::
The top row shows role-specific SDs and partner correlation. The bottom row shows the same information using dyad averages and half the differences between partners. Histograms show simulated values and red lines mark observed values. A red line far from most simulations suggests a mismatch.
Here, the restricted model produces a lower partner correlation than observed. It also overstates how much the partner differences vary across dyads. The plots do not identify which restriction causes the mismatch.
By default, the check subtracts the model's fixed-effect predictions from
observed and simulated responses. Good agreement is often expected for variances
and correlations estimated from the same data. See the
function help
for details and interpretation. Simulation studies show how often the check
detects omitted partner dependence
and incorrect covariance pooling,
and how often it raises false alarms. For intensive longitudinal data, also supply
member and time. A separate study covers this case.
Checking residuals and outcomes {#checking-residuals-and-outcomes}
We reuse apim_simulations with
check_dyad_residuals()
and
check_dyad_outcomes().
The residual check uses PIT residuals, which rank each outcome from 0 (low)
to 1 (high) among its simulated values. The outcome check compares outcome
distributions and summaries.
We check by gender (role = gender), even though the restricted model treats
partners as exchangeable. We do this because role differences the model misses
can cancel out when all observations are pooled: if one role's outcomes are
higher than simulated and the other's lower, a pooled check can look fine. Here,
the pooled QQ plot (not shown) stays inside its envelope. Any other column that
gives each member a fixed role, such as older and younger partner, can be
checked the same way.
With predictors, we add a page of PIT quartiles across the partner's provided
support. If the model captures the partner effect, the quartiles stay near their
dashed lines at 0.25, 0.5, and 0.75 for all levels of partner support. With a
missing or misspecified partner effect, they move away from these lines as
partner support changes.
::: {.no-invert-dark-mode}
dyadMLM::check_dyad_residuals( apim_simulations, dyad = coupleID, role = gender, predictors = ".provided_support_gmc_partner", data = apim_distinguishable_data )
:::
::: {.no-invert-dark-mode}
dyadMLM::check_dyad_outcomes( apim_simulations, dyad = coupleID, role = gender, data = apim_distinguishable_data )
:::
Look for flags (red curves outside their blue envelopes, or red lines outside their dashed limits) that form a pattern across plots. Here, women's observed closeness is higher than simulated and men's closeness is lower, in the QQ plots, the PIT quartiles, and the outcome ECDFs. For men, the gap grows at high predicted outcomes, as well as at high partner support. Mean PIT distance, response SDs, and largest absolute deviations stay within their limits. These panels mainly react to spread, so the model seems to reproduce how much closeness varies around its predictions. The pattern suggests that the model omits relevant gender differences, at least in mean closeness.
Adding gender-specific intercepts (0 + .is_female + .is_male instead of 1)
makes the actor estimate smaller and the partner estimate larger. Women provide
slightly more support than men here, so the omitted gender difference in mean
closeness had shifted both effects. Rechecking this model (not shown) flags only
men's PIT quartiles across predicted outcomes. The distinguishable APIM, which
also lets actor and partner effects and residual variances differ by gender,
shows no flags. The
simulation study
shows how often these panels flag correct and mismatched models.
If you know DHARMa: the uniform QQ plot and the PIT quartiles across predicted outcomes use the same principle as its residual plots and are read the same way. Two things differ:
- Our simulations redraw all random effects. DHARMa's default (since version
0.5.0) holds them at their fitted values and redraws only the residual
variation. In APIMs fitted with
dispformula = ~ 0, as here, the random effects absorb all residual variation, so these simulations would reproduce the data almost exactly. Redrawing all random effects also keeps partner dependence in the reference. - The blue envelopes are global envelopes [@myllymakiGlobalEnvelopeTests2017] built from whole simulated datasets, which keep the fitted partner and time dependence. They are an alternative to DHARMa's tests, which assume independent residuals.
See why we do not use DHARMa directly.
Intensive longitudinal APIMs
For longitudinal APIMs, time-varying predictors are decomposed into
within-person and between-person components before actor and partner variables
are constructed [@bolgerIntensiveLongitudinalMethods2013;
@gistelinckMultilevelAutoregressiveDyadic2020]. The default "auto" selects
"2l" when both time and predictors are supplied. The within-person (cwp)
component captures occasion-specific deviations from each member's observed
mean, whereas the between-person (cbp) component captures each member's
observed mean relative to the sample grand mean.
Note that observed person means used to construct the between-person (cbp)
predictors can be unreliable when each member contributes few occasions,
which can bias between-person estimates
[@gottfredsonStraightforwardApproachCoping2019].
Concurrent ILD Gaussian APIM for distinguishable dyads
The decomposition above can be combined with role-specific effects. We first retain the female-male distinction when preparing the data:
ild_distinguishable_data <- dyadMLM::prepare_dyad_data( dyads_ild, dyad = coupleID, member = personID, role = gender, time = diaryday, predictors = provided_support, model_types = "apim", keep_compositions = "female-male" ) |> dplyr::mutate( # we grand-mean center in this example to help convergence and # interpretation diaryday_gmc = diaryday - mean(diaryday) )
The following model allows the intercepts, time trends, and within- and
between-person actor and partner effects to differ between female and male
members. Centering diaryday makes the role-specific intercepts refer to the
average study day:
ild_distinguishable_model <- glmmTMB::glmmTMB( closeness ~ 0 + # Role-specific intercepts .is_female + .is_male + # Role-specific time trends .is_female:diaryday_gmc + .is_male:diaryday_gmc + # Role-specific within-person actor effects .is_female:.provided_support_cwp_actor + .is_male:.provided_support_cwp_actor + # Role-specific within-person partner effects .is_female:.provided_support_cwp_partner + .is_male:.provided_support_cwp_partner + # Role-specific between-person actor effects .is_female:.provided_support_cbp_actor + .is_male:.provided_support_cbp_actor + # Role-specific between-person partner effects .is_female:.provided_support_cbp_partner + .is_male:.provided_support_cbp_partner + # Stable dyad-level covariance us(0 + .is_female + .is_male | coupleID) + # Same-occasion covariance us(0 + .is_female + .is_male | coupleID:diaryday), dispformula = ~ 0, family = gaussian(), data = ild_distinguishable_data )
Adding role-specific residual AR(1) {#distinguishable-residual-ar1}
The same-occasion block above allows the partners' residuals to covary within a day, but it does not describe persistence across days. In glmmTMB, we can model independent persistent AR(1) components per member.
For distinguishable dyads, we add one term for each role while retaining the same-occasion covariance. We spell out the complete model so the concurrent and persistent components are visible together.
ild_distinguishable_ar_model <- glmmTMB::glmmTMB( closeness ~ 0 + # Role-specific intercepts .is_female + .is_male + # Role-specific time trends .is_female:diaryday_gmc + .is_male:diaryday_gmc + # Role-specific within-person actor effects .is_female:.provided_support_cwp_actor + .is_male:.provided_support_cwp_actor + # Role-specific within-person partner effects .is_female:.provided_support_cwp_partner + .is_male:.provided_support_cwp_partner + # Role-specific between-person actor effects .is_female:.provided_support_cbp_actor + .is_male:.provided_support_cbp_actor + # Role-specific between-person partner effects .is_female:.provided_support_cbp_partner + .is_male:.provided_support_cbp_partner + # Stable dyad-level covariance us(0 + .is_female + .is_male | coupleID) + # Same-occasion covariance us(0 + .is_female + .is_male | coupleID:diaryday) + # Role-specific residual persistence ar1(0 + .is_female:factor(diaryday) | coupleID) + ar1(0 + .is_male:factor(diaryday) | coupleID), dispformula = ~ 0, family = gaussian(), data = ild_distinguishable_data, # Non-default settings help this example converge. control = glmmTMB::glmmTMBControl(profile = TRUE) ) glmmTMB::VarCorr(ild_distinguishable_ar_model)
The two AR standard deviations and correlations describe each role's persistent component. They are not correlations of the complete residual and do not represent partner influence. The dyad-day block continues to describe the remaining same-day variances and partner covariance.
Random slopes
For example, role-specific within-person actor random slopes can be added by replacing the stable dyad-level block above with:
us( 0 + .is_female + .is_male + .is_female:.provided_support_cwp_actor + .is_male:.provided_support_cwp_actor | coupleID )
This block estimates the covariance among the two role-specific random intercepts and actor slopes.
Concurrent ILD Gaussian APIM for exchangeable dyads
In longitudinal Gaussian exchangeable APIMs, the sum-and-difference parametrization from @delrosarioPracticalGuideSpecifying2025 can be extended to the dyad-occasion level to represent same-occasion residual dependence.
We first prepare within-person (cwp) and between-person (cbp) actor and
partner predictors. We retain the female-female dyads as one substantively
exchangeable composition:
ild_apim_data <- dyadMLM::prepare_dyad_data( dyads_ild, dyad = coupleID, member = personID, role = gender, time = diaryday, predictors = provided_support, model_types = "apim", keep_compositions = "female-female", seed = 123 ) print(ild_apim_data, n = 4)
The example below estimates same-day associations between support and closeness
and includes diaryday to adjust for a linear time trend. Its stable and
same-occasion covariance blocks impose exchangeability through paired shared
and member-contrast terms.
ild_apim_model <- glmmTMB::glmmTMB( closeness ~ 1 + diaryday + # Within-person actor and partner effects .provided_support_cwp_actor + .provided_support_cwp_partner + # Between-person actor and partner effects .provided_support_cbp_actor + .provided_support_cbp_partner + # Stable exchangeable dyad-level covariance us(1 | coupleID) + us(0 + .member_contrast_arbitrary | coupleID) + # Same-occasion exchangeable covariance us(1 | coupleID:diaryday) + us(0 + .member_contrast_arbitrary | coupleID:diaryday) , dispformula = ~ 0 , family = gaussian() , data = ild_apim_data ) summary(ild_apim_model)
Adding pooled residual AR(1) {#exchangeable-residual-ar1}
Under exchangeability, each member still has a separate time series, but the AR standard deviation and correlation are pooled. We add one term grouped by the member nested within the dyad:
ild_exchangeable_ar_model <- glmmTMB::glmmTMB( closeness ~ 1 + diaryday + # Within-person actor and partner effects .provided_support_cwp_actor + .provided_support_cwp_partner + # Between-person actor and partner effects .provided_support_cbp_actor + .provided_support_cbp_partner + # Stable exchangeable dyad-level covariance us(1 | coupleID) + us(0 + .member_contrast_arbitrary | coupleID) + # Same-occasion exchangeable covariance us(1 | coupleID:diaryday) + us(0 + .member_contrast_arbitrary | coupleID:diaryday) + # Pooled residual persistence across member series ar1(0 + factor(diaryday) | coupleID:personID), dispformula = ~ 0, family = gaussian(), data = ild_apim_data, # Non-default settings help this example converge. control = glmmTMB::glmmTMBControl(profile = TRUE) ) glmmTMB::VarCorr(ild_exchangeable_ar_model)
The pooled AR parameters describe the persistent within-member component. The dyad-day blocks continue to describe same-day partner covariance.
For the fitted AR model, we can recover the member-level covariance matrices for both the stable dyad effects and the same-occasion residual dependence. The two matched block pairs are returned separately:
recovered_covariance <- dyadMLM::recover_exchangeable_covariance( ild_exchangeable_ar_model ) print(recovered_covariance)
The cwp terms estimate actor and partner associations for occasion-specific
deviations from each member's usual support. The cbp terms estimate actor and
partner associations involving members' usual support levels. These are
concurrent associations. They do not by themselves represent temporal
carryover.
Extension to exchangeable random slopes {#exchangeable-random-slope-back-transformation}
For a focused random-slope demonstration, we branch from the concurrent base model and replace its stable intercept blocks with shared and member-contrast intercept-and-slope blocks. This advanced fit does not include the AR term above, so its covariance estimates illustrate the parameterization rather than complete recovery of the simulation.
ild_apim_random_slope_model <- update( ild_apim_model, formula = . ~ . - us(1 | coupleID) - us(0 + .member_contrast_arbitrary | coupleID) + us(1 + .provided_support_cwp_actor | coupleID) + us(0 + .member_contrast_arbitrary + .member_contrast_arbitrary: .provided_support_cwp_actor | coupleID), # Non-default settings help this example converge. control = glmmTMB::glmmTMBControl( profile = TRUE, optimizer = stats::optim, optArgs = list(method = "BFGS") ) ) random_slope_covariance <- dyadMLM::recover_exchangeable_covariance( ild_apim_random_slope_model ) print(random_slope_covariance)
The same shared/difference back-transformation described above applies separately to every random slope. It also applies to covariances among the random intercept, actor slope, and partner slope.
Testing random-effect constraints {#fitted-constraints-and-omitted-blocks}
The random-slope model above estimates both a shared and a member-contrast actor random slope. We can first test a smaller model that omits only the actor random slope from the member-contrast block. The member-contrast random intercept and both same-occasion blocks remain in the model:
ild_apim_no_contrast_slope <- update( ild_apim_random_slope_model, formula = . ~ . - us(0 + .member_contrast_arbitrary + .member_contrast_arbitrary: .provided_support_cwp_actor | coupleID) + us(0 + .member_contrast_arbitrary | coupleID) ) dyadMLM::compare_nested_models( ild_apim_no_contrast_slope, ild_apim_random_slope_model )
Without the member-contrast slope, the two members have identical actor random slopes at the stable dyad level. We tell the back-transformation which fitted blocks represent this constraint.
Since we omitted terms, automatic matching is no longer possible and we
need to tell dyadMLM::recover_exchangeable_covariance() what blocks
belong together.
no_contrast_slope_covariance <- dyadMLM::recover_exchangeable_covariance( ild_apim_no_contrast_slope, block_pairings = list( dyad = list( shared_block = "us(1 + .provided_support_cwp_actor | coupleID)", difference_block = "us(0 + .member_contrast_arbitrary | coupleID)", difference_indicator = ".member_contrast_arbitrary" ) ) ) print(no_contrast_slope_covariance, representation = "sdcor")
We can impose the stronger constraint by omitting the full member-contrast block at the stable dyad level. The same-occasion member-contrast block again remains in the model:
ild_apim_no_contrast_block <- update( ild_apim_random_slope_model, formula = . ~ . - us(0 + .member_contrast_arbitrary + .member_contrast_arbitrary: .provided_support_cwp_actor | coupleID) ) dyadMLM::compare_nested_models( ild_apim_no_contrast_block, ild_apim_random_slope_model )
Here, both members have identical random intercepts and actor random slopes at
the stable dyad level. Because the full member-contrast block is absent, we
specify it as NULL:
no_contrast_block_covariance <- dyadMLM::recover_exchangeable_covariance( ild_apim_no_contrast_block, block_pairings = list( dyad = list( shared_block = "us(1 + .provided_support_cwp_actor | coupleID)", difference_block = NULL, difference_indicator = ".member_contrast_arbitrary" ) ) ) print(no_contrast_block_covariance, representation = "sdcor")
These are constraints on the stable dyad-level random effects, not on the same-occasion residual structure. Because variance constraints lie on the boundary of the parameter space, the usual chi-squared reference distribution for the likelihood-ratio tests should be interpreted cautiously.
When carryover is the research question {#dynamic-models}
Residual AR(1) and lagged-outcome models answer different questions. The AR(1) terms above retain the concurrent APIM mean model while describing residual persistence. When temporal carryover is itself the research question, the member's own and the partner's lagged outcomes can enter the mean model [@gistelinckMultilevelAutoregressiveDyadic2020; @delrosarioPracticalGuideSpecifying2025]. A partner-lag coefficient is an association. Interpreting it as influence requires stronger assumptions. Residual dependence should still be assessed after adding outcome lags.
By adding the outcome to predictors and selecting it with lag1_predictors,
dyadMLM::prepare_dyad_data() returns lag-1 raw and within-person scores alongside
the contemporaneous scores. Between-person scores are not lagged because they
describe stable differences between members.
Here, lag1_predictors creates the lagged actor and partner outcome columns:
ild_apim_data_dynamic <- dyadMLM::prepare_dyad_data( dyads_ild, dyad = coupleID, member = personID, role = gender, time = diaryday, predictors = closeness, lag1_predictors = closeness, model_types = "apim", keep_compositions = "female-female", seed = 123 )
This returns lag-1 raw and within-person closeness scores. Lags match exactly
diaryday - 1, so omitted diary days are not bridged.
These data were generated with residual AR(1), not a lagged-outcome mean model. The following unevaluated model demonstrates the specification only; its coefficients would not recover the data-generating process.
lagged_outcome_example <- glmmTMB::glmmTMB( closeness ~ 1 + # Own-outcome lag (stability association) .closeness_actor_lag1 + # Partner-outcome lag association .closeness_partner_lag1 + # Linear time trend diaryday + # Stable exchangeable dyad-level covariance us(1 | coupleID) + us(0 + .member_contrast_arbitrary | coupleID) + # Same-day exchangeable dyad-level covariance us(1 | coupleID:diaryday) + us(0 + .member_contrast_arbitrary | coupleID:diaryday) , dispformula = ~ 0 , family = gaussian() , data = ild_apim_data_dynamic # Non-default settings help this example converge. , control = glmmTMB::glmmTMBControl( profile = TRUE, optimizer = stats::optim, optArgs = list(method = "BFGS") ) )
The model uses raw outcome lags. Within-person-centered alternatives end in
_cwp_actor_lag1 and _cwp_partner_lag1, but person-mean centering can bias
carryover downward in short panels [@hamakerCenterNotCenter2015;
@nickellBiasesDynamicModels1981]. Raw lags avoid that centering bias but can
remain related to stable member levels
[@gistelinckMultilevelAutoregressiveSmall2021]. For short panels, an LD-APIM can
model the initial outcomes jointly with those stable levels
[@gistelinckMultilevelAutoregressiveDyadic2020].
Independent member AR(1) terms do not represent cross-partner residual
carryover. A full residual VAR or RDSEM is a different model and is not directly
available through the glmmTMB interface used here
[@asparouhovComparisonModelsAnalysis2020; @mcneishPrimerTwoLevel2020].
From here, choose the model-specific vignette that matches the research question:
- Dyad-Individual Model vignette for the exchangeable DIM parameterization; or
- Dyadic Score Model vignette for the distinguishable DSM parameterization.
Or return to the online package overview.
References
::: {#refs} :::
Try the dyadMLM package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.