Distinguishable and Exchangeable Dyads: Bayesian Multilevel Modelling

Cross-Sectional and Intensive Longitudinal APIM and DIM

Pascal M. Küng

2026-08-26

Overview

  • Distinguishable Dyads
    • Cross-Sectional APIM
    • Longitudinal APIM
  • Exchangeable Dyads
    • Cross-Sectional APIM
    • Cross-Sectional DIM
    • Equivalence between APIM and DIM
    • Longitudinal DIM
    • Longitudinal APIM

Quick Access

For readers who want the model code and reporting tables without the full tutorial:

Preparing Dyadic Data with dyadMLM

dyadMLM prepares dyadic data for multilevel models (Küng, 2026). It:

  • validates long-format dyadic data
  • supports cross-sectional and intensive longitudinal designs
  • creates model-ready variables for APIM, DIM, and DSM specifications
  • handles distinguishable, exchangeable, and mixed dyad compositions

Install the current R-universe version, using CRAN for dependencies:

install.packages("dyadMLM", repos = c(
  "https://pascal-kueng.r-universe.dev",
  "https://cloud.r-project.org"
))

Getting started · Function reference · Source code

Loading Libraries and setting CmdStan Backend

Some functions used in this presentation are sourced from the files available in the folder 00_R_Functions.

library(tidyverse)
library(dyadMLM)
library(brms)
library(glmmTMB)
library(easystats)
library(DiagrammeR)
library(DHARMa)
library(MASS)
library(purrr)
source(file.path('00_R_Functions', 'ReportModels.R'))
source(file.path('00_R_Functions', 'PrettyTables.R'))

print_apim_preview <- function(df) {
  preview_cols <- c(
    "userID", "coupleID", "gender", ".is_male", ".is_female",
    ".member_contrast_arbitrary", "satisfaction",
    ".communication_gmc_actor", ".communication_gmc_partner"
  )
  df <- df |>
    dplyr::select(dplyr::any_of(preview_cols))
  print_df(head(df))
}

options(
  brms.backend = 'cmdstanr',
  brms.file_refit = 'on_change'
)

Distinguishable Dyads

Distinguishable Dyads - Cross-Sectional APIM

Distinguishable Dyads - Data: Simulated Dyads

print_df(head(df))
userID coupleID gender communication satisfaction
1_1 1 female 4.851107 4.841332
1_2 1 male 6.004533 5.577165
2_1 2 female 5.881321 7.248741
2_2 2 male 6.433723 6.960293
3_1 3 female 4.283971 6.928699
3_2 3 male 5.516060 6.077688

Distinguishable Dyads - Cross-Sectional APIM

Distinguishable APIM (e.g., Kenny et al., 2006; Kenny & Cook, 1999; Kenny & Kashy, 2011)

Distinguishable Dyads - Preparing Data

df_apim <- prepare_dyad_data(
  data = df,
  dyad = coupleID,
  member = userID,
  role = gender,
  predictors = communication,
  model_types = "apim",
  keep_compositions = "female-male",
  include_arbitrary_member_contrast = TRUE,
  add_apim_gmc_predictors = TRUE,
  seed = 123
)

print_apim_preview(df_apim)

Distinguishable Dyads - Preparing Data

userID coupleID gender .is_male .is_female .member_contrast_arbitrary satisfaction .communication_gmc_actor .communication_gmc_partner
1_1 1 female 0 1 -1 4.841332 -0.1553604 0.9980657
1_2 1 male 1 0 1 5.577165 0.9980657 -0.1553604
2_1 2 female 0 1 -1 7.248741 0.8748537 1.4272563
2_2 2 male 1 0 1 6.960293 1.4272563 0.8748537
3_1 3 female 0 1 -1 6.928699 -0.7224960 0.5095933
3_2 3 male 1 0 1 6.077688 0.5095933 -0.7224960

For the likelihood-based constraint comparison, the glmmTMB models below use .communication_gmc_actor and .communication_gmc_partner. The Bayesian models retain the raw actor and partner scores because centering also changes the intercept’s prior interpretation.

Distinguishable Dyads - Fitting the Model in glmmTMB

  • glmmTMB models residual variance/dispersion via dispformula, but residual covariance is handled through random-effect covariance structures.

  • For Gaussian mixed models, dispformula = ~ 0 fixes the residual variance to a tiny value. This pushes the Gaussian residual variance into the random-effect block.

  • In this Gaussian dyadic model, the us(0 + .is_male + .is_female | coupleID) block therefore represents the 2x2 residual variance-covariance structure between partners.

  • For non-Gaussian models, this is not an ordinary residual covariance on the response scale. The same random-effect trick adds latent Gaussian heterogeneity/correlation on the linear predictor scale.

  • Do not use dispformula = ~ 0 as a general non-Gaussian residual-covariance trick. For families with an estimable dispersion parameter, dispformula models that family-specific parameter via a log link; for families without one, such as Poisson or binomial, it is ignored. Partner-specific dispersion, e.g. dispformula = ~ 0 + .is_male + .is_female, is only appropriate when the family has a dispersion parameter and this is substantively intended.

Distinguishable Dyads - Fitting the Model in glmmTMB

model_dist_apim_glmmtmb <- glmmTMB(
  satisfaction ~

    # Remove global intercept, introduce male and female intercepts.
    0 + .is_male + .is_female +

    # Actor effect for male and female
    .communication_gmc_actor:.is_male +
    .communication_gmc_actor:.is_female +

    # Partner effect for male and female
    .communication_gmc_partner:.is_male +
    .communication_gmc_partner:.is_female +

    # Modelling residual non-independence with an unstructured 2x2 covariance.
    us(0 + .is_male + .is_female | coupleID)

    # Gaussian response:
    # dispformula = ~ 0 fixes residual variance to a tiny value,
    # pushing the Gaussian residual variance into this random-effect block.
    , dispformula = ~ 0

    # Non-Gaussian responses:
    # dispformula is family-specific dispersion, or ignored for families
    # without a dispersion parameter. Do not use ~0 as a residual trick.
    # If dispersion exists and should differ by partner, use e.g.:
    # dispformula = ~ 0 + .is_male + .is_female

  , data = df_apim
  , family = gaussian
)

Distinguishable Dyads - Fitting the Model in glmmTMB

Note that ‘unstructured’ is the default in glmmTMB, so instead of

us(0 + .is_male + .is_female | coupleID)

we could also remove us and simply write

(0 + .is_male + .is_female | coupleID).

An unstructured covariance matrix gives .is_male and .is_female their own variance and freely estimates their correlation. With more than two members (e.g., triads), it would estimate a separate correlation for every pair.

For two partners, the covariance matrix is:

\[ \boldsymbol{\Sigma} = \begin{pmatrix} \sigma^2_{female} & \rho\,\sigma_{female}\sigma_{male} \\ \rho\,\sigma_{female}\sigma_{male} & \sigma^2_{male} \end{pmatrix}. \]

Distinguishable Dyads - Fitting the Model in glmmTMB

 Family: gaussian  ( identity )
Formula:          satisfaction ~ 0 + .is_male + .is_female + .communication_gmc_actor:.is_male +  
    .communication_gmc_actor:.is_female + .communication_gmc_partner:.is_male +  
    .communication_gmc_partner:.is_female + us(0 + .is_male +      .is_female | coupleID)
Dispersion:                    ~0
Data: df_apim

      AIC       BIC    logLik -2*log(L)  df.resid 
   3776.7    3821.7   -1879.3    3758.7      1091 

Random effects:

Conditional model:
 Groups   Name       Variance Std.Dev. Corr 
 coupleID .is_male   2.559    1.600         
          .is_female 1.682    1.297    0.51 
Number of obs: 1100, groups:  coupleID, 550

Conditional model:
                                      Estimate Std. Error z value Pr(>|z|)    
.is_male                               4.60635    0.06821   67.54  < 2e-16 ***
.is_female                             5.61086    0.05530  101.47  < 2e-16 ***
.is_male:.communication_gmc_actor      1.84363    0.04931   37.39  < 2e-16 ***
.is_female:.communication_gmc_actor    1.59247    0.04112   38.72  < 2e-16 ***
.is_male:.communication_gmc_partner    0.16765    0.05072    3.31 0.000949 ***
.is_female:.communication_gmc_partner  0.31084    0.03997    7.78 7.49e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

glmmTMB: Diagnostics with Near-Zero Gaussian Sigma

With dispformula = ~ 0, glmmTMB fixes the Gaussian observation-level residual SD near zero because the dyadic variance and covariance are represented by the random-effect block (Brooks et al., 2026). DHARMa’s default conditional simulation therefore contains almost no observation-level variation. For a marginal check, resimulate the random effects and interpret the result for that target (Hartig, 2026):

dharma_glmmtmb <- DHARMa::simulateResiduals(
  model_dist_apim_glmmtmb,
  n = 1000,
  simulateREs = "unconditional"
)

Distinguishable Dyads - Fitting the Model in brms

In brms, we can use brms-provided residual covariance structures directly. We do not have to set residuals internally to zero.

formula_dist_apim_b <- bf(
  satisfaction ~

    # Remove global intercept, introduce male and female intercepts.
    0 + .is_male + .is_female +

    # Actor effect for male and female
    .communication_actor:.is_male +
    .communication_actor:.is_female +

    # Partner effect for male and female
    .communication_partner:.is_male +
    .communication_partner:.is_female +

    # Gender indexes the two within-couple residual positions.
    unstr(time = gender, gr = coupleID),

  # Allowing two distinct residual variances for males vs females:
  sigma = ~ 0 + .is_male + .is_female

)

priors_dist_apim_b <- c(
  # prior(normal(2, 3), class = "Intercept"),
  prior(normal(0, 5), class = "b"),
  prior(student_t(3, 0, 1.5), class = "b", dpar = "sigma"),
  prior(lkj(2), class = "cortime")
)

model_dist_apim_b <- brm(
  formula = formula_dist_apim_b,
  data = df_apim,
  family = gaussian(link = identity),
  prior = priors_dist_apim_b,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000,
  seed = 123,
  file = file.path('brms_cache', 'example1_dist_apim_b') # Cache the model
)

Distinguishable Dyads - Fitting the Model in brms

 Family: gaussian 
  Links: mu = identity; sigma = log 
Formula: satisfaction ~ 0 + .is_male + .is_female + .communication_actor:.is_male + .communication_actor:.is_female + .communication_partner:.is_male + .communication_partner:.is_female + unstr(time = gender, gr = coupleID) 
         sigma ~ 0 + .is_male + .is_female
   Data: df_apim (Number of observations: 1100) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Correlation Structures:
                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
cortime(female,male)     0.51      0.03     0.44     0.57 1.00     2895     2470

Regression Coefficients:
                                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
.is_male                             -5.44      0.28    -5.99    -4.88 1.00     2458     2617
.is_female                           -3.90      0.23    -4.37    -3.47 1.00     2401     2641
.is_male:.communication_actor         1.84      0.05     1.75     1.94 1.00     2678     2997
.is_female:.communication_actor       1.59      0.04     1.51     1.67 1.00     2371     2549
.is_male:.communication_partner       0.17      0.05     0.07     0.26 1.00     2555     2726
.is_female:.communication_partner     0.31      0.04     0.23     0.39 1.00     2398     2724
sigma_.is_male                        0.47      0.03     0.42     0.53 1.00     3555     2827
sigma_.is_female                      0.26      0.03     0.21     0.33 1.00     3317     2788

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Distinguishable Dyads - Poisson Covariance in brms

Poisson has no separately estimated observation-level residual SD. For distinguishable dyads, we therefore use correlated role-specific dyad effects. The non-Gaussian unstr() formulation instead has one common latent SD and would impose equal role variances (Bürkner, 2025). We simulate an actual count outcome rather than discretize continuous satisfaction.

Distinguishable Dyads - Fitting the Poisson Model in brms

formula_dist_apim_c <- bf(
  satisfaction_count ~

    # Remove global intercept, introduce male and female intercepts.
    0 + .is_male + .is_female +

    # Actor effect for male and female
    .communication_actor:.is_male +
    .communication_actor:.is_female +

    # Partner effect for male and female
    .communication_partner:.is_male +
    .communication_partner:.is_female +

    # Model within-dyad dependence with correlated Gaussian random effects.
    # Poisson has no separately estimated residual-dispersion parameter.

    (0 + .is_male + .is_female | coupleID)

    # sigma ~ 0 + .is_male + .is_female
)

priors_dist_apim_c <- c(
    prior(normal(0, 5), class = "b"),
    prior(student_t(3, 0, 1.5), class = "sd"),
    prior(lkj(2), class = "cor")
  )

model_dist_apim_c <- brm(
  formula = formula_dist_apim_c,
  data = df_apim_count,
  family = poisson(),
  prior = priors_dist_apim_c,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000,
  seed = 123,
  file = file.path('brms_cache', 'example1_dist_apim_c') # Cache the model
)

Poisson Model in brms: Results

 Family: poisson 
  Links: mu = log 
Formula: satisfaction_count ~ 0 + .is_male + .is_female + .communication_actor:.is_male + .communication_actor:.is_female + .communication_partner:.is_male + .communication_partner:.is_female + (0 + .is_male + .is_female | coupleID) 
   Data: df_apim_count (Number of observations: 1100) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Multilevel Hyperparameters:
~coupleID (Number of levels: 550) 
                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(.is_male)                 0.38      0.03     0.33     0.43 1.00     1771     2547
