One possible workflow

1. Define dyads, distinguishability, and research questions
2. Prepare the data and validate the dyadic structure
3. Describe and visualize measures and observed relationships
4. Estimate the models
5. Verify model adequacy, distinguishability (if applicable), and robustness
6. Visualize, interpret, and report results

One possible workflow

1. Define dyads, distinguishability, and research questions
2. Prepare the data and validate the dyadic structure
3. Describe and visualize measures and observed relationships
4. Estimate the models
5. Verify model adequacy, distinguishability (if applicable), and robustness
6. Visualize, interpret, and report results

The Dataset

Same as before:

• 36 physically inactive romantic couples
• Both partners intended to become more physically active
• Smartphone-based planning interventions and JITAIs
• Daily diaries over 55 days

Aggregated for cross-sectional analyses

Measurement - predictor

Self-efficacy

Rough English translation of the original German item
“I am confident that I can be physically active tomorrow, even if it becomes difficult.”

Adapted from Sniehotta et al. (2005)

0 = not at all true today
5 = completely true today

Measurement - outcome

Moderate to vigorous physical activity (MVPA)

Rough English translation of the original German item
“How many minutes did you spend engaging in moderate-to-vigorous physical activity today?”

Adapted from Amireault & Godin (2015)

• Alone: _____ minutes
• Together with your partner: _____ minutes

Total MVPA = solo MVPA + joint MVPA.

Workflow

1. Define dyads, distinguishability, and research questions
2. Prepare the data and validate the dyadic structure
3. Describe and visualize measures and observed relationships
4. Estimate the models
5. Verify model adequacy, distinguishability, and robustness
6. Visualize, interpret, and report results

Load and inspect the raw data

library(dplyr) # for easier data-wrangling, exports functions like `select()`

raw_dyad_data <- readRDS("dyadic-person-means.rds") |>
  dplyr::select(
    couple_id, person_id, gender, efficacy, total_mvpa
  )

# usually use head(raw_dyad_data), slide_table is a custom function for html tables on these slides
slide_table(raw_dyad_data[1:4,])

Create the required format and validate the data

For the distinguishable APIM we need to add:

• efficacy_actor and efficacy_partner columns
• Dummy coded numeric variables for the two roles

Create the required format and validate the data

We can use the dyadMLM package to help with data preparation and validation.

library(dyadMLM)

apim_distinguishable_data <- dyadMLM::prepare_dyad_data( # exported function from dyadMLM
  data = raw_dyad_data,
  dyad = couple_id,
  member = person_id,
  role = gender,
  predictors = efficacy,
  model_types = 'apim',
  # if we want to fit an exchangeable model from the same data we can:
  include_arbitrary_member_contrast = TRUE,
  add_apim_gmc_predictors = TRUE,
  seed = 123
)

print(apim_distinguishable_data)

Create the required format and validate the data Manual code →

Looking at the columns immediately relevant to us:

head(apim_distinguishable_data) |>
  dplyr::select(
    couple_id, person_id, gender, .is_female, .is_male, .efficacy_gmc, .efficacy_gmc_actor, .efficacy_gmc_partner, total_mvpa
  ) |>
  slide_table()

Workflow

1. Define dyads, distinguishability, and research questions
2. Prepare the data and validate the dyadic structure
3. Describe and visualize measures and observed relationships
4. Estimate the models
5. Verify model adequacy, distinguishability, and robustness
6. Visualize, interpret, and report results

Create a wide view for role-specific exploration

The models will NOT use this dataset

apim_distinguishable_wide_data <- apim_distinguishable_data |>
  dplyr::select(couple_id, gender, efficacy, total_mvpa) |>
  tidyr::pivot_wider(
    names_from = gender,
    values_from = c(efficacy, total_mvpa)
  )

slide_table(head(apim_distinguishable_wide_data))

Describe focal variables per role

library(report)

apim_distinguishable_wide_data |>
  dplyr::select(-couple_id) |>
  report::report_table() |>
  slide_table()

Visualize focal variable by partner role

library(ggplot2)
library(see)

plot_role_histograms <- function(data, group, var) {
  var_label <- rlang::as_label(rlang::enquo(var))

  ggplot(
    data,
    aes(x = {{ var }}, y = after_stat(density))
  ) +
    geom_histogram(
      bins = 15, boundary = 0,
      fill = "steelblue", color = "white", alpha = 0.75
    ) +
    geom_density(
      bounds = c(0, Inf), color = "navy", linewidth = 1.1
    ) +
    facet_wrap(vars({{ group }}), nrow = 1) +
    labs(x = var_label, y = "Density") +
    see::theme_modern(base_size = 16)
}

Visualize focal variable by partner role - Efficacy

plot_role_histograms(apim_distinguishable_data, gender, efficacy)

Visualize focal variable by partner role - MVPA

plot_role_histograms(apim_distinguishable_data, gender, total_mvpa)

Observed dyadic correlations

The wide data ensure that each couple contributes once.

library(correlation)
library(insight)

apim_distinguishable_wide_data |>
  correlation::correlation(
    select = c("efficacy_female", "total_mvpa_female"),
    select2 = c("efficacy_male", "total_mvpa_male")
  ) |>
  summary() |>
  insight::print_md(footer = "")

Visualize partner similarity

