exchangeable-dyads-quick-guide.Rmd

---
title: "Exchangeable Dyads: Model Specification and Back-Transformation"
author: "Pascal Küng"
bibliography: Embed/Temp.bib
csl: Embed/apa.csl
format:
  pdf:
    pdf-engine: xelatex
    toc: false
    number-sections: true
    papersize: a4
    geometry: ["margin=1in"]
    colorlinks: true
    include-in-header:
      text: |
        \usepackage{fvextra}
        \DefineVerbatimEnvironment{Highlighting}{Verbatim}{fontsize=\scriptsize,breaklines,breakanywhere,commandchars=\\\{\}}
        \usepackage{pdflscape}
execute:
  warning: false
  message: false
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE,
  cache = FALSE,
  message = FALSE,
  warning = FALSE
)

library(brms)
library(glmmTMB)
library(dyadMLM)
library(dplyr)
library(tidyr)
library(tibble)
library(knitr)
library(kableExtra)
source(file.path("00_R_Functions", "ReportModels.R"))

options(
  brms.backend = "cmdstanr",
  brms.file_refit = "on_change"
)
```

> **Full GitHub repository:**
> The full repository contains the tutorial source files, reusable helper
> functions, model code, rendered PDFs, and additional information:
> <https://github.com/Pascal-Kueng/05DyadicDataAnalysis>

# Background

Exchangeable dyads do not have substantively meaningful partner roles. Partner
labels are arbitrary, and flipping the labels should not change the model
interpretation.

The full intensive-longitudinal specification used here models
[@bolgerIntensiveLongitudinalMethods2013; @delrosarioPracticalGuideSpecifying2025]:

- pooled fixed effects, because the two partner labels are not substantively
  different;
- pooled average time trends and pooled actor/partner effects;
- stable between-dyad differences in average levels and slopes;
- partner deviations from those dyad-level effects, because the two people in a
  dyad are still distinct repeated-observation units;
- correlations among the dyad-level random effects and among the partner
  deviation random effects;
- same-day residual interdependence between partners after fixed and random
  effects have been accounted for;
- a residual structure that remains invariant when the arbitrary partner labels
  are switched.

```{r exchangeable-apim-graphic, echo=FALSE, out.width="68%", fig.align="center"}
knitr::include_graphics(file.path("Embed", "APIM_ind_1.png"))
```

When both are requested, `prepare_dyad_data()` creates the APIM and DIM
predictors together. Omitting `role` states the exchangeability assumption,
while `seed` makes its arbitrary member contrast reproducible.

# Data structure

The 12 rows below only illustrate the prepared structure; the model
specifications assume study-sized data, and the output uses the cached fit to
the main tutorial's 11,000-row simulation.

```{r exchangeable-data-preparation}
exchangeable_raw <- tibble(
  coupleID = rep(31:32, each = 6),
  userID = rep(c("31_1", "31_2", "32_1", "32_2"), each = 3),
  diaryday = rep(0:2, 4),
  closeness = c(5.2, 5.5, 5.1, 4.8, 5.0, 5.3,
                4.9, 5.1, 5.0, 5.4, 5.2, 5.5),
  provided_support = c(4.1, 4.5, 4.3, 3.9, 4.2, 4.4,
                       4.0, 4.3, 4.1, 4.6, 4.5, 4.7)
)

df_long_dim <- prepare_dyad_data(
  data = exchangeable_raw,
  dyad = coupleID,
  member = userID,
  time = diaryday,
  predictors = provided_support,
  model_types = c("apim", "dim"),
  seed = 123
) |>
  mutate(
    diaryday_c = diaryday - mean(diaryday)
  )
```

```{r member-contrast-data-structure-table, echo=FALSE}
member_contrast_example <- df_long_dim |>
  select(
    userID, diaryday, .member_contrast_arbitrary,
    .provided_support_cwp_actor, .provided_support_cwp_partner
  )

kable(
  member_contrast_example,
  caption = "Prepared rows for an exchangeable dyad",
  booktabs = TRUE
) |>
  kable_styling(latex_options = "scale_down", font_size = 7)