sd(.is_female)               0.27      0.03     0.22     0.32 1.00     1868     2480
cor(.is_male,.is_female)     0.30      0.11     0.08     0.52 1.00     1107     2003

Regression Coefficients:
                                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
.is_male                              0.47      0.11     0.27     0.68 1.00     2639     2753
.is_female                            0.64      0.09     0.46     0.81 1.00     3087     2953
.is_male:.communication_actor         0.21      0.02     0.17     0.25 1.00     2587     3041
.is_female:.communication_actor       0.17      0.02     0.14     0.20 1.00     2933     2656
.is_male:.communication_partner       0.01      0.02    -0.03     0.04 1.00     2519     2720
.is_female:.communication_partner     0.06      0.01     0.03     0.09 1.00     3308     2926

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Distinguishable Dyads - Fitting the Model in glmmTMB

model_glmmTMB_dist_apim_c <- glmmTMB(
  satisfaction_count ~

    # Remove global intercept, introduce male and female intercepts.
    0 + .is_male + .is_female +

    # Actor effect for male and female
    .communication_actor:.is_male +
    .communication_actor:.is_female +

    # Partner effect for male and female
    .communication_partner:.is_male +
    .communication_partner:.is_female +

    # Model within-dyad dependence with correlated Gaussian random effects.
    # Poisson has no separately estimated residual-dispersion parameter.

    (0 + .is_male + .is_female | coupleID)
  , data = df_apim_count
  , family = poisson()
)

Poisson Model in glmmTMB: Results

 Family: poisson  ( log )
Formula:          satisfaction_count ~ 0 + .is_male + .is_female + .communication_actor:.is_male +  
    .communication_actor:.is_female + .communication_partner:.is_male +  
    .communication_partner:.is_female + (0 + .is_male + .is_female |      coupleID)
Data: df_apim_count

      AIC       BIC    logLik -2*log(L)  df.resid 
     5434      5479     -2708      5416      1091 

Random effects:

Conditional model:
 Groups   Name       Variance Std.Dev. Corr 
 coupleID .is_male   0.14072  0.3751        
          .is_female 0.07157  0.2675   0.30 
Number of obs: 1100, groups:  coupleID, 550

Conditional model:
                                  Estimate Std. Error z value Pr(>|z|)    
.is_male                           0.47400    0.10744   4.412 1.02e-05 ***
.is_female                         0.63559    0.08969   7.087 1.37e-12 ***
.is_male:.communication_actor      0.20833    0.01815  11.476  < 2e-16 ***
.is_female:.communication_actor    0.17450    0.01541  11.322  < 2e-16 ***
.is_male:.communication_partner    0.00901    0.01854   0.486 0.626980    
.is_female:.communication_partner  0.05710    0.01481   3.856 0.000115 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Distinguishable Dyads - Check Model Convergence and Fit

Check Rhats and Effective Sample Sizes (ESS_tail and ESS_bulk) directly from the brms summary. Additionally you can (using Model B as an example):

# In brms 2.23, draws from either Stan backend are stored as a stanfit here.
rstan::check_hmc_diagnostics(model_dist_apim_b$fit)

Divergences:
0 of 4000 iterations ended with a divergence.

Tree depth:
0 of 4000 iterations saturated the maximum tree depth of 10.

Energy:
E-BFMI indicated no pathological behavior.

Distinguishable Dyads - Check Model Convergence and Fit

plot(model_dist_apim_b, ask = FALSE)

Distinguishable Dyads - Check Model Convergence and Fit

pp_check(model_dist_apim_b, 'dens_overlay_grouped', group = '.is_male')

Distinguishable Dyads - Check Model Convergence and Fit

pp_check(model_dist_apim_b, 'ecdf_overlay_grouped', group = '.is_male')

Distinguishable Dyads - Check Model Convergence and Fit

# Custom function to make DHARMa work with brms (see file 'Functions')
DHARMa.check_brms(model_dist_apim_b)

Distinguishable Dyads - Check Model Convergence and Fit

  • Diagnose convergence, residual fit, and predictive accuracy separately; check coding, scaling, priors, and identifiability before changing the model.
  • Choose the family and link from the outcome and substantive assumptions—for example, cumulative() for ordered outcomes or bernoulli() for binary data.
  • If chains are healthy, more post-warmup draws can raise ESS. More warmup helps adaptation, while higher adapt_delta specifically targets divergences (Bürkner, 2025).

Excursion: Generalized Dyadic Models

For Gaussian outcomes, residual covariance has a direct response-scale interpretation:

  • partner-specific residual variances
  • partner residual correlations
  • same-day residual covariance in ILD data

For non-Gaussian outcomes, the same syntax usually does not estimate an ordinary residual covariance.

Instead, random-effect covariance terms represent latent Gaussian heterogeneity on the link scale.

Excursion: Count vs Binary Outcomes

For count outcomes, latent partner covariance is often useful and can be estimable:

  • Poisson or negative-binomial models can include dyad-level latent covariance terms
  • estimates live on the log-link scale
  • covariance parameters need diagnostics and sensitivity checks

For dichotomous outcomes, the same idea is much harder:

  • each partner contributes only 0 or 1 per occasion
  • many latent variance-correlation combinations can imply similar pair probabilities
  • same-day binary covariance is especially weakly identified

Excursion: Practical Warning

Binary dyadic covariance models can converge mechanically while still being unreliable.

Watch for:

  • variance estimates near zero or extremely large
  • correlations near -1 or 1
  • low ESS, divergent transitions, or strong prior sensitivity
  • different software giving very different covariance estimates

Excursion: Practical Strategy

Practical default:

  • estimate simpler structures first
  • use informative priors if fitting latent binary covariance in brms
  • treat same-day binary partner covariance as a specialized model, not a routine extension
  • if the main goal is fixed effects, consider simpler random effects or robust/clustered inference

Distinguishable Dyads - Results

Est. 95% CI Rhat Bulk_ESS Tail_ESS
Fixed Effects
.is_male -5.44 [-5.99, -4.88] 1.003 2458 2617
.is_female -3.90 [-4.37, -3.47] 1.001 2401 2641
.is_male:.communication_actor 1.84 [ 1.75, 1.94] 1.000 2678 2997
.is_female:.communication_actor 1.59 [ 1.51, 1.67] 1.001 2371 2549
.is_male:.communication_partner 0.17 [ 0.07, 0.26] 1.001 2555 2726
.is_female:.communication_partner 0.31 [ 0.23, 0.39] 1.000 2398 2724
Residual Structure
cortime(female,male) 0.51 [0.44, 0.57] 1.000 2895 2470
sigma_.is_male 1.60 [1.51, 1.71] 1.000 3555 2827
sigma_.is_female 1.30 [1.23, 1.39] 1.001 3317 2788

NOTE: brms estimates sigma on the log scale. In the tables, sigmas are exponentiated to show residual SDs on the response scale, like glmmTMB output.

Distinguishable Dyads

Distinguishable Dyads - Intensive Longitudinal APIM (L-APIM)

Distinguishable Dyads - L-APIM - Raw Data

print_df(head(df_long))
userID coupleID diaryday gender closeness provided_support
1_1 1 0 female 5.659395 4.733482
1_1 1 1 female 7.249842 5.019950
1_1 1 2 female 6.040587 5.365179
1_1 1 3 female 6.898253 5.148704
1_1 1 4 female 5.350049 3.702039
1_1 1 5 female 6.191282 4.844602

Distinguishable Dyads - L-APIM - Data Preparation

df_long_apim <- prepare_dyad_data(
  data = df_long,
  dyad = coupleID,
  member = userID,
  role = gender,
  time = diaryday,
  predictors = provided_support,
  model_types = "apim"
) %>%
  mutate(
    diaryday_c = diaryday - mean(diaryday)
  )

Distinguishable Dyads - L-APIM - Prepared Data

userID coupleID diaryday gender closeness provided_support .composition .composition_role .is_female .is_male .provided_support_cwp .provided_support_cbp .provided_support_actor .provided_support_partner .provided_support_cwp_actor .provided_support_cwp_partner .provided_support_cbp_actor .provided_support_cbp_partner diaryday_c
31_1 31 0 female 5.47 6.02 female_x_male female_x_male_female 1 0 0.41 0.65 6.02 7.09 0.41 1.54 0.65 0.58 -27
31_1 31 1 female 7.69 6.13 female_x_male female_x_male_female 1 0 0.51 0.65 6.13 6.95 0.51 1.4 0.65 0.58 -26
31_1 31 2 female 5.83 6.1 female_x_male female_x_male_female 1 0 0.49 0.65 6.1 6.34 0.49 0.79 0.65 0.58 -25
31_1 31 3 female 6.55 6.99 female_x_male female_x_male_female 1 0 1.38 0.65 6.99 5.39 1.38 -0.15 0.65 0.58 -24
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
31_2 31 0 male 5.19 7.09 female_x_male female_x_male_male 0 1 1.54 0.58 7.09 6.02 1.54 0.41 0.58 0.65 -27
31_2 31 1 male 6.96 6.95 female_x_male female_x_male_male 0 1 1.4 0.58 6.95 6.13 1.4 0.51 0.58 0.65 -26
31_2 31 2 male 6.9 6.34 female_x_male female_x_male_male 0 1 0.79 0.58 6.34 6.1 0.79 0.49 0.58 0.65 -25
31_2 31 3 male 7.59 5.39 female_x_male female_x_male_male 0 1 -0.15 0.58 5.39 6.99 -0.15 1.38 0.58 0.65 -24

Distinguishable Dyads - L-APIM - Centering

With repeated measures we need to center variables in order to not conflate levels of analysis (Bolger & Laurenceau, 2013):

  • Between-Person: The person-mean in relation to the grand mean
  • Within-Person: Daily fluctuations of an individual from their person-mean

Although the data contain couple, person, and day structure, the APIM is written at the person-day level after decomposing predictors into between-person and within-person components. Dyad-level dependence is still modelled explicitly through partner-specific random effects and residual covariance structures.

Distinguishable Dyads - L-APIM - Centering

prepare_dyad_data() performs this predictor decomposition before creating the actor and partner columns. The centered diaryday_c variable above is a separate time transformation for the models, not part of the dyadic predictor construction.

Distinguishable Dyads - L-APIM - Model nlme

nlme model of this data for reference (simplified random effects structure to achieve convergence).

This package is old, but still often used and cited. So here, we briefly show it to demonstrate that we can achieve the same model with much more modern and powerful methods.

library(nlme)

model1 <- lme(
  fixed =
    # Intercepts
    closeness ~ 0 + .is_male + .is_female +

    # Time Slopes
    .is_male:diaryday_c +
    .is_female:diaryday_c +

    # Between-Person APIM
    .is_male:.provided_support_cbp_actor +
    .is_male:.provided_support_cbp_partner +

    .is_female:.provided_support_cbp_actor +
    .is_female:.provided_support_cbp_partner +

    # Within-Person APIM
    .is_male:.provided_support_cwp_actor +
    .is_male:.provided_support_cwp_partner +

    .is_female:.provided_support_cwp_actor +
    .is_female:.provided_support_cwp_partner,

  random =
    ~ 0 + .is_male + .is_female

    + .is_male:diaryday_c + .is_female:diaryday_c

  # in nlme we need to remove many random slopes to achieve convergence
  # alternatively, we could try to remove correlations between the slopes.
    #+ .is_male:.provided_support_cwp_actor
    #+ .is_male:.provided_support_cwp_partner
    #+ .is_female:.provided_support_cwp_actor
    #+ .is_female:.provided_support_cwp_partner
  | coupleID,

  weights = varIdent(form = ~ 1 | gender), # heterogeneous residual variances
  corr = corCompSymm(form = ~ 1 | coupleID/diaryday), # compound symmetry (partner's values on each day)

  data = df_long_apim,
  control = list(maxIter=1000)
)

summary(model1)

Distinguishable Dyads - L-APIM - Model glmmTMB

model_dist_apim_long_glmmtmb <- glmmTMB(
  closeness ~ 0 + .is_male + .is_female +

    # Time Slopes
    .is_male:diaryday_c +
    .is_female:diaryday_c +

    # Between-Person APIM
    .is_male:.provided_support_cbp_actor +
    .is_male:.provided_support_cbp_partner +

    .is_female:.provided_support_cbp_actor +
    .is_female:.provided_support_cbp_partner +

    # Within-Person APIM
    .is_male:.provided_support_cwp_actor +
    .is_male:.provided_support_cwp_partner +

    .is_female:.provided_support_cwp_actor +
    .is_female:.provided_support_cwp_partner +

  # Accounting for non-independence between partners' means and trajectories
  # and effect sensitivities via random effects:
  (0 + .is_male + .is_female
    + .is_male:diaryday_c + .is_female:diaryday_c
  | coupleID) +

  # Accounting for daily non-independence
  us(0 + .is_male + .is_female | coupleID:diaryday)

  # heteroscedastic residuals
  , dispformula = ~ 0
  , data = df_long_apim
)

Distinguishable Dyads - L-APIM - Model glmmTMB

 Family: gaussian  ( identity )
Formula:          closeness ~ 0 + .is_male + .is_female + .is_male:diaryday_c +  
    .is_female:diaryday_c + .is_male:.provided_support_cbp_actor +  
    .is_male:.provided_support_cbp_partner + .is_female:.provided_support_cbp_actor +  
    .is_female:.provided_support_cbp_partner + .is_male:.provided_support_cwp_actor +  
    .is_male:.provided_support_cwp_partner + .is_female:.provided_support_cwp_actor +  
    .is_female:.provided_support_cwp_partner + (0 + .is_male +  
    .is_female + .is_male:diaryday_c + .is_female:diaryday_c |  
    coupleID) + us(0 + .is_male + .is_female | coupleID:diaryday)