plot_partner_similarity <- function(data, x, y, id, title, n_mahalanobis = 3) {
  x_name <- rlang::as_name(rlang::enquo(x))
  y_name <- rlang::as_name(rlang::enquo(y))
  xy <- cbind(data[[x_name]], data[[y_name]])
  labels <- data |>
    dplyr::mutate(
      .distance = mahalanobis(xy, colMeans(xy), cov(xy))
    ) |>
    dplyr::distinct({{ id }}, .keep_all = TRUE) |>
    dplyr::slice_max(.distance, n = n_mahalanobis, with_ties = FALSE)

  ggplot(data, aes(x = {{ x }}, y = {{ y }})) +
    geom_point(color = "steelblue") +
    geom_smooth(method = "lm", se = FALSE, color = "navy") +
    geom_text(
      data = labels,
      aes(label = {{ id }}),
      vjust = -0.6, check_overlap = TRUE, size = 4
    ) +
    labs(title = title, x = "Female", y = "Male") +
    scale_y_continuous(expand = expansion(mult = c(0.05, 0.12))) +
    see::theme_modern(base_size = 16)
}

Visualize partner similarity

n_mahalanobis <- 4

see::plots(
  plot_partner_similarity(
    apim_distinguishable_wide_data,
    efficacy_female, efficacy_male, couple_id, "Efficacy",
    n_mahalanobis = n_mahalanobis
  ),
  plot_partner_similarity(
    apim_distinguishable_wide_data,
    total_mvpa_female, total_mvpa_male, couple_id, "Daily MVPA",
    n_mahalanobis = n_mahalanobis
  )
)

Workflow

1. Define dyads, distinguishability, and research questions
2. Prepare the data and validate the dyadic structure
3. Describe and visualize measures and observed relationships
4. Estimate the models
5. Verify model adequacy, distinguishability, and robustness
6. Visualize, interpret, and report results

Distinguishable APIM

The APIM starts with two member-specific equations

For each dyad \(i\), each member with role \(r \in \{\mathrm{M},\mathrm{F}\}\) has their own predictor \(X^{(\mathrm{A})}_{ri}\) and their partner’s predictor \(X^{(\mathrm{P})}_{ri}\):