```

`.member_contrast_arbitrary` is constant within person and opposite within
dyad. Its sign has no substantive meaning.

# glmmTMB model

```{r glmmTMB-model, eval=FALSE}
model_exch_long_apim_glmmtmb <- glmmTMB(
  closeness ~ 1 +

    # Pooled fixed effects: no gender-specific intercepts or slopes.
    diaryday_c +
    .provided_support_cwp_actor +
    .provided_support_cwp_partner +
    .provided_support_cbp_actor +
    .provided_support_cbp_partner +

    # Dyad-level means and slopes.
    (1 + diaryday_c +
       .provided_support_cwp_actor +
       .provided_support_cwp_partner | coupleID) +

    # Partner deviations from the dyad-level effects.
    # This block is kept separate from the dyad-level block so that the
    # common/couple-mean and deviation components are independent.
    (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) +

    # Same-day residual interdependence:
    # Independent common and member-deviation components imply homogeneous
    # residual variance and covariance.
    (1 | coupleID:diaryday) +
      (0 + .member_contrast_arbitrary | coupleID:diaryday),

  # The Gaussian residual covariance is represented by the two blocks above.
  # Without dispformula = ~ 0, glmmTMB would add a second independent
  # residual variance on top of the dyadic residual covariance structure.
  # For non-Gaussian families, dispformula controls family-specific dispersion
  # when one exists and may be ignored otherwise; do not use ~0 as this trick.
  dispformula = ~ 0,
  family = gaussian(),
  data = df_long_dim
)
```

The common and member-deviation components are rotated to member-level quantities in the
Back-transformation section.

# brms model

```{r brms-model, eval=FALSE}
formula_exch_long_apim <- bf(
  closeness ~ 1 +

    # Pooled fixed effects.
    diaryday_c +
    .provided_support_cwp_actor +
    .provided_support_cwp_partner +
    .provided_support_cbp_actor +
    .provided_support_cbp_partner +

    # Dyad-level means and slopes.
    (1 + diaryday_c +
       .provided_support_cwp_actor +
       .provided_support_cwp_partner | coupleID) +

    # Partner deviations from the dyad-level effects.
    (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) +

    # Native same-day residual correlation, allowed to be positive or negative.
    # With exactly two members per group, unstr() estimates the single pairwise
    # correlation. For triads or larger groups it estimates a fully
    # unstructured correlation matrix.
    unstr(time = .member_contrast_arbitrary, gr = coupleID:diaryday)

  # No sigma formula is specified here. Thus, sigma ~ 1 is implied,
  # meaning one homogeneous residual variance for both partners.
  # If sigma is modelled with its own formula, brms estimates the sigma
  # coefficients on the log scale, so those coefficients must be exponentiated
  # before reporting them as residual SDs.
)

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")
)
```

```{r load-cached-fit, include=FALSE}
cached_fit <- file.path("brms_cache", "model_apim_ind_long_apim.rds")
if (!file.exists(cached_fit)) {
  stop(
    "Cached brms fit not found: ", cached_fit,
    ". Render DyadicDataAnalysis.Rmd first, then render this file again."
  )
}
model_exch_long_apim <- readRDS(cached_fit)
```

The cached illustrative fit has no HMC warnings, but a few parameters, including the
intercept, mix slowly (maximum $\widehat{R} = 1.049$); refit longer before substantive
inference.

# Back-transformation

The model estimates a common/couple-mean block and a member-deviation block. For exchangeable dyads,
the full APIM random-effects matrix is recovered by rotating these two blocks:

$$
\operatorname{Var}(A) = \operatorname{Var}(B) =
\Sigma_{\text{common}} + \Sigma_{\text{deviation}}
$$

$$
\operatorname{Cov}(A, B) =
\Sigma_{\text{common}} - \Sigma_{\text{deviation}}
$$

This exact rotation assumes that the common/couple-mean and member-deviation random-effect blocks are
independent.

The helper function `summarize_exchangeable_apim_brms()` used below is defined in
`00_R_Functions/ReportModels.R` on GitHub. It reconstructs the common/couple-mean and member-deviation covariance
matrices for every posterior draw, rotates each draw into the full APIM
partner-level matrix, and then summarizes the transformed covariance, SD,
correlation, fixed-effect, and residual quantities.

For `glmmTMB`, `dyadMLM::recover_exchangeable_covariance()` provides the
algebraic point back-transformation. Interval estimates require a separate
parametric bootstrap or transformation of asymptotic parameter draws.

```{r backtransform-code}
apim_results <- summarize_exchangeable_apim_brms(
  model_exch_long_apim,
  gr = "coupleID",
  deviation_term = ".member_contrast_arbitrary",
  term_labels = c("int", "day", "actor_wp", "partner_wp"),
  partner_labels = c("A", "B")
)
```

# Variance-covariance matrix

```{r covariance-table}
kable(
  round(apim_results$full_covariance_matrix, 3),
  caption = "Back-transformed full APIM random-effects variance-covariance matrix",
  booktabs = TRUE
) |>
  kable_styling(latex_options = "scale_down", font_size = 7)
```

# SD table

```{r sd-table}
kable(
  apim_results$sd_summary,
  digits = 3,
  caption = "Back-transformed APIM random-effect standard deviations with 95% credible intervals"
)
```

# Correlation matrix

```{r correlation-table}
kable(
  round(apim_results$full_correlation_matrix, 3),
  caption = "Back-transformed APIM random-effect correlation matrix",
  booktabs = TRUE
) |>
  kable_styling(latex_options = "scale_down", font_size = 7)
```

# Final reporting summary

```{r final-summary-table, echo=FALSE}
apim_results$reporting_summary |>
  mutate(across(where(is.numeric), ~ round(.x, 3))) |>
  kable(
    caption = "Compact reporting summary",
    booktabs = TRUE
  ) |>
  kable_styling(latex_options = c("scale_down", "repeat_header"), font_size = 7)
```