Dispersion:                 ~0
Data: df_long_apim

      AIC       BIC    logLik -2*log(L)  df.resid 
  30001.8   30184.4  -14975.9   29951.8     10975 

Random effects:

Conditional model:
 Groups            Name                  Variance  Std.Dev. Corr              
 coupleID          .is_male              0.6813649 0.82545                    
                   .is_female            0.7929829 0.89050   0.19             
                   .is_male:diaryday_c   0.0002637 0.01624   0.55  0.03       
                   .is_female:diaryday_c 0.0002571 0.01603   0.21  0.54 -0.01 
 coupleID:diaryday .is_male              1.2583395 1.12176                    
                   .is_female            0.5206603 0.72157  0.04              
Number of obs: 11000, groups:  coupleID, 100; coupleID:diaryday, 5500

Conditional model:
                                          Estimate Std. Error z value Pr(>|z|)    
.is_male                                  4.486375   0.085960   52.19  < 2e-16 ***
.is_female                                5.736720   0.091763   62.52  < 2e-16 ***
.is_male:diaryday_c                      -0.003258   0.001883   -1.73  0.08355 .  
.is_female:diaryday_c                     0.010606   0.001717    6.18 6.46e-10 ***
.is_male:.provided_support_cbp_actor      1.128018   0.101922   11.07  < 2e-16 ***
.is_male:.provided_support_cbp_partner    0.298902   0.098817    3.02  0.00249 ** 
.is_female:.provided_support_cbp_actor    1.380753   0.105607   13.07  < 2e-16 ***
.is_female:.provided_support_cbp_partner  0.587254   0.108927    5.39 7.00e-08 ***
.is_male:.provided_support_cwp_actor      0.184479   0.020564    8.97  < 2e-16 ***
.is_male:.provided_support_cwp_partner    0.124548   0.020571    6.05 1.41e-09 ***
.is_female:.provided_support_cwp_actor    0.419162   0.013277   31.57  < 2e-16 ***
.is_female:.provided_support_cwp_partner  0.170127   0.013270   12.82  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Comparing outputs between nlme and glmmTMB

Let’s compare the random effects outputs to see if the models are truly equivalent.

Residual variance in nlme must first be constructed from the output by taking the residual SD and multiplying it by the parameter estimate under the variance function for each partner.

nlme: 
                     StdDev     Corr
.is_male              0.83645552 is_mal is_fml is_m:_
.is_female            0.90199634  0.189
.is_male:diaryday_c   0.01634796  0.544  0.033
.is_female:diaryday_c 0.01612665  0.204  0.534 -0.012
Residual             0.72170702

Correlation Structure: Compound symmetry
 Formula: ~1 | coupleID/diaryday
 Parameter estimate(s):
       Rho
0.04094751
Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | gender
 Parameter estimates:
       1        2
1.000000 1.554608



Manually computed by multiplying the residual SD by the multipliers:

Residual SDs by partner:
  Partner 1: 0.721707
  Partner 2: 1.121972
glmmTMB: 
Random effects:
Conditional model:
 Groups            Name                 Variance  Std.Dev. Corr
 coupleID          .is_male              0.6814498 0.82550
                   .is_female            0.7930407 0.89053   0.19
                   .is_male:diaryday_c   0.0002637 0.01624   0.55  0.03
                   .is_female:diaryday_c 0.0002571 0.01603   0.21  0.54 -0.01
 coupleID:diaryday .is_male              1.2583552 1.12176
                   .is_female            0.5206608 0.72157  0.04

Distinguishable Dyads - L-APIM - Model brms

Equivalent brms model.

In contrast to glmmTMB, brms provides native variance-covariance structures to model residuals directly. Therefore, we do NOT set the residuals to zero!

formula_dist_apim_long_simple <- bf(
  closeness ~ 0 + .is_male + .is_female +

    # Time Slopes
    .is_male:diaryday_c +
    .is_female:diaryday_c +

    # Between-Person APIM
    .is_male:.provided_support_cbp_actor +
    .is_male:.provided_support_cbp_partner +

    .is_female:.provided_support_cbp_actor +
    .is_female:.provided_support_cbp_partner +

    # Within-Person APIM
    .is_male:.provided_support_cwp_actor +
    .is_male:.provided_support_cwp_partner +

    .is_female:.provided_support_cwp_actor +
    .is_female:.provided_support_cwp_partner +

  # Accounting for non-independence between partners' means and trajectories
  # and effect sensitivities via random effects:
  (0 + .is_male + .is_female
    + .is_male:diaryday_c + .is_female:diaryday_c
  | coupleID) +
    
  # Gender indexes the two partner residual positions within each couple-day.
  unstr(time = gender, gr = coupleID:diaryday)

  # heteroscedastic residuals
  , sigma ~ 0 + .is_male + .is_female
)

# brms' defaults for group-level SDs and correlations are sensible weak defaults
# (half Student-t priors for SDs and LKJ priors for correlation matrices). Here we
# set weakly regularizing priors explicitly so the modelling assumptions are visible.

priors_dist_apim_long_simple <- c(
  prior(normal(4, 2), class = "b", coef = ".is_male"), # male intercept
  prior(normal(4, 2), class = "b", coef = ".is_female"), # female intercept
  prior(normal(0, 2), class = "b"), # other beta coefficients
  prior(exponential(1), class = "sd"),
  prior(lkj(2), class = "cor"),
  prior(normal(0, 0.7), class = "b", dpar = "sigma"),
  prior(lkj(2), class = "cortime")
)

model_dist_apim_long_simple <- brm(
  formula = formula_dist_apim_long_simple,
  data = df_long_apim,
  family = gaussian(link = identity), # student() is also often a nice robust option.
  prior = priors_dist_apim_long_simple,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000,
  seed = 123,
  file = file.path('brms_cache', 'example1_dist_apim_long_simpel') # Cache the model
)

Distinguishable Dyads - L-APIM - Results (scrollable)

Est. 95% CI Rhat Bulk_ESS Tail_ESS
Fixed Effects
.is_male 4.49 [ 4.32, 4.65] 1.007 324 794
.is_female 5.74 [ 5.56, 5.91] 1.009 355 656
.is_male:diaryday_c 0.00 [-0.01, 0.00] 1.007 1072 2113
.is_female:diaryday_c 0.01 [ 0.01, 0.01] 1.002 880 1778
.is_male:.provided_support_cbp_actor 1.14 [ 0.91, 1.34] 1.006 489 791
.is_male:.provided_support_cbp_partner 0.30 [ 0.10, 0.51] 1.004 553 716
.is_female:.provided_support_cbp_actor 1.38 [ 1.18, 1.58] 1.005 507 1077
.is_female:.provided_support_cbp_partner 0.60 [ 0.36, 0.81] 1.005 440 851
.is_male:.provided_support_cwp_actor 0.18 [ 0.14, 0.23] 1.000 5328 3057
.is_male:.provided_support_cwp_partner 0.12 [ 0.08, 0.16] 1.003 7228 3389
.is_female:.provided_support_cwp_actor 0.42 [ 0.39, 0.45] 1.002 5990 3068
.is_female:.provided_support_cwp_partner 0.17 [ 0.14, 0.20] 1.001 5829 3105
Random Effects
sd(.is_male) 0.84 [ 0.73, 0.97] 1.010 520 1077
sd(.is_female) 0.92 [ 0.80, 1.06] 1.006 556 998
sd(.is_male:diaryday_c) 0.02 [ 0.01, 0.02] 1.001 1744 2637
sd(.is_female:diaryday_c) 0.02 [ 0.01, 0.02] 1.001 1358 2043
cor(.is_male,.is_female) 0.18 [-0.05, 0.35] 1.010 321 688
cor(.is_male,.is_male:diaryday_c) 0.51 [ 0.32, 0.67] 1.004 1528 1804
cor(.is_female,.is_male:diaryday_c) 0.03 [-0.21, 0.24] 1.002 1315 2204
cor(.is_male,.is_female:diaryday_c) 0.19 [-0.01, 0.38] 1.003 753 1691
cor(.is_female,.is_female:diaryday_c) 0.51 [ 0.33, 0.65] 1.002 1073 1832
cor(.is_male:diaryday_c,.is_female:diaryday_c) -0.02 [-0.25, 0.22] 1.008 867 1748
cortime(female,male) 0.04 [ 0.01, 0.07] 1.000 6904 3067
Residual Structure
sigma_.is_male 1.12 [ 1.10, 1.14] 1.001 6996 2514
sigma_.is_female 0.72 [ 0.71, 0.74] 0.999 6880 2868

Note: Sigma rows are shown on the response scale; they were exponentiated from brms’ log-sigma parameterization.

Distinguishable Dyads - L-APIM - AR(1)?

In this model, unstr(time = gender, gr = coupleID:diaryday) handles same-day partner interdependence. Current univariate brms models cannot combine this residual correlation structure with ar() in the same formula (Bürkner, 2025).

Standard glmmTMB can add a same-day dyadic block and independent person-specific AR(1) processes:

us(0 + gender | coupleID:diaryday) +
  ar1(0 + diaryday_factor | userID)

This additive model represents same-day dyadic shocks plus person-specific persistence. It splits lag-zero variance across the two blocks. These residual blocks contribute zero cross-partner covariance across different days, although dyad-level random effects can still induce cross-day covariance; this is not one separable partner-by-time residual process (Brooks et al., 2026).

Excursion: Dynamic Workaround

A dynamic MLM can account for temporal dependence by including lagged within-person centered outcomes as predictors (include both actor inertia and partner spillover).

But:

  • The APIM coefficients become effects conditional on prior closeness, not total within-person associations.
  • Lagged-outcome models can also suffer from dynamic-panel bias, often called Nickell bias (Nickell, 1981), because the lagged outcome is related to person/couple effects and prior residuals. This is especially problematic with short time series.
    • One solution for the Nickell bias is to use multilevel DSEM through MPLUS.

Maximal Model Specification (non-dynamic/no AR(1))

We add all random slopes and their correlations. This is neither converging in nlme nor glmmTMB.

formula_dist_apim_long_complex <- bf(
  closeness ~

    # Intercepts
    0 + .is_male + .is_female +

    # Time Slopes
    .is_male:diaryday_c +
    .is_female:diaryday_c +

    # Between-Person APIM
    .is_male:.provided_support_cbp_actor +
    .is_male:.provided_support_cbp_partner +

    .is_female:.provided_support_cbp_actor +
    .is_female:.provided_support_cbp_partner +

    # Within-Person APIM
    .is_male:.provided_support_cwp_actor +
    .is_male:.provided_support_cwp_partner +

    .is_female:.provided_support_cwp_actor +
    .is_female:.provided_support_cwp_partner +

    # Accounting for non-independence between partners' means and trajectories
    # and effect sensitivities via random effects:
    (0 + .is_male + .is_female +

       .is_male:diaryday_c + .is_female:diaryday_c +

       .is_male:.provided_support_cwp_actor +
       .is_male:.provided_support_cwp_partner +
       .is_female:.provided_support_cwp_actor +
       .is_female:.provided_support_cwp_partner
    | coupleID ) +

    # Accounting for daily non-independence
    unstr(time = gender, gr = coupleID:diaryday),

  # No more AR1. You could also test MA or ARMA.

  sigma ~ 0 + .is_male + .is_female   # heteroscedastic residuals
)

priors_dist_apim_long_complex <- c(
  prior(normal(4, 2), class = "b", coef = ".is_male"),
  prior(normal(4, 2), class = "b", coef = ".is_female"),
  prior(normal(0, 2), class = "b"),
  prior(exponential(1), class = "sd"),
  prior(lkj(2), class = "cor"),
  prior(normal(0, 0.7), class = "b", dpar = "sigma"),
  prior(lkj(2), class = "cortime")
)


model_dist_apim_long_complex <- brm(
  formula = formula_dist_apim_long_complex,
  data = df_long_apim,
  family = gaussian(link = identity),
  prior = priors_dist_apim_long_complex,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000,
  seed = 123,
  file = file.path('brms_cache', 'example1_dist_apim_long_complex') # Cache the model
)

Maximal Model Specification: Results