\[ \begin{gathered} Y_{\mathrm{M},i} = b_{0,\mathrm{M}} + \color{#1769AA}{a_{\mathrm{M}} \times X^{(\mathrm{A})}_{\mathrm{M},i}} + \color{#7B3FA1}{p_{\mathrm{M}} \times X^{(\mathrm{P})}_{\mathrm{M},i}} + \epsilon_{\mathrm{M},i} \\ Y_{\mathrm{F},i} = b_{0,\mathrm{F}} + \color{#1769AA}{a_{\mathrm{F}} \times X^{(\mathrm{A})}_{\mathrm{F},i}} + \color{#7B3FA1}{p_{\mathrm{F}} \times X^{(\mathrm{P})}_{\mathrm{F},i}} + \epsilon_{\mathrm{F},i} \end{gathered} \]

- Actor effects - Partner effects

In long format, each equation applies to one member row of the dyad.

How can we combine these two equations into one model?

Role indicators combine the equations

The indicators \(D_{\mathrm{M},ri}\) and \(D_{\mathrm{F},ri}\) identify the male or female row:

\[ \begin{aligned} Y_{ri} ={}& \underbrace{ b_{0,\mathrm{M}} \times D_{\mathrm{M},ri} + b_{0,\mathrm{F}} \times D_{\mathrm{F},ri} }_{\text{role-specific intercepts}} \\[0.35em] &+ \color{#1769AA}{ \underbrace{ a_{\mathrm{M}} \times D_{\mathrm{M},ri} \times X^{(\mathrm{A})}_{ri} + a_{\mathrm{F}} \times D_{\mathrm{F},ri} \times X^{(\mathrm{A})}_{ri} }_{\text{role-specific actor effects}} } \\[0.35em] &+ \color{#7B3FA1}{ \underbrace{ p_{\mathrm{M}} \times D_{\mathrm{M},ri} \times X^{(\mathrm{P})}_{ri} + p_{\mathrm{F}} \times D_{\mathrm{F},ri} \times X^{(\mathrm{P})}_{ri} }_{\text{role-specific partner effects}} } + \epsilon_{ri}. \end{aligned} \]

Residuals remain correlated within dyads

The two member rows from the same dyad retain a joint residual structure:

\[ \operatorname{Cov} \begin{pmatrix} \epsilon_{\mathrm{F}i} \\ \epsilon_{\mathrm{M}i} \end{pmatrix} = \boldsymbol{\Sigma}_{\epsilon} = \begin{bmatrix} \sigma_{\epsilon_\mathrm{F}}^{2} & \rho_{\epsilon_\mathrm{F}\epsilon_\mathrm{M}} \sigma_{\epsilon_\mathrm{F}}\sigma_{\epsilon_\mathrm{M}} \\ \rho_{\epsilon_\mathrm{F}\epsilon_\mathrm{M}} \sigma_{\epsilon_\mathrm{F}}\sigma_{\epsilon_\mathrm{M}} & \sigma_{\epsilon_\mathrm{M}}^{2} \end{bmatrix} \]

Fitting the distinguishable APIM in glmmTMB

library(glmmTMB)

apim_distinguishable_model <- glmmTMB(
  total_mvpa ~
    # Role-specific intercepts
    0 + .is_female + .is_male +

    # Role-specific actor effects
    .is_female:.efficacy_gmc_actor +
    .is_male:.efficacy_gmc_actor +

    # Role-specific partner effects
    .is_female:.efficacy_gmc_partner +
    .is_male:.efficacy_gmc_partner +

    # Within-dyad residual covariance
    (0 + .is_female + .is_male | couple_id),
  dispformula = ~ 0,
  family = gaussian(),
  data = apim_distinguishable_data
)

summary(apim_distinguishable_model)

Workflow

1. Define dyads, distinguishability, and research questions
2. Prepare the data and validate the dyadic structure
3. Describe and visualize measures and observed relationships
4. Estimate the models
5. Verify model adequacy, distinguishability, and robustness
6. Visualize, interpret, and report results

Verify model adequacy, distinguishability, and robustness

1. Did estimation succeed?
2. Does the model reproduce the important data patterns?
3. Is distinguishability supported?
4. Are conclusions robust or driven by few dyads or consequential modeling choices? (sensitivity analyses)

Did estimation succeed / converge?

The model converged without warnings… a good first sign!

Residual diagnostics / reproduction of data patterns

Residual diagnostics are not as straightforward anymore

• Linear models: straightforward checks
• Mixed models: residuals arise at multiple levels
• Generalized models: expectations depend on the family and link
• Dyadic models: partners’ residuals are correlated

Solution: Compare the data with simulations from the fitted model.

DHARMa: one diagnostic scale across models

The DHARMa package (Hartig, 2026):

• Simulates outcomes from the fitted model
• Checks whether the observed residuals behave as expected under simulations from the fitted model.
• Adequate fit \(\rightarrow\) approximately uniform, pattern-free residuals

One function works across many model families!

Simulate the full dyadic residual structure

library(DHARMa)

apim_distinguishable_dharma_residuals <- DHARMa::simulateResiduals(
  fittedModel = apim_distinguishable_model,
  n = 10000,

  # Since our model uses dispformula = ~ 0, we cannot use DHARMa's default
  # conditional simulations. With unconditional, the random-effects block representing
  # our residuals is re-simulated.
  simulateREs = "unconditional",

  # Male and female residuals are correlated. DHARMa estimates this correlation
  # from the simulations and rotates the residuals to "remove" it.
  rotation = 'estimated',
  seed = 123
)

Inspect DHARMa residuals by role

apim_distinguishable_female_dharma_residuals <- DHARMa::recalculateResiduals(
  apim_distinguishable_dharma_residuals,

  # Recalculate from the same simulations for one role.
  # One observation per couple: no rotation needed.
  sel = ~ .is_female == 1,
  rotation = NULL,
  seed = 123
)

apim_distinguishable_male_dharma_residuals <- DHARMa::recalculateResiduals(
  apim_distinguishable_dharma_residuals,
  sel = ~ .is_male == 1,
  rotation = NULL,
  seed = 123
)

Transform—or change the outcome family?

• Transformation may help for a nonnegative, right-skewed outcome
  • Try sqrt() or log1p() (handles zeros)
  • Refit and rerun DHARMa; the histogram alone is not sufficient
• A non-Gaussian family may be preferable
  • Particularly for counts, proportions, bounded, or zero-heavy outcomes
  • Family and link choices require outcome-specific care
• Workshop scope: generalized models are not covered here

Test distinguishability

• Fit model with exchangeability constraints

  • Partial constraints can be very informative!

• Compare models

• Pool only if theoretically defensible, aligned with RQ and theory, and compatible with the data

  • Depending on goals you can also report the distinguishable model and report results from this analysis

However, non-significant comparisons are not proof that roles are identical.

Gistelinck et al. (2018)

Estimate and compare:

dyadMLM::compare_nested_models(
  apim_distinguishable_model, apim_exchangeable_comparison_model
)

Likelihood-ratio test for nested models fitted to equivalent data
Assumes mathematical nesting and an appropriate chi-squared reference distribution.

                                   Df    AIC    BIC  logLik deviance  Chisq
apim_exchangeable_comparison_model  5 604.25 615.64 -297.13   594.25       
apim_distinguishable_model          9 605.44 625.93 -293.72   587.44 6.8155
                                   Chi Df Pr(>Chisq)
apim_exchangeable_comparison_model                  
apim_distinguishable_model              4      0.146

Conclusion (5% level): The likelihood-ratio test finds no clear improvement from `apim_exchangeable_comparison_model` to `apim_distinguishable_model` (p = 0.146). This does not establish equal fit.

Test robustness - are conclusions driven by few dyads?

plot_actor_association <- function(data, model, role, title,
                                   n_mahalanobis = 3) {
  xy <- cbind(data$.efficacy_gmc_actor, data$total_mvpa)
  labels <- data |>
    dplyr::mutate(.distance = mahalanobis(xy, colMeans(xy), cov(xy))) |>
    dplyr::slice_max(.distance, n = n_mahalanobis)
  coefficients <- glmmTMB::fixef(model)$cond

  ggplot(data, aes(.efficacy_gmc_actor, total_mvpa)) +
    geom_point(color = "steelblue") +
    geom_abline(intercept = coefficients[role],
                slope = coefficients[paste0(role, ":.efficacy_gmc_actor")],
                color = "navy", linewidth = 16 / 11) +
    geom_text(
      data = labels, aes(label = couple_id),
      vjust = -0.6, check_overlap = TRUE, size = 4
    ) +
    labs(
      title = title,
      x = "Actor efficacy (grand-mean centered)",
      y = "Daily MVPA (minutes)"
    ) +
    scale_y_continuous(expand = expansion(mult = c(0.05, 0.12))) +
    see::theme_modern(base_size = 16) +
    theme(plot.title.position = "panel",
          plot.title = element_text(hjust = 0.5))
}

Test robustness - are conclusions driven by few dyads?

n_mahalanobis <- 4

see::plots(
  plot_actor_association(
    dplyr::filter(apim_distinguishable_data, gender == "female"),
    apim_distinguishable_model, ".is_female", "Women",
    n_mahalanobis = n_mahalanobis
  ),
  plot_actor_association(
    dplyr::filter(apim_distinguishable_data, gender == "male"),
    apim_distinguishable_model, ".is_male", "Men",
    n_mahalanobis = n_mahalanobis
  )
) +
  patchwork::plot_layout(
    axes = "collect",
    axis_titles = "collect"
  )

Omitting Couple 28 changes the associations

apim_without_couple_28 <- update(
  apim_distinguishable_model,
  data = filter(apim_distinguishable_data, couple_id != 28)
)

actor_sensitivity <- parameters::compare_parameters(
  `Full data` = apim_distinguishable_model,
  `Without Couple 28` = apim_without_couple_28,
  keep = "\\.efficacy_gmc_actor"
)

Other sensitivity analyses you could do

Fit a small set of (prespecified) models to test:

• coding decisions
• different plausible outcome families or transformations
• missing-data handling
• covariates with a substantive confounding or precision rationale

Workflow

1. Define dyads, distinguishability, and research questions
2. Prepare the data and validate the dyadic structure
3. Describe and visualize measures and observed relationships
4. Estimate the models
5. Verify model adequacy, distinguishability, and robustness
6. Visualize, interpret, and report results

Reporting results and answering the research question

library(parameters)

parameters::model_parameters(apim_distinguishable_model, effects = 'fixed') |>
  insight::print_md(footer = "")

Reporting results and answering the research question

Random effects in our model represent the residual structure

parameters::model_parameters(apim_distinguishable_model, effects = 'random') |>
  insight::print_md(footer = "")

Visualize results

parameters::model_parameters(apim_distinguishable_model, effects = "fixed") |>
  plot()

What should be reported?

• Analysis sample: dyads and observations, exclusions, aggregation, and missing-data handling
• Model specification: outcome family and link, fixed effects, dyadic covariance structure, and equality constraints
• Estimation and model decision: estimator, software and version, convergence, boundary estimates, distinguishability comparison, and retained model
• Results: focal estimates with 95% confidence intervals in substantive units and the fitted variance-covariance parameters
• Checks and interpretation: residual diagnostics, influential dyads, prespecified sensitivity analyses, material changes, and causal limits

Workflow

1. Define dyads, distinguishability, and research questions
2. Prepare the data and validate the dyadic structure
3. Describe and visualize measures and observed relationships
4. Estimate the models
5. Verify model adequacy, distinguishability, and robustness
6. Visualize, interpret, and report results

Define dyads, distinguishability, and research questions

Same data and RQ, now with exchangeability constraints for illustration.

Workflow

1. Define dyads, distinguishability, and research questions
2. Prepare the data and validate the dyadic structure
3. Describe and visualize measures and observed relationships
4. Estimate the models
5. Verify model adequacy, distinguishability, and robustness
6. Visualize, interpret, and report results

Create the required format and validate the data

You can use dyadMLM::prepare_dyad_data() to prepare and validate the data.

For exchangeable dyads, we can just omit the role argument.

apim_exchangeable_data <- dyadMLM::prepare_dyad_data(
  data = raw_dyad_data,
  dyad = couple_id,
  member = person_id,
  predictors = efficacy,
  model_types = 'apim',
  add_apim_gmc_predictors = TRUE,
  seed = 123
)

print(apim_exchangeable_data)

Create the required format and validate the data

Looking at the relevant columns:

head(apim_exchangeable_data) |>
  dplyr::select(
    couple_id, person_id, .member_contrast_arbitrary, .efficacy_gmc,
    .efficacy_gmc_actor, .efficacy_gmc_partner, total_mvpa
  ) |>
  slide_table()

Workflow

1. Define dyads, distinguishability, and research questions
2. Prepare the data and validate the dyadic structure
3. Describe and visualize measures and observed relationships
4. Estimate the models
5. Verify model adequacy, distinguishability, and robustness
6. Visualize, interpret, and report results

Describe focal variables (pooled)

No wide-view of the data is needed for exchangeable dyads. We can describe the data directly to obtain pooled estimates:

apim_exchangeable_data |>
  dplyr::select(efficacy, total_mvpa) |>
  report::report_table() |>
  slide_table()

Visualize focal variables

plot_pooled_histogram <- function(data, var) {
  var_label <- rlang::as_label(rlang::enquo(var))

  ggplot(
    data,
    aes(x = {{ var }}, y = after_stat(density))
  ) +
    geom_histogram(
      bins = 15, boundary = 0,
      fill = "steelblue", color = "white", alpha = 0.75
    ) +
    geom_density(
      bounds = c(0, Inf), color = "navy", linewidth = 1.1
    ) +
    see::theme_modern(base_size = 16)
}

Visualize focal variables

see::plots(
  plot_pooled_histogram(apim_exchangeable_data, efficacy),
  plot_pooled_histogram(apim_exchangeable_data, total_mvpa)
)

Observed correlations

For exchangeable dyads, ICCs summarize partner similarity:

library(wbCorr)

apim_exchangeable_correlations <- wbCorr(
  apim_exchangeable_data[,c('efficacy','total_mvpa')],
  cluster = apim_exchangeable_data$couple_id
)

summary(apim_exchangeable_correlations, 'wb')

$merged_wb
           efficacy total_mvpa
efficacy     [0.60]     0.44**
total_mvpa  0.56***     [0.79]

$note
[1] "***p < 0.001, **p < 0.01, *p < 0.05"

Diagonal ICCs: expected correlation between partners on each variable

Off-diagonals: top right = within-couple correlation; bottom left = between-couple correlation

Visualize partner similarity

partner_similarity_data <- apim_exchangeable_data |>
  dplyr::group_by(couple_id) |>
  dplyr::mutate(
    efficacy_partner = rev(efficacy),
    total_mvpa_partner = rev(total_mvpa)
  ) |>
  dplyr::ungroup()

see::plots(
  plot_partner_similarity(
    partner_similarity_data,
    efficacy, efficacy_partner, couple_id, "Efficacy",
    n_mahalanobis = 4
  ) +
    labs(x = "Focal member", y = "Partner"),
  plot_partner_similarity(
    partner_similarity_data,
    total_mvpa, total_mvpa_partner, couple_id, "Daily MVPA",
    n_mahalanobis = 4
  ) +
    labs(x = "Focal member", y = "Partner")
)

Workflow

1. Define dyads, distinguishability, and research questions
2. Prepare the data and validate the dyadic structure
3. Describe and visualize measures and observed relationships
4. Estimate the models
5. Verify model adequacy, distinguishability, and robustness
6. Visualize, interpret, and report results

Exchangeable APIM

In MLM we constrain / pool by omitting terms

One single equation:

For each dyad \(i\), each member \(j\) has their own predictor \(X^{(\mathrm{A})}_{j,i}\) and their partner’s predictor \(X^{(\mathrm{P})}_{j,i}\):

\[ Y_{j,i} = b_{0} + \color{#1769AA}{a \times X^{(\mathrm{A})}_{j,i}} + \color{#7B3FA1}{p \times X^{(\mathrm{P})}_{j,i}} + \epsilon_{j,i} \]

- Actor effects - Partner effects

In long format, each equation applies to one member row of the dyad.

The tricky part: residuals need to remain correlated within dyads

• The two member rows from the same dyad retain a joint residual structure

• Both members now have the same residual variance (pooled)

• Position 1 and 2 within the dyad are arbitrary

\[ \operatorname{Cov} \begin{pmatrix} \epsilon_{1i} \\ \epsilon_{2i} \end{pmatrix} = \boldsymbol{\Sigma}_{\epsilon} = \begin{bmatrix} \sigma_{\epsilon}^{2} & \rho_{\epsilon_1\epsilon_2}\sigma_{\epsilon}^{2} \\ \rho_{\epsilon_1\epsilon_2}\sigma_{\epsilon}^{2} & \sigma_{\epsilon}^{2} \end{bmatrix} \]

Fitting the exchangeable APIM in glmmTMB

apim_exchangeable_model <- glmmTMB(
  total_mvpa ~
    # Pooled intercept
    1 +

    # Pooled actor effect
    .efficacy_gmc_actor +

    # Pooled partner effect
    .efficacy_gmc_partner +

    #------ Within-dyad residual covariance (Sum/Deviation) ------

    # Residuals of partners might be similar (positively correlated)
    (1 | couple_id) +

    # Residuals of partners might move opposite (negatively correlated)
    (0 + .member_contrast_arbitrary | couple_id),

  dispformula = ~ 0,
  family = gaussian(),
  data = apim_exchangeable_data
)

summary(apim_exchangeable_model)

The mean-deviation parameterization

• The DIM is equivalent to the APIM

• The DIM represents member residuals using two uncorrelated components:

  • Mean block (1 | couple_id): \(r_{\mathrm{M}i}=(\epsilon_{1i}+\epsilon_{2i})/2\); both members move together

  • Deviation block (0 + .member_contrast_arbitrary | couple_id): \(r_{\mathrm{D}i}=(\epsilon_{1i}-\epsilon_{2i})/2\); members move in opposite directions

Using this parametrization allows us to represent random effects and later in ILD even random slopes with the correct constraints.

We can later backtransform to the familiar person-level structure.

Workflow

1. Define dyads, distinguishability, and research questions
2. Prepare the data and validate the dyadic structure
3. Describe and visualize measures and observed relationships
4. Estimate the models
5. Verify model adequacy, distinguishability, and robustness
6. Visualize, interpret, and report results

Verify model adequacy, distinguishability, and robustness

Same as before:

1. Did estimation succeed?
2. Does the model reproduce the important data patterns?
3. Does a distinguishable model fit clearly better?
4. Are conclusions robust or driven by few dyads or consequential modeling choices? (sensitivity analyses)

Workflow

1. Define dyads, distinguishability, and research questions
2. Prepare the data and validate the dyadic structure
3. Describe and visualize measures and observed relationships
4. Estimate the models
5. Verify model adequacy, distinguishability, and robustness
6. Visualize, interpret, and report results

Reporting results and answering the research questions

parameters::model_parameters(apim_exchangeable_model, effects = 'fixed') |>
  insight::print_md(footer = "")

Reporting results and answering the research questions

parameters::model_parameters(apim_exchangeable_model, effects = 'random') |>
  insight::print_md(footer = "")

Reporting results and answering the research questions

No need to do so by hand. Rotating back with dyadMLM:

dyadMLM::recover_exchangeable_covariance(apim_exchangeable_model)

Recovered exchangeable member-level covariance

Pair `pair_1`
Shared:     us(1 | couple_id)
Difference: us(0 + .member_contrast_arbitrary | couple_id)

Variance-covariance:
                       1       2      
1 member1: (Intercept) 343.216 259.269
2 member2: (Intercept) 259.269 343.216

Standard deviations and correlations:
                       1      2     
1 member1: (Intercept) 18.526 0.755 
2 member2: (Intercept) 0.755  18.526

Visualize results

parameters::model_parameters(apim_exchangeable_model, effects = 'fixed') |>
  plot(show_intercept = TRUE)

Intensive longitudinal APIM - Workflow

We can follow the same workflow. Additional possible ILD considerations:

• Check compliance / missingness
• Inspect and visualize paired trajectories and patterns
• Separate within- and between-person variation
  • Inspect ICCs and multilevel reliability where applicable
  • Model dyads on multiple levels
• Check and account for residual temporal dependence

Workflow

1. Define dyads, distinguishability, and research questions
2. Prepare the data and validate the dyadic structure
3. Describe and visualize measures and observed relationships
4. Estimate the models
5. Verify model adequacy, distinguishability, and robustness
6. Visualize, interpret, and report results

The Dataset

Same as before:

• 36 physically inactive romantic couples
• Both partners intended to become more physically active
• Smartphone-based planning interventions and JITAIs
• Daily diaries over 55 days

Not aggregated anymore

Measurement - predictor

Perceived behavioral control (PBC) over next-day physical activity

Rough English translation of the original German item
“Being physically active tomorrow would/will be for me …”

Adapted from Ajzen (1991)

0 = impossible—crucial prerequisites are not met
3 = possible with some effort—several prerequisites are not met
6 = easily possible—all prerequisites are met

• Lagged predictor: perceived behavioral control on day \(t-1\), MVPA on day \(t\)

Measurement - outcome

Device-based MVPA

Accelerometer-derived measure
Daily minutes of moderate-to-vigorous physical activity

• Days with less than 10 hours of awake wear time \(\rightarrow\) NA
• Analyze \(\log(\text{MVPA})\) to reduce right skew

Workflow

1. Define dyads, distinguishability, and research questions
2. Prepare the data and validate the dyadic structure
3. Describe and visualize measures and observed relationships
4. Estimate the models
5. Verify model adequacy, distinguishability, and robustness
6. Visualize, interpret, and report results

Load the daily measures

We set MVPA to NA on days with less than 10 hours of weartime, center time at the middle of the study in full-study units, and select the needed variables.

raw_ild_dyad_data <- readRDS("dyadic-data.rds") |>
  dplyr::mutate(
    device_based_mvpa = dplyr::if_else(
      awake_wear_minutes >= 600,
      device_based_mvpa,
      NA
    ),
    log_device_based_mvpa = log(device_based_mvpa),
    diaryday_c = (diaryday - 27) / 54
  )|>
  dplyr::select(
    couple_id, person_id, diaryday, diaryday_c, gender,
    pbc = perceived_behavioral_control,
    device_based_mvpa, log_device_based_mvpa, awake_wear_minutes
  ) |>
  dplyr::arrange(couple_id, diaryday, person_id)

slide_table(head(raw_ild_dyad_data))

Create the required format and validate the data

For ILD we decompose the predictor(s):

• Between-person component: Person mean centered on the mean of person means (stable)
  • \(\bar X_{ri} - \bar X\)
• Within-person component: Daily deviations from the person mean
  • \(X_{ri,t} - \bar X_{ri}\)

In this example we then lag the within-person component (\(X_{ri,t-1}\))

Create the required format and validate the data Manual code →

The dyadMLM package can handle centering alongside the construction and validation of the dyadic dataset:

ild_apim_data <- dyadMLM::prepare_dyad_data(
  data = raw_ild_dyad_data,
  dyad = couple_id,
  member = person_id,
  role = gender,
  time = diaryday,
  predictors = pbc,
  lag1_predictors = pbc,
  model_types = "apim",
  temporal_decomposition = "2l",
  include_arbitrary_member_contrast = TRUE,
  seed = 123
)

ild_apim_data <- ild_apim_data |>
  dplyr::mutate(
    diaryday_f = factor(diaryday, levels = 0:54)
  )

print(ild_apim_data)

Workflow

1. Define dyads, distinguishability, and research questions
2. Prepare the data and validate the dyadic structure
3. Describe and visualize measures and observed relationships
4. Estimate the models
5. Verify model adequacy, distinguishability, and robustness
6. Visualize, interpret, and report results

Report measures

Among other things you may check and report:

• Missingness and general summary statistics
• ICCs
• Within- and between- person correlations (e.g., with wbCorr)
• Visualize outcome trajectories over time

Plot paired raw MVPA trajectories

raw_ild_dyad_data |>
  dplyr::filter(couple_id %in% c(2, 7, 17, 26)) |>
  ggplot(aes(diaryday, device_based_mvpa, color = gender,
             group = person_id)) +
  geom_line(linewidth = 0.65, na.rm = TRUE) +
  geom_point(size = 1.1, na.rm = TRUE) +
  facet_wrap(~ couple_id, ncol = 2) +
  scale_color_manual(values = c(female = "#004D40", male = "#00A6B2")) +
  scale_y_continuous(breaks = seq(0, 300, by = 100)) +
  labs(x = "Diary day", y = "MVPA (minutes)", color = NULL) +
  see::theme_modern(base_size = 16) +
  theme(legend.position = "top")

Workflow

1. Define dyads, distinguishability, and research questions
2. Prepare the data and validate the dyadic structure
3. Describe and visualize measures and observed relationships
4. Estimate the models
5. Verify model adequacy, distinguishability, and robustness
6. Visualize, interpret, and report results

Fitting the ILD APIM: Two member-specific equations

For dyad \(i\) at occasion \(t\), superscripts \(\mathrm{W}\) and \(\mathrm{B}\) denote within-person and between-person components.

\[ \begin{aligned} Y_{\mathrm{M},i,t} ={}& b_{0,\mathrm{M}} + \color{#1769AA}{ a^{\mathrm{W}}_{\mathrm{M}}X^{(\mathrm{A,W})}_{\mathrm{M},i,t-1} + a^{\mathrm{B}}_{\mathrm{M}}X^{(\mathrm{A,B})}_{\mathrm{M},i} } + \color{#7B3FA1}{ p^{\mathrm{W}}_{\mathrm{M}}X^{(\mathrm{P,W})}_{\mathrm{M},i,t-1} + p^{\mathrm{B}}_{\mathrm{M}}X^{(\mathrm{P,B})}_{\mathrm{M},i} } + u_{0,\mathrm{M}i} + \epsilon_{\mathrm{M},i,t} \\[0.6em] Y_{\mathrm{F},i,t} ={}& b_{0,\mathrm{F}} + \color{#1769AA}{ a^{\mathrm{W}}_{\mathrm{F}}X^{(\mathrm{A,W})}_{\mathrm{F},i,t-1} + a^{\mathrm{B}}_{\mathrm{F}}X^{(\mathrm{A,B})}_{\mathrm{F},i} } + \color{#7B3FA1}{ p^{\mathrm{W}}_{\mathrm{F}}X^{(\mathrm{P,W})}_{\mathrm{F},i,t-1} + p^{\mathrm{B}}_{\mathrm{F}}X^{(\mathrm{P,B})}_{\mathrm{F},i} } + u_{0,\mathrm{F}i} + \epsilon_{\mathrm{F},i,t}. \end{aligned} \]

Actor effects · Partner effects

Repeated measures separate two covariance levels

Fitting the distinguishable ILD APIM with glmmTMB

ild_apim_distinguishable_model <- glmmTMB::glmmTMB(
  log_device_based_mvpa ~

    # Role-specific fixed intercepts
    0 + .is_female + .is_male +

    # Role-specific fixed time-slope
    .is_female:diaryday_c + .is_male:diaryday_c +

    # --------------- WITHIN PERSON APIM ---------------
    # Role-specific fixed actor- and partner effects

    .is_female:.pbc_cwp_actor_lag1 +
    .is_male:.pbc_cwp_actor_lag1 +

    .is_female:.pbc_cwp_partner_lag1 +
    .is_male:.pbc_cwp_partner_lag1 +

    # --------------- BETWEEN PERSON APIM ---------------
    # Role-specific fixed actor-and partner effects

    .is_female:.pbc_cbp_actor +
    .is_male:.pbc_cbp_actor +

    .is_female:.pbc_cbp_partner +
    .is_male:.pbc_cbp_partner +

    # Stable dyad/member random intercepts (correlated)
    (0 + .is_female + .is_male | couple_id) +

    # Same-occasion Gaussian residual covariance
    (0 + .is_female + .is_male | couple_id:diaryday),
  dispformula = ~ 0,
  family = gaussian(),
  data = ild_apim_data
)

summary(ild_apim_distinguishable_model)

Workflow

1. Define dyads, distinguishability, and research questions
2. Prepare the data and validate the dyadic structure
3. Describe and visualize measures and observed relationships
4. Estimate the models
5. Verify model adequacy, distinguishability, and robustness
6. Visualize, interpret, and report results

Testing for residual autocorrelation

ild_apim_no_ar_dharma <- simulate_ild_dharma(
  model = ild_apim_distinguishable_model,
  dyad = ild_apim_data$couple_id,
  member = ild_apim_data$person_id,
  time = ild_apim_data$diaryday,
  n = 1000,
  seed = 123
)

diagnostic <- test_ild_lag1(
  ild_apim_no_ar_dharma,
  plot = TRUE
)

Add role-specific residual AR(1) processes

• Separate AR(1) series for the female and male partner
• Role-specific residual variability and persistence \(\rho\)
• Same-day partner covariance retained separately

ild_apim_distinguishable_ar1_model <- update(
  ild_apim_distinguishable_model,
  . ~ . +
    ar1(0 + .is_female:diaryday_f | couple_id) +
    ar1(0 + .is_male:diaryday_f | couple_id)
)

glmmTMB::VarCorr(ild_apim_distinguishable_ar1_model)

Simulate residuals from the AR(1) model

Does the fitted AR(1) reproduce lag-1 dependence?

ild_apim_ar1_lag1_test <- test_ild_lag1(
  ild_apim_dharma_residuals,
  plot = TRUE
)

Fit the exchangeable model with pooled AR(1)

ild_apim_exchangeable_ar1_model <- glmmTMB::glmmTMB(
  log_device_based_mvpa ~
    1 +
    diaryday_c +
    .pbc_cwp_actor_lag1 +
    .pbc_cwp_partner_lag1 +
    .pbc_cbp_actor +
    .pbc_cbp_partner +

    # Couple-level variance-covariance structure
    (1 | couple_id) +
    (0 + .member_contrast_arbitrary | couple_id) +

    # Occasion-level residual variance-covariance structure
    (1 | couple_id:diaryday) +
    (0 + .member_contrast_arbitrary | couple_id:diaryday) +
    ar1(0 + diaryday_f | couple_id:person_id),
  dispformula = ~ 0,
  family = gaussian(),
  data = ild_apim_data
)

summary(ild_apim_exchangeable_ar1_model)

Backtransform to member-level variance-covariance structure

dyadMLM::recover_exchangeable_covariance(ild_apim_exchangeable_ar1_model)

Recovered exchangeable member-level covariances (2 block pairs)

Pair `pair_1`
Shared:     us(1 | couple_id)
Difference: us(0 + .member_contrast_arbitrary | couple_id)

Variance-covariance:
                       1     2    
1 member1: (Intercept) 0.075 0.012
2 member2: (Intercept) 0.012 0.075

Standard deviations and correlations:
                       1     2    
1 member1: (Intercept) 0.274 0.160
2 member2: (Intercept) 0.160 0.274

Pair `pair_2`
Shared:     us(1 | couple_id:diaryday)
Difference: us(0 + .member_contrast_arbitrary | couple_id:diaryday)

Variance-covariance:
                       1     2    
1 member1: (Intercept) 0.134 0.051
2 member2: (Intercept) 0.051 0.134

Standard deviations and correlations:
                       1     2    
1 member1: (Intercept) 0.366 0.380
2 member2: (Intercept) 0.380 0.366

Compare distinguishable and exchangeable models

This is a global test, partial constraints can be informative too!

dyadMLM::compare_nested_models(
  ild_apim_distinguishable_ar1_model,
  ild_apim_exchangeable_ar1_model
)

Likelihood-ratio test for nested models fitted to equivalent data
Assumes mathematical nesting and an appropriate chi-squared reference distribution.

                                   Df    AIC    BIC  logLik deviance  Chisq
ild_apim_exchangeable_ar1_model    12 3052.4 3123.8 -1514.2   3028.4       
ild_apim_distinguishable_ar1_model 22 3040.9 3171.8 -1498.5   2996.9 31.487
                                   Chi Df Pr(>Chisq)    
ild_apim_exchangeable_ar1_model                         
ild_apim_distinguishable_ar1_model     10  0.0004873 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Conclusion (5% level): The likelihood-ratio test provides evidence that `ild_apim_distinguishable_ar1_model` fits better than `ild_apim_exchangeable_ar1_model` (p < 0.001).

Workflow

1. Define dyads, distinguishability, and research questions
2. Prepare the data and validate the dyadic structure
3. Describe and visualize measures and observed relationships
4. Estimate the models
5. Verify model adequacy, distinguishability, and robustness
6. Visualize, interpret, and report results

Report the retained distinguishable model

We exponentiate estimates from the log-transformed outcome to obtain geometric mean ratios (GMRs; representing multiplicative differences).

parameters::model_parameters(
  ild_apim_distinguishable_ar1_model,
  effects = 'fixed',
  exponentiate = TRUE
  ) |>
  insight::print_md(footer = "")

Report the retained distinguishable model

see::plots(
  parameters::model_parameters(
    ild_apim_distinguishable_ar1_model,
    effects = "fixed",
    exponentiate = TRUE,
    keep = "_cwp_"
  ) |> plot(),
  parameters::model_parameters(
    ild_apim_distinguishable_ar1_model,
    effects = "fixed",
    exponentiate = TRUE,
    keep = "_cbp_"
  ) |> plot(),
  n_columns = 2,
  tags = c("Within person", "Between persons")
) &
  ggplot2::labs(x = "GMR")

Report the retained distinguishable model

Within person: Men’s MVPA was 3.5% higher when their own previous-day PBC was one point above usual, and 2.5% higher when their partner’s was.

Between persons: There was no clear evidence of actor or partner associations.

Many studies lack complete reporting

Systematic review of 73 longitudinal APIM studies published from 2008–2025 (Mohammadzadeh Yazd et al., 2026).

Make every dyadic analysis decision visible

1. Preregister actor and partner hypotheses and key analysis decisions.
2. Explain missing data: describe missingness and the handling strategy.
3. Specify the full model: roles and distinguishability, actor and partner paths, constraints, and covariance structure.
4. Report reliability where applicable: within- and between-level omega for multilevel item scales, e.g., multilevelTools::omegaSEM() (Geldhof et al., 2014; Wiley, 2025).

Mohammadzadeh Yazd et al. (2026)