Est. 95% CI Rhat Bulk_ESS Tail_ESS
Fixed Effects
.is_male 4.49 [ 4.31, 4.67] 1.006 334 722
.is_female 5.73 [ 5.54, 5.92] 1.006 271 770
.is_male:diaryday_c 0.00 [-0.01, 0.00] 1.002 1249 1965
.is_female:diaryday_c 0.01 [ 0.01, 0.01] 1.005 840 2033
.is_male:.provided_support_cbp_actor 1.12 [ 0.92, 1.32] 1.004 688 1529
.is_male:.provided_support_cbp_partner 0.31 [ 0.10, 0.51] 1.002 656 1267
.is_female:.provided_support_cbp_actor 1.37 [ 1.16, 1.59] 1.002 491 987
.is_female:.provided_support_cbp_partner 0.60 [ 0.35, 0.82] 1.003 432 740
.is_male:.provided_support_cwp_actor 0.18 [ 0.14, 0.23] 1.000 5755 3046
.is_male:.provided_support_cwp_partner 0.13 [ 0.08, 0.18] 1.000 6028 3158
.is_female:.provided_support_cwp_actor 0.42 [ 0.39, 0.45] 1.000 5549 3243
.is_female:.provided_support_cwp_partner 0.17 [ 0.14, 0.20] 1.000 5548 3014
Random Effects
sd(.is_male) 0.85 [ 0.74, 0.98] 1.002 857 1449
sd(.is_female) 0.91 [ 0.79, 1.05] 1.006 515 1265
sd(.is_male:diaryday_c) 0.02 [ 0.01, 0.02] 1.000 1973 2634
sd(.is_female:diaryday_c) 0.02 [ 0.01, 0.02] 1.001 1361 2175
sd(.is_male:.provided_support_cwp_actor) 0.09 [ 0.02, 0.15] 1.002 1162 1085
sd(.is_male:.provided_support_cwp_partner) 0.14 [ 0.06, 0.20] 1.000 1012 947
sd(.is_female:.provided_support_cwp_actor) 0.08 [ 0.02, 0.12] 1.004 923 905
sd(.is_female:.provided_support_cwp_partner) 0.06 [ 0.01, 0.10] 1.001 1002 1720
cor(.is_male,.is_female) 0.15 [-0.04, 0.34] 1.022 350 819
cor(.is_male,.is_male:diaryday_c) 0.48 [ 0.28, 0.63] 1.001 2211 2867
cor(.is_female,.is_male:diaryday_c) 0.00 [-0.20, 0.22] 1.003 1801 3156
cor(.is_male,.is_female:diaryday_c) 0.17 [-0.03, 0.36] 1.005 940 2321
cor(.is_female,.is_female:diaryday_c) 0.49 [ 0.30, 0.64] 1.008 1172 2078
cor(.is_male:diaryday_c,.is_female:diaryday_c) -0.03 [-0.26, 0.19] 1.003 1108 1730
cor(.is_male,.is_male:.provided_support_cwp_actor) 0.37 [-0.06, 0.70] 1.001 4857 2732
cor(.is_female,.is_male:.provided_support_cwp_actor) 0.19 [-0.22, 0.57] 1.003 5018 2783
cor(.is_male:diaryday_c,.is_male:.provided_support_cwp_actor) 0.33 [-0.13, 0.68] 1.001 5222 3270
cor(.is_female:diaryday_c,.is_male:.provided_support_cwp_actor) 0.27 [-0.17, 0.64] 1.000 4899 2758
cor(.is_male,.is_male:.provided_support_cwp_partner) 0.06 [-0.26, 0.38] 1.000 5278 3217
cor(.is_female,.is_male:.provided_support_cwp_partner) 0.10 [-0.22, 0.43] 1.001 4991 2571
cor(.is_male:diaryday_c,.is_male:.provided_support_cwp_partner) 0.04 [-0.31, 0.40] 1.000 4298 2455
cor(.is_female:diaryday_c,.is_male:.provided_support_cwp_partner) 0.11 [-0.22, 0.45] 1.000 4766 2466
cor(.is_male:.provided_support_cwp_actor,.is_male:.provided_support_cwp_partner) 0.24 [-0.27, 0.68] 1.005 1124 2137
cor(.is_male,.is_female:.provided_support_cwp_actor) -0.06 [-0.44, 0.29] 1.001 4620 2478
cor(.is_female,.is_female:.provided_support_cwp_actor) 0.21 [-0.17, 0.55] 1.000 4730 2276
cor(.is_male:diaryday_c,.is_female:.provided_support_cwp_actor) -0.14 [-0.51, 0.25] 1.001 4056 3026
cor(.is_female:diaryday_c,.is_female:.provided_support_cwp_actor) -0.07 [-0.43, 0.30] 1.000 5099 3132
cor(.is_male:.provided_support_cwp_actor,.is_female:.provided_support_cwp_actor) -0.17 [-0.61, 0.37] 1.002 1540 2245
cor(.is_male:.provided_support_cwp_partner,.is_female:.provided_support_cwp_actor) -0.10 [-0.56, 0.40] 1.001 2566 2969
cor(.is_male,.is_female:.provided_support_cwp_partner) 0.27 [-0.21, 0.63] 1.000 5732 2872
cor(.is_female,.is_female:.provided_support_cwp_partner) 0.07 [-0.37, 0.48] 1.000 5615 2735
cor(.is_male:diaryday_c,.is_female:.provided_support_cwp_partner) 0.12 [-0.34, 0.53] 1.000 5813 3459
cor(.is_female:diaryday_c,.is_female:.provided_support_cwp_partner) 0.15 [-0.31, 0.56] 1.002 5869 2966
cor(.is_male:.provided_support_cwp_actor,.is_female:.provided_support_cwp_partner) 0.16 [-0.39, 0.64] 1.001 2540 3043
cor(.is_male:.provided_support_cwp_partner,.is_female:.provided_support_cwp_partner) 0.06 [-0.48, 0.57] 1.002 3413 3154
cor(.is_female:.provided_support_cwp_actor,.is_female:.provided_support_cwp_partner) 0.11 [-0.44, 0.60] 1.001 3345 3240
cortime(female,male) 0.04 [ 0.01, 0.07] 1.002 8477 2860
Residual Structure
sigma_.is_male 1.11 [ 1.09, 1.14] 1.002 6681 2536
sigma_.is_female 0.72 [ 0.70, 0.73] 1.001 6554 2959

Maximal Model: Diagnostics and Prediction Target

  • Check convergence and model fit before comparing predictive accuracy.
  • loo() holds out one observation while retaining its partner and the dyad’s other days. This evaluates conditional observation-level prediction, not a completely new dyad.
  • For Gaussian unstr() models, brms computes the required conditional density \(p(y_i\mid y_{-i},\theta)\) for non-factorizable LOO (Bürkner et al., 2021).
  • New-dyad prediction would instead require a grouped validation design (Vehtari et al., 2017).

Maximal Model: PSIS-LOO Diagnostics

loo1 <- loo(model_dist_apim_long_simple)
loo2 <- loo(model_dist_apim_long_complex)

loo::pareto_k_table(loo1)
loo::pareto_k_table(loo2)
All Pareto k estimates are good (k < 0.7).

All Pareto k estimates are good (k < 0.7).

Pareto \(k\) diagnoses the importance-sampling approximation; high values are not themselves evidence of overfitting (Vehtari et al., 2017).

Maximal Model: LOO Comparison

a <- loo_compare(loo1, loo2)
print(a)
                             elpd_diff se_diff
model_dist_apim_long_complex   0.0       0.0
model_dist_apim_long_simple  -15.4       6.7

The complex model predicts the held-out observations better by 15.4 ELPD units (SE 6.7). This supports its additional structure for this prediction target, but does not prove that every variance or correlation is necessary; ELPD differences should be interpreted with their uncertainty (Sivula et al., 2025). If performance were similar, theory, identifiability, and parsimony should guide simplification (Barr et al., 2013; del Rosario & West, 2025).

Maximal Model Specification: Comparing Estimates to Simulated Values (Ground Truth)

Group Parameter True (Simulated) Estimate 95% CI Coverage Abs. Error Recovery
Fixed Intercept (male) 4.5659893 4.4890395 [4.305, 4.668] Yes 0.077 Covered
Fixed Intercept (female) 5.7823087 5.7316403 [5.537, 5.919] Yes 0.051 Covered
Fixed Time slope (male) -0.0050000 -0.0029689 [-0.007, 0.001] Yes 0.002 Covered
Fixed Time slope (female) 0.0100000 0.0104637 [0.007, 0.014] Yes 0.000 Covered
Between-person Actor BP (male) 1.1855503 1.1204791 [0.917, 1.322] Yes 0.065 Covered
Between-person Partner BP (male) 0.2707719 0.3078625 [0.104, 0.511] Yes 0.037 Covered
Between-person Actor BP (female) 1.4298780 1.3697538 [1.156, 1.594] Yes 0.060 Covered
Between-person Partner BP (female) 0.5291829 0.5924361 [0.351, 0.817] Yes 0.063 Covered
Within-person Actor WP (male) 0.2000000 0.1812234 [0.135, 0.227] Yes 0.019 Covered
Within-person Partner WP (male) 0.1000000 0.1270781 [0.079, 0.176] Yes 0.027 Covered
Within-person Actor WP (female) 0.4000000 0.4205947 [0.390, 0.451] Yes 0.021 Covered
Within-person Partner WP (female) 0.2000000 0.1713545 [0.143, 0.200] Yes 0.029 Covered
Random effects (SD) Couple intercept SD (male) 0.7919045 0.8495693 [0.740, 0.984] Yes 0.058 Covered
Random effects (SD) Couple intercept SD (female) 0.8979616 0.9154796 [0.793, 1.055] Yes 0.018 Covered
Random effects (SD) Time slope SD (male) 0.0149848 0.0167471 [0.014, 0.020] Yes 0.002 Covered
Random effects (SD) Time slope SD (female) 0.0150338 0.0164231 [0.014, 0.019] Yes 0.001 Covered
Random effects (SD) Actor WP slope SD (male) 0.1006336 0.0899573 [0.021, 0.153] Yes 0.011 Covered
Random effects (SD) Partner WP slope SD (male) 0.0806085 0.1344609 [0.055, 0.201] Yes 0.054 Covered
Random effects (SD) Actor WP slope SD (female) 0.0955935 0.0759606 [0.019, 0.120] Yes 0.020 Covered
Random effects (SD) Partner WP slope SD (female) 0.0678586 0.0542037 [0.005, 0.102] Yes 0.014 Covered
Random effects (cor) cor(Intercept (male), Intercept (female)) 0.1803231 0.1521969 [-0.042, 0.335] Yes 0.028 Covered
Random effects (cor) cor(Intercept (male), Time slope (male)) 0.5083773 0.4708544 [0.283, 0.632] Yes 0.038 Covered
Random effects (cor) cor(Intercept (male), Time slope (female)) 0.1921701 0.1708417 [-0.033, 0.365] Yes 0.021 Covered
Random effects (cor) cor(Intercept (male), Actor WP slope (male)) 0.2262539 0.3543439 [-0.059, 0.696] Yes 0.128 Covered
Random effects (cor) cor(Intercept (male), Partner WP slope (male)) 0.0953425 0.0605446 [-0.260, 0.378] Yes 0.035 Covered
Random effects (cor) cor(Intercept (male), Actor WP slope (female)) 0.0398182 -0.0614619 [-0.436, 0.294] Yes 0.101 Covered
Random effects (cor) cor(Intercept (male), Partner WP slope (female)) 0.3479821 0.2520320 [-0.206, 0.629] Yes 0.096 Covered
Random effects (cor) cor(Intercept (female), Time slope (male)) 0.0526005 0.0061449 [-0.201, 0.221] Yes 0.046 Covered
Random effects (cor) cor(Intercept (female), Time slope (female)) 0.6138595 0.4831640 [0.303, 0.636] Yes 0.131 Covered
Random effects (cor) cor(Intercept (female), Actor WP slope (male)) 0.1319990 0.1866849 [-0.225, 0.569] Yes 0.055 Covered
Random effects (cor) cor(Intercept (female), Partner WP slope (male)) -0.0241164 0.1010175 [-0.221, 0.429] Yes 0.125 Covered
Random effects (cor) cor(Intercept (female), Actor WP slope (female)) 0.1385319 0.2046311 [-0.172, 0.549] Yes 0.066 Covered
Random effects (cor) cor(Intercept (female), Partner WP slope (female)) -0.0263927 0.0667585 [-0.371, 0.485] Yes 0.093 Covered
Random effects (cor) cor(Time slope (male), Time slope (female)) 0.1339757 -0.0304431 [-0.257, 0.190] Yes 0.164 Covered
Random effects (cor) cor(Time slope (male), Actor WP slope (male)) 0.2400640 0.3130248 [-0.133, 0.682] Yes 0.073 Covered
Random effects (cor) cor(Time slope (male), Partner WP slope (male)) 0.0918211 0.0427317 [-0.313, 0.397] Yes 0.049 Covered
Random effects (cor) cor(Time slope (male), Actor WP slope (female)) 0.0024579 -0.1337158 [-0.512, 0.253] Yes 0.136 Covered
Random effects (cor) cor(Time slope (male), Partner WP slope (female)) 0.2635152 0.1105994 [-0.343, 0.527] Yes 0.153 Covered
Random effects (cor) cor(Time slope (female), Actor WP slope (male)) 0.0711586 0.2634154 [-0.166, 0.641] Yes 0.192 Covered
Random effects (cor) cor(Time slope (female), Partner WP slope (male)) 0.0441926 0.1093859 [-0.225, 0.448] Yes 0.065 Covered
Random effects (cor) cor(Time slope (female), Actor WP slope (female)) 0.1380171 -0.0687735 [-0.428, 0.300] Yes 0.207 Covered
Random effects (cor) cor(Time slope (female), Partner WP slope (female)) 0.0279326 0.1470529 [-0.313, 0.558] Yes 0.119 Covered
Random effects (cor) cor(Actor WP slope (male), Partner WP slope (male)) 0.1128802 0.2345585 [-0.266, 0.682] Yes 0.122 Covered
Random effects (cor) cor(Actor WP slope (male), Actor WP slope (female)) 0.0478755 -0.1584750 [-0.610, 0.370] Yes 0.206 Covered
Random effects (cor) cor(Actor WP slope (male), Partner WP slope (female)) 0.1135981 0.1524683 [-0.386, 0.639] Yes 0.039 Covered
Random effects (cor) cor(Partner WP slope (male), Actor WP slope (female)) -0.0107839 -0.0987572 [-0.562, 0.404] Yes 0.088 Covered
Random effects (cor) cor(Partner WP slope (male), Partner WP slope (female)) 0.1384787 0.0602770 [-0.484, 0.573] Yes 0.078 Covered
Random effects (cor) cor(Actor WP slope (female), Partner WP slope (female)) 0.1494376 0.1005532 [-0.442, 0.599] Yes 0.049 Covered
Residuals Same-day residual correlation 0.0366176 0.0413748 [0.015, 0.068] Yes 0.005 Covered
Residuals Residual SD (male) (exp) 1.1210517 1.1147509 [1.093, 1.136] Yes 0.006 Covered
Residuals Residual SD (female) (exp) 0.7173353 0.7181891 [0.705, 0.732] Yes 0.001 Covered

Exchangeable Dyads

Exchangeable Dyads - Cross-Sectional

Exchangeable Dyads - Cross-Sectional APIM

Exchangeable APIM (e.g., del Rosario & West, 2025; Kenny & Ackerman, 2023)

Exchangeable Dyads - Cross-Sectional APIM

Exchangeable APIM (e.g., del Rosario & West, 2025; Kenny & Ackerman, 2023)

Exchangeable Dyads - Main Assumptions

Partners are exchangeable, i.e., not systematically different.

  • Equal actor effects
  • Equal partner effects
  • Equal means
  • Equal residual variances

But they should still be allowed to vary within each couple, while being correlated.

(e.g., del Rosario & West, 2025; Kenny & Ackerman, 2023)

Exchangeable Dyads - Cross-Sectional APIM: Assumption

Exchangeability is a modelling constraint, not a claim that the dyads are truly exchangeable. To test that constraint, fit both models to the same prepared observations.

Exchangeable Dyads - Cross-Sectional APIM: Preparation

# Reuse the distinguishable preparation, including its arbitrary
# -1/+1 member contrast.
df_exchangeable_apim <- df_apim

Exchangeable Dyads - Cross-Sectional APIM: Prepared Data

print_apim_preview(df_exchangeable_apim)
userID coupleID gender .is_male .is_female .member_contrast_arbitrary satisfaction .communication_gmc_actor .communication_gmc_partner
1_1 1 female 0 1 -1 4.841332 -0.1553604 0.9980657
1_2 1 male 1 0 1 5.577165 0.9980657 -0.1553604
2_1 2 female 0 1 -1 7.248741 0.8748537 1.4272563
2_2 2 male 1 0 1 6.960293 1.4272563 0.8748537
3_1 3 female 0 1 -1 6.928699 -0.7224960 0.5095933
3_2 3 male 1 0 1 6.077688 0.5095933 -0.7224960

Exchangeable Dyads - Cross-Sectional APIM: glmmTMB

  • For the fixed effects, we can pool the intercepts .is_male and .is_female by simply reintroducing a single shared intercept 1.

  • For the actor and partner effects, we can remove the interactions to get a pooled estimate.

  • For the residual, we want to constrain both partners’ variances to be equal. Often, this is attempted via: (1 | coupleID).

    • This gives both partners the same dyad-level random intercept variance.
    • With the usual Gaussian residual variance still present, the implied partner correlation is positive but not freely estimated: \(\frac{\tau^2}{\tau^2 + \sigma^2}\).
    • If residual variance is removed with dispformula = ~ 0, the shared random intercept would imply a correlation near 1, which is too restrictive.

Exchangeable Dyads - Cross-Sectional APIM: glmmTMB

One solution in glmmTMB is the following:

model_exch_apim_glmmtmb <- glmmTMB(
  satisfaction ~

    # Pooled single intercept
    1 +

    # Pooled actor and partner effects
    .communication_gmc_actor + .communication_gmc_partner +

    # Independent common and member-deviation terms parameterize a homogeneous
    # two-member covariance using dyadMLM's arbitrary contrast.
    (1 | coupleID) + (0 + .member_contrast_arbitrary | coupleID)

    # The dyadic covariance is represented by the two random-effect blocks, so
    # remove the additional independent Gaussian residual variance.
    , dispformula = ~ 0

    , data = df_exchangeable_apim
    , family = gaussian()
)

Exchangeable APIM - Covariance Recovery

recover_exchangeable_covariance() back-transforms the common and deviation blocks to the homogeneous member variance and partner correlation:

Recovered exchangeable member-level covariance

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

Variance-covariance:
                       1     2    
1 member1: (Intercept) 2.403 0.827
2 member2: (Intercept) 0.827 2.403

Standard deviations and correlations:
                       1     2    
1 member1: (Intercept) 1.550 0.344
2 member2: (Intercept) 0.344 1.550

Exchangeable Dyads - Cross-Sectional APIM: Interpreting the Covariance

Let \(\tau_c^2\) be the variance of the common intercept and \(\tau_d^2\) the variance of the independent member deviation. With .member_contrast_arbitrary = -1/+1, the back-transformation is

\[ \operatorname{Var}(e_1)=\operatorname{Var}(e_2)=\tau_c^2+\tau_d^2, \]

\[ \operatorname{Cov}(e_1,e_2)=\tau_c^2-\tau_d^2, \qquad \rho=\frac{\tau_c^2-\tau_d^2}{\tau_c^2+\tau_d^2}. \]

The model can therefore represent negative partner correlation when \(\tau_d^2>\tau_c^2\) (del Rosario & West, 2025). For repeated measures, index the same structure at the couple-day level.

Exchangeable Dyads - Cross-Sectional APIM: brms

In brms, this issue can be solved more easily with a native residual correlation structure and one pooled residual variance:

formula_exch_apim <- bf(
  satisfaction ~ 1 +
    .communication_actor + .communication_partner +

    # To achieve a homogeneous two-partner covariance structure we use an
    # unstructured residual correlation over the two member positions.
    unstr(time = .member_contrast_arbitrary, gr = coupleID)

  # ... but now, we do not allow individual variances anymore.
  # We want one pooled residual variance. So we do not specify a sigma formula.
  # This implies:
  # , sigma = ~ 1
)

priors_exch_apim <- c(
  prior(normal(2, 10), class = "Intercept"),
  prior(normal(0, 5), class = "b"),
  prior(student_t(3, 0, 1.5), class = "sigma"),
  prior(lkj(2), class = "cortime")
)

model_exch_apim <- brm(
  formula = formula_exch_apim,
  data = df_exchangeable_apim,
  family = gaussian(link = identity),
  prior = priors_exch_apim,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000,
  seed = 123,
  file = file.path('brms_cache', 'example1_ind_apim') # Cache the model
)

Exchangeable Dyads - Cross-Sectional APIM: Results

Est. 95% CI Rhat Bulk_ESS Tail_ESS
Fixed Effects
Intercept -4.71 [-5.13, -4.26] 1.004 4681 2769
.communication_actor 1.72 [ 1.66, 1.78] 1.002 4425 2895
.communication_partner 0.24 [ 0.18, 0.30] 1.000 4971 3242
Residual Structure
cortime(-1,1) 0.34 [0.26, 0.42] 1.000 4665 3140
sigma 1.55 [1.49, 1.63] 1.001 4114 3191

Constraints vs Prediction

compare_nested_models() tests the imposed exchangeability constraints with a likelihood-ratio comparison:

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 Chi Df Pr(>Chisq)    
model_exch_apim_glmmtmb  5 4026.6 4051.7 -2008.3   4016.6                             
model_dist_apim_glmmtmb  9 3776.7 3821.7 -1879.3   3758.7 257.95      4  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

The following observation-level PSIS-LOO comparison instead evaluates conditional predictive performance (Bürkner et al., 2021). It is not a formal test of distinguishability or prediction of a new dyad.

Exchangeable Dyads - APIM: LOO Comparison

Model elpd_diff se_diff
distinguishable 0.0 0.0
exchangeable -172.9 16.5

For this observation-level target, the distinguishable model predicts better. The ELPD difference is descriptive predictive evidence; se_diff summarizes its uncertainty (Sivula et al., 2025; Vehtari et al., 2017).

Exchangeable Dyads - Cross-Sectional DIM

Exchangeable Dyads - Cross-Sectional DIM: Data

For the DIM parameterization, prepare the raw observations under the exchangeability assumption.

userID coupleID satisfaction communication .communication_dyad_mean_gmc .communication_within_dyad_dev .member_contrast_arbitrary
1_1 1 4.841332 4.851107 0.4213527 -0.5767130 -1
1_2 1 5.577165 6.004533 0.4213527 0.5767130 1
2_1 2 7.248741 5.881321 1.1510550 -0.2762013 -1
2_2 2 6.960293 6.433723 1.1510550 0.2762013 1
3_1 3 6.928699 4.283971 -0.1064514 -0.6160447 -1
3_2 3 6.077688 5.516060 -0.1064514 0.6160447 1

Exchangeable Dyads - Cross-Sectional DIM: Centering

Decompose communication variance into:

  • Between-couple (cbc): Couple-mean communication skills in relation to other couples
  • Within-couple (cwc): Individuals’ communication skills in relation to their couple-mean

Same assumption about exchangeability as in the APIM

Exchangeable Dyads - Cross-Sectional DIM: Centering

  • .communication_dyad_mean_gmc: the grand-mean-centered dyad mean
  • .communication_within_dyad_dev: each member’s deviation from that dyad mean

Exchangeable Dyads - Cross-Sectional DIM: Centering

userID coupleID satisfaction .communication_dyad_mean_gmc .communication_within_dyad_dev .member_contrast_arbitrary
1_1 1 4.841332 0.4213527 -0.5767130 -1
1_2 1 5.577165 0.4213527 0.5767130 1
2_1 2 7.248741 1.1510550 -0.2762013 -1
2_2 2 6.960293 1.1510550 0.2762013 1
3_1 3 6.928699 -0.1064514 -0.6160447 -1
3_2 3 6.077688 -0.1064514 0.6160447 1

Exchangeable Dyads - Cross-Sectional DIM: Model

formula_exch_dim <- bf(
  satisfaction ~ 1 +
    .communication_dyad_mean_gmc + .communication_within_dyad_dev +
    unstr(time = .member_contrast_arbitrary, gr = coupleID)

    # sigma ~ 1 is implied.
)

priors_exch_dim <- c(
  prior(normal(2, 10), class = "Intercept"),
  prior(normal(0, 5), class = "b"),
  prior(student_t(3, 0, 1.5), class = "sigma"),
  prior(lkj(2), class = "cortime")
)

model_exch_dim <- brm(
  formula = formula_exch_dim,
  data = df_exchangeable_dim,
  family = gaussian(link = identity),
  prior = priors_exch_dim,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000,
  seed = 123,
  file = file.path('brms_cache', 'example1_ind_dim') # Cache the model
)

Exchangeable Dyads - Cross-Sectional DIM: Results

Est. 95% CI Rhat Bulk_ESS Tail_ESS
Fixed Effects
Intercept 5.11 [5.00, 5.22] 1.001 4390 2769
.communication_dyad_mean_gmc 1.96 [1.87, 2.05] 1.003 4658 3120
.communication_within_dyad_dev 1.48 [1.39, 1.57] 1.002 4285 2998
Residual Structure
cortime(-1,1) 0.34 [0.27, 0.42] 1.001 4164 2873
sigma 1.55 [1.49, 1.62] 1.000 4560 3396

Equivalence APIM and DIM

dyadMLM grand-mean centers the DIM dyad mean, whereas these APIM predictors are raw. APIM and DIM slopes require the sum/difference transformations shown later; their residual covariance is directly comparable. The printed intercepts also use different origins, so shift the APIM intercept to average communication before comparing it (Bolger et al., 2025):

\[ b_{0,\,APIM(GMC)} = b_{0,\,APIM(raw)} + \bar{x} (b_{actor} + b_{partner}). \]

Est. 95% CI
Fixed Effects
Intercept -4.71 [-5.13, -4.26]
.communication_actor 1.72 [ 1.66, 1.78]
.communication_partner 0.24 [ 0.18, 0.30]
Residual Structure
cortime(-1,1) 0.34 [0.26, 0.42]
sigma 1.55 [1.49, 1.63]
Est. 95% CI
Fixed Effects
Intercept 5.11 [5.00, 5.22]
.communication_dyad_mean_gmc 1.96 [1.87, 2.05]
.communication_within_dyad_dev 1.48 [1.39, 1.57]
Residual Structure
cortime(-1,1) 0.34 [0.27, 0.42]
sigma 1.55 [1.49, 1.62]

Equivalence APIM and DIM

Est. 95% CI
Fixed Effects
Intercept -4.71 [-5.13, -4.26]
.communication_actor 1.72 [ 1.66, 1.78]
.communication_partner 0.24 [ 0.18, 0.30]
Est. 95% CI
Fixed Effects
Intercept 5.11 [5.00, 5.22]
.communication_dyad_mean_gmc 1.96 [1.87, 2.05]
.communication_within_dyad_dev 1.48 [1.39, 1.57]

Top sliders: grand-mean-centered Communication for Actor and Partner (APIM).
Bottom sliders: reparameterized communication — Centered Between-Couple (xcbc) and Centered Within-Couple (xcwc) (DIM).

0.00 0.00 0.00 0.00

Equivalence APIM and DIM

Coordinate Inputs

APIM: individual coordinates
0.0
0.0
DIM: mean and deviation coordinates
0.0
0.0
xcbc = (xactor + xpartner) / 2
xcwc = (xactor − xpartner) / 2
APIM
DIM
Coefficients
Between-couple effect: the mean. If both partners increase by 1, the couple mean increases by 1 and the deviation stays 0.
bcbc = bactor + bpartner
Within-couple effect: the deviation. If the actor increases by 1 and the partner decreases by 1, the couple mean stays 0 and the deviation increases by 1.
bcwc = bactor − bpartner
The dot lives in APIM space: horizontal = actor predictor, vertical = partner predictor. The diagonal axes show the same point in DIM space.

Equivalence APIM and DIM

  • If the couple mean goes up by 1 and the within-couple deviation from the couple mean stays fixed, this means that both partners’ predictors must go up by 1 unit. Thus, the effects are linear combinations:

\[b_{cbc} = b_{actor} + b_{partner}\]

  • Similarly, if a person’s deviation from their couple mean is one unit higher and the couple mean stays fixed, this means their partners’ deviation must be one unit lower. Thus, the effects are a linear combination:

\[b_{cwc} = b_{actor} - b_{partner}\]

(Bolger et al., 2025)

Equivalence APIM and DIM

  • Conversely, to obtain the actor and partner effects given the DIM estimates:

\[b_{actor} = \frac{b_{cbc} + b_{cwc}}{2}\]

\[b_{partner} = \frac{b_{cbc} - b_{cwc}}{2}\]

Equivalence APIM and DIM

Example: using APIM coefficients to obtain DIM coefficients and computing Credible Intervals of DIM coefficients.

a <- hypothesis(
  model_exch_apim,
  ".communication_actor + .communication_partner = 0"
)

round(a$hypothesis[,c(2,3,4,5)], 2)
  Estimate Est.Error CI.Lower CI.Upper
1     1.96      0.04     1.87     2.04
Est. 95% CI
Fixed Effects
Intercept -4.71 [-5.13, -4.26]
.communication_actor 1.72 [ 1.66, 1.78]
.communication_partner 0.24 [ 0.18, 0.30]
Est. 95% CI
Fixed Effects
Intercept 5.11 [5.00, 5.22]
.communication_dyad_mean_gmc 1.96 [1.87, 2.05]
.communication_within_dyad_dev 1.48 [1.39, 1.57]

Equivalence APIM and DIM: Takeaway

  • APIM and DIM are likelihood reparametrizations of the same model (aside from prior impact in Bayesian fits)
  • APIM: intuitive actor/partner framing
  • DIM: clean separation of between vs within components
  • Estimating one model allows for directly obtaining estimates of the other
  • Same random-effects structure at dyad level and same assumptions

Side Note:

In the distinguishable case, the full Dyadic Score Model (DSM) (e.g., Stadler et al., 2023) is equivalent to the corresponding full APIM (aside from prior impact in Bayesian fits) (Bolger et al., 2025; Iida et al., 2018).

Exchangeable Dyads - L-DIM/L-APIM: Simulated Data

We start again from the raw longitudinal data and omit role, which tells dyadMLM to treat the two members as exchangeable. Requesting both model types creates the APIM and DIM predictors in the same model-ready object.

df_long_dim <- prepare_dyad_data(
  data = df_long,
  dyad = coupleID,
  member = userID,
  time = diaryday,
  predictors = provided_support,
  model_types = c("apim", "dim"),
  seed = 123
) %>%
  mutate(
    diaryday_c = diaryday - mean(diaryday)
  )
userID coupleID diaryday gender closeness provided_support .composition .composition_role .is_exchangeable .member_contrast_arbitrary .provided_support_cwp .provided_support_cbp .provided_support_actor .provided_support_partner .provided_support_cwp_actor .provided_support_cwp_partner .provided_support_cbp_actor .provided_support_cbp_partner .provided_support_dyad_mean_gmc .provided_support_within_dyad_dev .provided_support_cwp_dyad_mean .provided_support_cwp_within_dyad_dev .provided_support_cbp_dyad_mean .provided_support_cbp_within_dyad_dev diaryday_c
31_1 31 0 female 5.47 6.02 assumed_exchangeable assumed_exchangeable 1 -1 0.41 0.65 6.02 7.09 0.41 1.54 0.65 0.58 1.59 -0.53 0.97 -0.57 0.61 0.04 -27
31_1 31 1 female 7.69 6.13 assumed_exchangeable assumed_exchangeable 1 -1 0.51 0.65 6.13 6.95 0.51 1.4 0.65 0.58 1.57 -0.41 0.96 -0.44 0.61 0.04 -26
31_1 31 2 female 5.83 6.1 assumed_exchangeable assumed_exchangeable 1 -1 0.49 0.65 6.1 6.34 0.49 0.79 0.65 0.58 1.25 -0.12 0.64 -0.15 0.61 0.04 -25
31_1 31 3 female 6.55 6.99 assumed_exchangeable assumed_exchangeable 1 -1 1.38 0.65 6.99 5.39 1.38 -0.15 0.65 0.58 1.23 0.8 0.61 0.77 0.61 0.04 -24
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
31_2 31 0 male 5.19 7.09 assumed_exchangeable assumed_exchangeable 1 1 1.54 0.58 7.09 6.02 1.54 0.41 0.58 0.65 1.59 0.53 0.97 0.57 0.61 -0.04 -27
31_2 31 1 male 6.96 6.95 assumed_exchangeable assumed_exchangeable 1 1 1.4 0.58 6.95 6.13 1.4 0.51 0.58 0.65 1.57 0.41 0.96 0.44 0.61 -0.04 -26
31_2 31 2 male 6.9 6.34 assumed_exchangeable assumed_exchangeable 1 1 0.79 0.58 6.34 6.1 0.79 0.49 0.58 0.65 1.25 0.12 0.64 0.15 0.61 -0.04 -25
31_2 31 3 male 7.59 5.39 assumed_exchangeable assumed_exchangeable 1 1 -0.15 0.58 5.39 6.99 -0.15 1.38 0.58 0.65 1.23 -0.8 0.61 -0.77 0.61 -0.04 -24

Exchangeable Dyads - Intensive Longitudinal

Exchangeable Dyads - Longitudinal DIM/APIM

  • We need an APIM or DIM on both the between person level and the within-person level.
  • Due to equivalence, we could have a DIM on one level and an APIM on the other.

The generated .provided_support_cwp_dyad_mean and .provided_support_cwp_within_dyad_dev columns form the within-person DIM. The corresponding .provided_support_cbp_dyad_mean and .provided_support_cbp_within_dyad_dev columns form the between-person DIM (Bolger & Laurenceau, 2013).

Exchangeable Dyads - L-DIM/L-APIM: Preparing the Member Contrast

dyadMLM generates .member_contrast_arbitrary as an arbitrary -1/+1 member contrast. The fixed seed makes that assignment reproducible.

Exchangeable Dyads - L-DIM/L-APIM: Preparing the Member Contrast

userID coupleID diaryday gender closeness provided_support .composition .composition_role .is_exchangeable .member_contrast_arbitrary .provided_support_cwp .provided_support_cbp .provided_support_actor .provided_support_partner .provided_support_cwp_actor .provided_support_cwp_partner .provided_support_cbp_actor .provided_support_cbp_partner .provided_support_dyad_mean_gmc .provided_support_within_dyad_dev .provided_support_cwp_dyad_mean .provided_support_cwp_within_dyad_dev .provided_support_cbp_dyad_mean .provided_support_cbp_within_dyad_dev diaryday_c
31_1 31 0 female 5.47 6.02 assumed_exchangeable assumed_exchangeable 1 -1 0.41 0.65 6.02 7.09 0.41 1.54 0.65 0.58 1.59 -0.53 0.97 -0.57 0.61 0.04 -27
31_1 31 1 female 7.69 6.13 assumed_exchangeable assumed_exchangeable 1 -1 0.51 0.65 6.13 6.95 0.51 1.4 0.65 0.58 1.57 -0.41 0.96 -0.44 0.61 0.04 -26
31_1 31 2 female 5.83 6.1 assumed_exchangeable assumed_exchangeable 1 -1 0.49 0.65 6.1 6.34 0.49 0.79 0.65 0.58 1.25 -0.12 0.64 -0.15 0.61 0.04 -25
31_1 31 3 female 6.55 6.99 assumed_exchangeable assumed_exchangeable 1 -1 1.38 0.65 6.99 5.39 1.38 -0.15 0.65 0.58 1.23 0.8 0.61 0.77 0.61 0.04 -24
31_1 31 4 female 7.03 6.7 assumed_exchangeable assumed_exchangeable 1 -1 1.08 0.65 6.7 5.19 1.08 -0.36 0.65 0.58 0.97 0.75 0.36 0.72 0.61 0.04 -23
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
31_2 31 0 male 5.19 7.09 assumed_exchangeable assumed_exchangeable 1 1 1.54 0.58 7.09 6.02 1.54 0.41 0.58 0.65 1.59 0.53 0.97 0.57 0.61 -0.04 -27
31_2 31 1 male 6.96 6.95 assumed_exchangeable assumed_exchangeable 1 1 1.4 0.58 6.95 6.13 1.4 0.51 0.58 0.65 1.57 0.41 0.96 0.44 0.61 -0.04 -26
31_2 31 2 male 6.9 6.34 assumed_exchangeable assumed_exchangeable 1 1 0.79 0.58 6.34 6.1 0.79 0.49 0.58 0.65 1.25 0.12 0.64 0.15 0.61 -0.04 -25
31_2 31 3 male 7.59 5.39 assumed_exchangeable assumed_exchangeable 1 1 -0.15 0.58 5.39 6.99 -0.15 1.38 0.58 0.65 1.23 -0.8 0.61 -0.77 0.61 -0.04 -24
31_2 31 4 male 7.55 5.19 assumed_exchangeable assumed_exchangeable 1 1 -0.36 0.58 5.19 6.7 -0.36 1.08 0.58 0.65 0.97 -0.75 0.36 -0.72 0.61 -0.04 -23

Exchangeable Dyads - L-DIM/L-APIM: glmmTMB

model_exch_long_apim_glmmtmb <- glmmTMB(
  closeness ~ 1 +

    diaryday_c +

    # Within-person APIM
    .provided_support_cwp_actor +
    .provided_support_cwp_partner +
    # Equivalent to within-person DIM:
    # .provided_support_cwp_dyad_mean +
    # .provided_support_cwp_within_dyad_dev +

    # Between-person APIM
    .provided_support_cbp_actor +
    .provided_support_cbp_partner +
    # Equivalent to between-person DIM
    # .provided_support_cbp_within_dyad_dev +
    # .provided_support_cbp_dyad_mean +

    # Dyad-Level intercept and slopes for time-varying predictors
    (1 + diaryday_c + .provided_support_cwp_actor + .provided_support_cwp_partner | coupleID) +

    # Both partners' deviations from these dyad-level means and slopes.
    # Note: separate random-effects block to make them uncorrelated with dyad-level REs.
    (0 + .member_contrast_arbitrary +
         I(.member_contrast_arbitrary * diaryday_c) +
         I(.member_contrast_arbitrary * .provided_support_cwp_actor) +
         I(.member_contrast_arbitrary * .provided_support_cwp_partner) | coupleID) +

    # After accounting for stable differences between partners, they may still be
    # non-independent each day due to external occurrences and daily exchange processes
    # not captured in the fixed effects. Thus, residuals are not yet independent.
    # Homogeneous same-day residual covariance between partners:
    (1 | coupleID:diaryday) +
      (0 + .member_contrast_arbitrary | coupleID:diaryday)

  # For this Gaussian model, remove the additional independent residual
  # variance because the two couple-day blocks represent it.
  , dispformula = ~ 0
  , data = df_long_dim
)

Exchangeable Dyads - L-DIM/L-APIM: glmmTMB

In this example, the glmmTMB model does not converge properly. The random-effects structure is too complex for the information available in these data, so we would have to simplify it. Instead, we use brms below.

Exchangeable Dyads - L-DIM/L-APIM: brms

In brms, the same model can be estimated as:

Unlike the distinguishable longitudinal model, this exchangeable specification keeps one pooled residual sigma (sigma ~ 1 is implied). Its residual SD prior therefore targets class = "sigma" directly.

formula_exch_long_apim <- bf(
  closeness ~ 1 +

    diaryday_c +

    # Within-person APIM
    .provided_support_cwp_actor +
    .provided_support_cwp_partner +

    # Between-person APIM
    .provided_support_cbp_actor +
    .provided_support_cbp_partner +

    # Dyad-Level intercept and slopes for time-varying predictors
    (1 + diaryday_c + .provided_support_cwp_actor + .provided_support_cwp_partner | coupleID) +

    # Both partners' deviations from these dyad-level means and slopes.
    # Note: separate random-effects block to make them uncorrelated with dyad-level REs.
    (0 + .member_contrast_arbitrary +
         I(.member_contrast_arbitrary * diaryday_c) +
         I(.member_contrast_arbitrary * .provided_support_cwp_actor) +
         I(.member_contrast_arbitrary * .provided_support_cwp_partner) | coupleID) +

    # Modelling Level 1 residuals in brms:
    unstr(time = .member_contrast_arbitrary, gr = coupleID:diaryday)

  # Again, no need to model heterogeneous residual variances (sigma)
  # Implied: sigma = ~ 1
)

priors_exch_long_apim <- c(
  prior(normal(4, 2), class = "Intercept"),
  prior(normal(0, 2), class = "b"),
  prior(exponential(1), class = "sd"),
  prior(lkj(2), class = "cor"),
  prior(student_t(3, 0, 1.5), class = "sigma"),
  prior(lkj(2), class = "cortime")
)

model_exch_long_apim <- brm(
  formula = formula_exch_long_apim,
  data = df_long_dim,
  family = gaussian(link = identity),
  prior = priors_exch_long_apim,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000,
  seed = 123,
  file = file.path('brms_cache', 'model_apim_ind_long_apim')
)

Exchangeable Dyads - L-DIM/L-APIM: brms

The equivalent DIM specification:

# DIM Version
formula_exch_long_dim <- bf(
  closeness ~ 1 +

    diaryday_c +

    .provided_support_cwp_dyad_mean +
    .provided_support_cwp_within_dyad_dev +

    .provided_support_cbp_dyad_mean +
    .provided_support_cbp_within_dyad_dev +

    (1 + diaryday_c + .provided_support_cwp_actor + .provided_support_cwp_partner | coupleID) +
    (0 + .member_contrast_arbitrary +
         I(.member_contrast_arbitrary * diaryday_c) +
         I(.member_contrast_arbitrary * .provided_support_cwp_actor) +
         I(.member_contrast_arbitrary * .provided_support_cwp_partner) | coupleID) +

    unstr(time = .member_contrast_arbitrary, gr = coupleID:diaryday)
)

priors_exch_long_dim <- c(
  prior(normal(4, 2), class = "Intercept"),
  prior(normal(0, 2), class = "b"),
  prior(exponential(1), class = "sd"),
  prior(lkj(2), class = "cor"),
  prior(student_t(3, 0, 1.5), class = "sigma"),
  prior(lkj(2), class = "cortime")
)

model_exch_long_dim <- brm(
  formula = formula_exch_long_dim,
  data = df_long_dim,
  family = gaussian(link = identity),
  prior = priors_exch_long_dim,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000,
  seed = 123,
  file = file.path('brms_cache', 'model_apim_ind_long_dim')
)

Exchangeable Dyads - L-DIM/L-APIM: Member-Contrast Coding

For appropriate exchangeable random effects, we can use the generated .member_contrast_arbitrary column (del Rosario & West, 2025; Kenny & Ackerman, 2023):

  • Randomly give one partner a constant of -1 and the other a constant of 1
  • The couple-level intercept represents the average/common intercept of the two partners. (1 | coupleID)
  • A column of 1s and -1s represents each partner’s deviation from the couple-level effect, with partners contributing equally in opposite directions (+1 and -1)

Exchangeability is retained: if partners are flipped, the result is the same.

(del Rosario & West (2025) provide practical guidance on how to reduce the random effects structure in case of non-convergence.)

Exchangeable Dyads - L-DIM/L-APIM: Common and Deviation Components

For two partners coded with .member_contrast_arbitrary = +1 and .member_contrast_arbitrary = -1, the glmmTMB structure (1 | coupleID:diaryday) + (0 + .member_contrast_arbitrary | coupleID:diaryday) uses independent common and member-deviation components to imply a homogeneous two-partner residual covariance.

This requires exactly two rows per couple-day and independent common and deviation blocks, and dispformula = ~ 0. If their variances are \(\tau_c^2\) and \(\tau_d^2\), member variance is \(\tau_c^2+\tau_d^2\), covariance is \(\tau_c^2-\tau_d^2\), and their ratio is the correlation. No manual partner indicators are needed (del Rosario & West, 2025).

Exchangeable Dyads - L-DIM/L-APIM: Results

Comparing APIM (left) and DIM (right) output of fixed effects.

Est. 95% CI
Fixed Effects
Intercept 5.10 [4.97, 5.24]
diaryday_c 0.00 [0.00, 0.01]
.provided_support_cwp_actor 0.30 [0.27, 0.33]
.provided_support_cwp_partner 0.15 [0.12, 0.18]
.provided_support_cbp_actor 1.32 [1.14, 1.48]
.provided_support_cbp_partner 0.38 [0.21, 0.56]
Est. 95% CI
Fixed Effects
Intercept 5.11 [4.98, 5.23]
diaryday_c 0.00 [0.00, 0.01]
.provided_support_cwp_dyad_mean 0.45 [0.41, 0.49]
.provided_support_cwp_within_dyad_dev 0.15 [0.11, 0.19]
.provided_support_cbp_dyad_mean 1.70 [1.55, 1.86]
.provided_support_cbp_within_dyad_dev 0.93 [0.61, 1.21]
  • Equivalence holds on both levels. Numeric deviations arise due to priors, sampling noise, and rounding.

Exchangeable Dyads - L-DIM/L-APIM: Results

APIM Results in Detail:

Est. 95% CI Rhat Bulk_ESS Tail_ESS
Fixed Effects
Intercept 5.10 [4.97, 5.24] 1.017 175 225
diaryday_c 0.00 [0.00, 0.01] 1.005 463 1553
.provided_support_cwp_actor 0.30 [0.27, 0.33] 1.001 2845 2572
.provided_support_cwp_partner 0.15 [0.12, 0.18] 1.000 3249 2819
.provided_support_cbp_actor 1.32 [1.14, 1.48] 1.017 171 286
.provided_support_cbp_partner 0.38 [0.21, 0.56] 1.005 212 492
Random Effects
sd(Intercept) 0.67 [ 0.58, 0.77] 1.022 389 702
sd(diaryday_c) 0.01 [ 0.01, 0.01] 1.001 1581 2780
sd(.provided_support_cwp_actor) 0.05 [ 0.01, 0.09] 1.004 990 953
sd(.provided_support_cwp_partner) 0.08 [ 0.04, 0.12] 1.002 1158 1197
sd(.member_contrast_arbitrary) 0.84 [ 0.74, 0.97] 1.049 218 539
sd(I.member_contrast_arbitraryMUdiaryday_c) 0.01 [ 0.01, 0.02] 1.010 1152 1947
sd(I.member_contrast_arbitraryMU.provided_support_cwp_actor) 0.14 [ 0.11, 0.18] 1.004 1154 2284
sd(I.member_contrast_arbitraryMU.provided_support_cwp_partner) 0.07 [ 0.01, 0.11] 1.004 754 627
cor(Intercept,diaryday_c) 0.56 [ 0.38, 0.70] 1.004 1129 1909
cor(Intercept,.provided_support_cwp_actor) 0.49 [ 0.02, 0.82] 1.002 2958 2405
cor(diaryday_c,.provided_support_cwp_actor) 0.37 [-0.10, 0.76] 1.002 2789 2668
cor(Intercept,.provided_support_cwp_partner) 0.21 [-0.12, 0.53] 1.001 2870 2182
cor(diaryday_c,.provided_support_cwp_partner) 0.19 [-0.19, 0.54] 1.003 2961 2724
cor(.provided_support_cwp_actor,.provided_support_cwp_partner) 0.35 [-0.27, 0.80] 1.002 687 1343
cor(.member_contrast_arbitrary,I.member_contrast_arbitraryMUdiaryday_c) 0.60 [ 0.44, 0.73] 1.014 788 1455
cor(.member_contrast_arbitrary,I.member_contrast_arbitraryMU.provided_support_cwp_actor) 0.66 [ 0.46, 0.81] 1.004 1299 1965
cor(I.member_contrast_arbitraryMUdiaryday_c,I.member_contrast_arbitraryMU.provided_support_cwp_actor) 0.42 [ 0.15, 0.63] 1.008 1619 2359
cor(.member_contrast_arbitrary,I.member_contrast_arbitraryMU.provided_support_cwp_partner) 0.12 [-0.29, 0.53] 1.000 2816 1776
cor(I.member_contrast_arbitraryMUdiaryday_c,I.member_contrast_arbitraryMU.provided_support_cwp_partner) 0.11 [-0.32, 0.55] 1.003 3148 1966
cor(I.member_contrast_arbitraryMU.provided_support_cwp_actor,I.member_contrast_arbitraryMU.provided_support_cwp_partner) 0.34 [-0.17, 0.76] 1.000 1647 2123
cortime(-1,1) 0.04 [ 0.01, 0.07] 1.001 4908 2263
Residual Structure
sigma 0.94 [ 0.92, 0.95] 1.000 6231 2668

Extract Full APIM Random Effects Variance-Covariance Matrix

We can rotate the random effects structure back to obtain the full APIM Random effects variance-covariance matrix, based on (Kenny & Ackerman, 2023). Just as we would in a SEM model with equality constraints!

For example, the within-partner variance of the random intercept is:

\[ \operatorname{Var}(u_{\text{common intercept}}) + \operatorname{Var}(u_{\text{deviation intercept}}) \]

and the cross-partner covariance of both partners’ random intercepts is:

\[ \operatorname{Var}(u_{\text{common intercept}}) - \operatorname{Var}(u_{\text{deviation intercept}}) \]

In the model output, the common-intercept variance is labelled Intercept, and the deviation-intercept variance is labelled .member_contrast_arbitrary. This rotation only holds when the common/couple-mean and member-deviation random effects are uncorrelated. This is why we put them in separate random-effects blocks.

Extract Full APIM Random Effects Variance-Covariance Matrix

a <- summarize_exchangeable_apim_brms(
  model_exch_long_apim,
  gr = "coupleID",
  deviation_term = ".member_contrast_arbitrary",
  partner_labels = c("PartnerA", "PartnerB")
)

Extract Full APIM Random Effects Variance-Covariance Matrix: Within-Person Matrix

print_df(a$within_person_covariance_matrix)
Intercept diaryday_c .provided_support_cwp_actor .provided_support_cwp_partner
Intercept 1.1802004 0.0112332 0.0955767 0.0179927
diaryday_c 0.0112332 0.0003159 0.0010075 0.0002727
.provided_support_cwp_actor 0.0955767 0.0010075 0.0238464 0.0042531
.provided_support_cwp_partner 0.0179927 0.0002727 0.0042531 0.0118398

Extract Full APIM Random Effects Variance-Covariance Matrix: Cross-Person Matrix

print_df(a$cross_person_covariance_matrix)
Intercept diaryday_c .provided_support_cwp_actor .provided_support_cwp_partner
Intercept -0.2635988 -0.0024112 -0.0630951 0.0048753
diaryday_c -0.0024112 -0.0000434 -0.0005769 0.0000857
.provided_support_cwp_actor -0.0630951 -0.0005769 -0.0179709 -0.0015373
.provided_support_cwp_partner 0.0048753 0.0000857 -0.0015373 0.0024947

Extract Full APIM Random Effects Variance-Covariance Matrix: Full Matrix

print_wide_matrix(a$full_covariance_matrix, font_size = 6)
PartnerA_Intercept PartnerA_diaryday_c PartnerA_.provided_support_cwp_actor PartnerA_.provided_support_cwp_partner PartnerB_Intercept PartnerB_diaryday_c PartnerB_.provided_support_cwp_actor PartnerB_.provided_support_cwp_partner
PartnerA_Intercept 1.180 0.011 0.096 0.018 -0.264 -0.002 -0.063 0.005
PartnerA_diaryday_c 0.011 0.000 0.001 0.000 -0.002 0.000 -0.001 0.000
PartnerA_.provided_support_cwp_actor 0.096 0.001 0.024 0.004 -0.063 -0.001 -0.018 -0.002
PartnerA_.provided_support_cwp_partner 0.018 0.000 0.004 0.012 0.005 0.000 -0.002 0.002
PartnerB_Intercept -0.264 -0.002 -0.063 0.005 1.180 0.011 0.096 0.018
PartnerB_diaryday_c -0.002 0.000 -0.001 0.000 0.011 0.000 0.001 0.000
PartnerB_.provided_support_cwp_actor -0.063 -0.001 -0.018 -0.002 0.096 0.001 0.024 0.004
PartnerB_.provided_support_cwp_partner 0.005 0.000 -0.002 0.002 0.018 0.000 0.004 0.012

The helper function summarize_exchangeable_apim_brms() is defined in 00_R_Functions/ReportModels.R. It rotates every posterior draw, so the interval estimates for the transformed SDs, correlations, covariances, and residual quantities are computed on the transformed scale.

For glmmTMB, the same back-transformation is algebraically possible, but interval estimation is less direct because there are no posterior draws. One would usually transform parametric bootstrap draws or asymptotic draws from the fitted parameter covariance matrix.

Random effects SDs

print_df(
  a$sd_summary %>%
    mutate(across(where(is.numeric), ~ round(.x, 3)))
)
parameter estimate est_error q_low q_high
PartnerA_.provided_support_cwp_actor 0.153 0.018 0.119 0.191
PartnerA_.provided_support_cwp_partner 0.107 0.020 0.066 0.144
PartnerA_Intercept 1.085 0.056 0.987 1.206
PartnerA_diaryday_c 0.018 0.001 0.016 0.020
PartnerB_.provided_support_cwp_actor 0.153 0.018 0.119 0.191
PartnerB_.provided_support_cwp_partner 0.107 0.020 0.066 0.144
PartnerB_Intercept 1.085 0.056 0.987 1.206
PartnerB_diaryday_c 0.018 0.001 0.016 0.020

Random effect correlation matrix

print_wide_matrix(a$full_correlation_matrix, font_size = 6)
PartnerA_Intercept PartnerA_diaryday_c PartnerA_.provided_support_cwp_actor PartnerA_.provided_support_cwp_partner PartnerB_Intercept PartnerB_diaryday_c PartnerB_.provided_support_cwp_actor PartnerB_.provided_support_cwp_partner
PartnerA_Intercept 1.000 0.582 0.570 0.152 -0.223 -0.125 -0.376 0.041
PartnerA_diaryday_c 0.582 1.000 0.367 0.141 -0.125 -0.137 -0.210 0.044
PartnerA_.provided_support_cwp_actor 0.570 0.367 1.000 0.253 -0.376 -0.210 -0.754 -0.091
PartnerA_.provided_support_cwp_partner 0.152 0.141 0.253 1.000 0.041 0.044 -0.091 0.211
PartnerB_Intercept -0.223 -0.125 -0.376 0.041 1.000 0.582 0.570 0.152
PartnerB_diaryday_c -0.125 -0.137 -0.210 0.044 0.582 1.000 0.367 0.141
PartnerB_.provided_support_cwp_actor -0.376 -0.210 -0.754 -0.091 0.570 0.367 1.000 0.253
PartnerB_.provided_support_cwp_partner 0.041 0.044 -0.091 0.211 0.152 0.141 0.253 1.000

Full reporting summary

section_runs <- rle(a$reporting_summary$section)
section_ends <- cumsum(section_runs$lengths)
section_starts <- section_ends - section_runs$lengths + 1
section_labels <- c(
  "Fixed effects" = "Fixed Effects",
  "Back-transformed random-effect SDs" = "Back-Transformed SDs",
  "Back-transformed random-effect correlations" = "Back-Transformed Correlations",
  "Residual structure" = "Residual Structure"
)
reporting_rows_to_pack <- Map(c, section_starts, section_ends)
names(reporting_rows_to_pack) <- unname(section_labels[section_runs$values])

print_df(
  a$reporting_summary %>%
    mutate(across(where(is.numeric), ~ round(.x, 3))),
  rows_to_pack = reporting_rows_to_pack,
  scroll_height = '50vh'
)
section parameter estimate est_error q_low q_high
Fixed Effects
Fixed effects Intercept 5.102 0.070 4.966 5.240
Fixed effects diaryday_c 0.004 0.001 0.001 0.006
Fixed effects .provided_support_cwp_actor 0.300 0.014 0.273 0.328
Fixed effects .provided_support_cwp_partner 0.149 0.015 0.120 0.177
Fixed effects .provided_support_cbp_actor 1.316 0.088 1.140 1.483
Fixed effects .provided_support_cbp_partner 0.382 0.086 0.212 0.557
Back-Transformed SDs
Back-transformed random-effect SDs PartnerA_.provided_support_cwp_actor 0.153 0.018 0.119 0.191
Back-transformed random-effect SDs PartnerA_.provided_support_cwp_partner 0.107 0.020 0.066 0.144
Back-transformed random-effect SDs PartnerA_Intercept 1.085 0.056 0.987 1.206
Back-transformed random-effect SDs PartnerA_diaryday_c 0.018 0.001 0.016 0.020
Back-transformed random-effect SDs PartnerB_.provided_support_cwp_actor 0.153 0.018 0.119 0.191
Back-transformed random-effect SDs PartnerB_.provided_support_cwp_partner 0.107 0.020 0.066 0.144
Back-transformed random-effect SDs PartnerB_Intercept 1.085 0.056 0.987 1.206
Back-transformed random-effect SDs PartnerB_diaryday_c 0.018 0.001 0.016 0.020
Back-Transformed Correlations
Back-transformed random-effect correlations PartnerA_.provided_support_cwp_actor with PartnerA_Intercept 0.572 0.083 0.403 0.723
Back-transformed random-effect correlations PartnerA_.provided_support_cwp_actor with PartnerA_diaryday_c 0.369 0.101 0.169 0.558
Back-transformed random-effect correlations PartnerA_.provided_support_cwp_partner with PartnerA_.provided_support_cwp_actor 0.266 0.152 -0.028 0.558
Back-transformed random-effect correlations PartnerA_.provided_support_cwp_partner with PartnerA_Intercept 0.157 0.117 -0.069 0.382
Back-transformed random-effect correlations PartnerA_.provided_support_cwp_partner with PartnerA_diaryday_c 0.146 0.129 -0.114 0.400
Back-transformed random-effect correlations PartnerA_diaryday_c with PartnerA_Intercept 0.581 0.056 0.470 0.681
Back-transformed random-effect correlations PartnerB_.provided_support_cwp_actor with PartnerA_.provided_support_cwp_actor -0.755 0.151 -0.983 -0.412
Back-transformed random-effect correlations PartnerB_.provided_support_cwp_actor with PartnerA_.provided_support_cwp_partner -0.096 0.160 -0.412 0.215
Back-transformed random-effect correlations PartnerB_.provided_support_cwp_actor with PartnerA_Intercept -0.377 0.103 -0.572 -0.163
Back-transformed random-effect correlations PartnerB_.provided_support_cwp_actor with PartnerA_diaryday_c -0.211 0.111 -0.415 0.019
Back-transformed random-effect correlations PartnerB_.provided_support_cwp_actor with PartnerB_Intercept 0.572 0.083 0.403 0.723
Back-transformed random-effect correlations PartnerB_.provided_support_cwp_actor with PartnerB_diaryday_c 0.369 0.101 0.169 0.558
Back-transformed random-effect correlations PartnerB_.provided_support_cwp_partner with PartnerA_.provided_support_cwp_actor -0.096 0.160 -0.412 0.215
Back-transformed random-effect correlations PartnerB_.provided_support_cwp_partner with PartnerA_.provided_support_cwp_partner 0.230 0.394 -0.570 0.960
Back-transformed random-effect correlations PartnerB_.provided_support_cwp_partner with PartnerA_Intercept 0.043 0.125 -0.199 0.281
Back-transformed random-effect correlations PartnerB_.provided_support_cwp_partner with PartnerA_diaryday_c 0.046 0.138 -0.231 0.310
Back-transformed random-effect correlations PartnerB_.provided_support_cwp_partner with PartnerB_.provided_support_cwp_actor 0.266 0.152 -0.028 0.558
Back-transformed random-effect correlations PartnerB_.provided_support_cwp_partner with PartnerB_Intercept 0.157 0.117 -0.069 0.382
Back-transformed random-effect correlations PartnerB_.provided_support_cwp_partner with PartnerB_diaryday_c 0.146 0.129 -0.114 0.400
Back-transformed random-effect correlations PartnerB_Intercept with PartnerA_.provided_support_cwp_actor -0.377 0.103 -0.572 -0.163
Back-transformed random-effect correlations PartnerB_Intercept with PartnerA_.provided_support_cwp_partner 0.043 0.125 -0.199 0.281
Back-transformed random-effect correlations PartnerB_Intercept with PartnerA_Intercept -0.221 0.095 -0.402 -0.039
Back-transformed random-effect correlations PartnerB_Intercept with PartnerA_diaryday_c -0.124 0.084 -0.288 0.040
Back-transformed random-effect correlations PartnerB_diaryday_c with PartnerA_.provided_support_cwp_actor -0.211 0.111 -0.415 0.019
Back-transformed random-effect correlations PartnerB_diaryday_c with PartnerA_.provided_support_cwp_partner 0.046 0.138 -0.231 0.310
Back-transformed random-effect correlations PartnerB_diaryday_c with PartnerA_Intercept -0.124 0.084 -0.288 0.040
Back-transformed random-effect correlations PartnerB_diaryday_c with PartnerA_diaryday_c -0.136 0.118 -0.362 0.094
Back-transformed random-effect correlations PartnerB_diaryday_c with PartnerB_Intercept 0.581 0.056 0.470 0.681
Residual Structure
Residual structure residual_variance 0.878 0.012 0.854 0.903
Residual structure same_day_residual_correlation 0.037 0.014 0.009 0.066
Residual structure same_day_residual_covariance 0.033 0.012 0.008 0.058
Residual structure sigma 0.937 0.007 0.924 0.950

References

Barr, D. J., Levy, R., Scheepers, C., & Tily, H. J. (2013). Random effects structure for confirmatory hypothesis testing: Keep it maximal. Journal of Memory and Language, 68(3), 255–278. https://doi.org/10.1016/j.jml.2012.11.001
Bolger, N., & Laurenceau, J.-P. (2013). Intensive longitudinal methods: An introduction to diary and experience sampling research. Guilford Press. https://www.guilford.com/books/Intensive-Longitudinal-Methods/Bolger-Laurenceau/9781462506781
Bolger, N., Laurenceau, J.-P., & DiGiovanni, A. (2025). Unified analysis model for indistinguishable and distinguishable dyads. Innovations in Interpersonal Relationships and Health Research: Advancing the Integration of Interdisciplinary Approaches to Dyadic Behavior Change. https://doi.org/10.17605/OSF.IO/WYDCJ
Brooks, M., Bolker, B., Kristensen, K., Maechler, M., Magnusson, A., Skaug, H., Nielsen, A., Berg, C., & van Bentham, K. (2026). glmmTMB: Generalized linear mixed models using template model builder. https://github.com/glmmTMB/glmmTMB
Bürkner, P.-C. (2025). Brms: Bayesian regression models using stan. https://github.com/paul-buerkner/brms
Bürkner, P.-C., Gabry, J., & Vehtari, A. (2021). Efficient leave-one-out cross-validation for bayesian non-factorized normal and student-t models. Computational Statistics, 36, 1243–1261. https://doi.org/10.1007/s00180-020-01045-4
del Rosario, K. S., & West, T. V. (2025). A Practical Guide to Specifying Random Effects in Longitudinal Dyadic Multilevel Modeling. Advances in Methods and Practices in Psychological Science, 8(3), 25152459251351286. https://doi.org/10.1177/25152459251351286
Gabry, J., Češnovar, R., Johnson, A., & Bronder, S. (2025). Cmdstanr: R interface to CmdStan. https://mc-stan.org/cmdstanr/
Guo, J., Gabry, J., Goodrich, B., Johnson, A., Weber, S., & Badr, H. S. (2025). Rstan: R interface to stan. https://mc-stan.org/rstan/
Hartig, F. (2026). DHARMa: Residual diagnostics for hierarchical (multi-level / mixed) regression models. http://florianhartig.github.io/DHARMa/
Iida, M., Seidman, G., & Shrout, P. E. (2018). Models of interdependent individuals versus dyadic processes in relationship research. Journal of Social and Personal Relationships, 35(1), 59–88. https://doi.org/10.1177/0265407517725407
Kenny, D. A., & Ackerman, R. A. (2023). Estimation of random effects with over-time dyadic data using multilevel modeling: The sum and difference method. OSF Preprints. https://doi.org/10.31219/osf.io/fju72
Kenny, D. A., & Cook, W. (1999). Partner effects in relationship research: Conceptual issues, analytic difficulties, and illustrations. Personal Relationships, 6(4), 433–448. https://doi.org/10.1111/j.1475-6811.1999.tb00202.x
Kenny, D. A., & Kashy, D. A. (2011). Dyadic Data Analysis Using Multilevel Modeling. In Handbook of Advanced Multilevel Analysis. Routledge.
Kenny, D. A., Kashy, D. A., & Cook, W. L. (2006). Dyadic data analysis. Guilford Press.
Küng, P. (2026). dyadMLM: Tools for dyadic multilevel models. https://doi.org/10.5281/zenodo.22047083
Lüdecke, D., Makowski, D., Ben-Shachar, M. S., Patil, I., Wiernik, B. M., Bacher, E., & Thériault, R. (2026). Easystats: Framework for easy statistical modeling, visualization, and reporting. https://easystats.github.io/easystats/
Nickell, S. (1981). Biases in dynamic models with fixed effects. Econometrica, 49(6), 1417–1426. https://doi.org/10.2307/1911408
Pinheiro, J., Bates, D., & R Core Team. (2026). Nlme: Linear and nonlinear mixed effects models. https://svn.r-project.org/R-packages/trunk/nlme/
Ripley, B., & Venables, B. (2025). MASS: Support functions and datasets for venables and ripley’s MASS. http://www.stats.ox.ac.uk/pub/MASS4/
Sivula, T., Magnusson, M., Matamoros, A. A., & Vehtari, A. (2025). Uncertainty in bayesian leave-one-out cross-validation based model comparison. Bayesian Analysis.
Stadler, G., Scholz, U., Bolger, N., Shrout, P. E., Knoll, N., & Lüscher, J. (2023). How is companionship related to romantic partners’ affect, relationship satisfaction, and health behavior? Using a longitudinal dyadic score model to understand daily and couple-level effects of a dyadic predictor. Applied Psychology: Health and Well-Being, 15(4), 1530–1554. https://doi.org/10.1111/aphw.12450
Vehtari, A., Gabry, J., Magnusson, M., Yao, Y., Bürkner, P.-C., Paananen, T., & Gelman, A. (2025). Loo: Efficient leave-one-out cross-validation and WAIC for bayesian models. https://mc-stan.org/loo/
Vehtari, A., Gelman, A., & Gabry, J. (2017). Practical bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing, 27, 1413–1432. https://doi.org/10.1007/s11222-016-9696-4
Vuorre, M. (2023). Bmlm: Bayesian multilevel mediation. https://github.com/mvuorre/bmlm/
Wickham, H. (2023). Tidyverse: Easily install and load the tidyverse. https://tidyverse.tidyverse.org