DyadicDataAnalysis.Rmd

---
title: "Distinguishable and Exchangeable Dyads: Bayesian Multilevel Modelling"
subtitle: "Cross-Sectional and Intensive Longitudinal APIM and DIM"
author: "Pascal M. Küng"
date: "last-modified"
abstract: >
  This tutorial presents applied Bayesian multilevel modelling for dyadic data,
  covering both distinguishable and exchangeable dyads. We develop cross-sectional
  and intensive longitudinal variants of the Actor–Partner Interdependence Model
  (APIM) and the Dyad–Individual Model (DIM), demonstrate their algebraic
  equivalence, and discuss implications for interpretation. Using simulated data,
  we provide fully reproducible code to rebuild the models and assess parameter
  recovery. The tutorial offers comprehensive {brms} implementations, including
  random-effects structures, heterogeneous residual variances, compound-symmetry
  residuals, and optional AR(1) processes. We illustrate diagnostics and model
  comparison and show how to extract dyad-level variance–covariance components.
keywords:
  - dyadic data
  - APIM
  - DIM
  - Bayesian multilevel modelling
  - brms
  - cmdstanr
  - longitudinal diary data
  - compound symmetry
  - leave-one-out cross-validation
format:
  revealjs:
    embed-resources: true
    mathjax: embed
    theme: default
    slide-number: "c/t"
    progress: true
    hash: true
    preview-links: auto
    auto-animate: true
    code-line-numbers: true
    incremental: false
    code-overflow: wrap
    code-copy: true
    code-tools: false
    toc: false
    df-print: default
    transition: slide
    transition-speed: fast
    link-external-newwindow: true
    footer: "<span>Dyadic Multilevel Modelling • brms/cmdstanr • Pascal Küng</span>"
    css: Embed/styles.css
    include-after-body: Embed/fit-scroll-blocks.html
  pdf:
    pdf-engine: xelatex
    include-in-header:
      text: |
        % Code/output wrapping
        \usepackage{fvextra}
        \DefineVerbatimEnvironment{Highlighting}{Verbatim}{fontsize=\small,breaklines,breakanywhere,commandchars=\\\{\}}
        \fvset{breaklines,breakanywhere}

        % Tighter tables and more flexible line breaking
        \usepackage{etoolbox}
        \AtBeginEnvironment{longtable}{\setlength{\tabcolsep}{4pt}\renewcommand{\arraystretch}{0.95}\small}
        \AtBeginEnvironment{CSLReferences}{\small}
        \emergencystretch=2em
        \usepackage{pdflscape}
        \usepackage{needspace}
    include-before-body:
      text: |
        \RecustomVerbatimEnvironment{Highlighting}{Verbatim}{fontsize=\scriptsize,breaklines,breakanywhere,commandchars=\\\{\}}
    toc: false
    toc-depth: 2
    number-sections: true
    papersize: a4
    geometry: ["margin=1in"]
    colorlinks: true
    df-print: default
lang: en
bibliography:
  - Embed/Temp.bib
  - Embed/Rpackages.bib
csl: Embed/apa.csl
nocite: |
  @R-brms, @R-bmlm, @R-easystats, @R-DHARMa, @R-tidyverse, @R-cmdstanr, @R-loo, @R-rstan, @R-MASS, @R-nlme

---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE,
  cache.lazy = FALSE,
  cache.extra = paste(
    "dyadMLM-canonical-column-schema-v3",
    if (knitr::is_latex_output()) "pdf-output-font-v1" else "html-output"
  )
)
rm(list = ls())
gc()
default_output_hook <- knitr::knit_hooks$get("output")
knitr::knit_hooks$set(output = function(x, options) {
  out <- default_output_hook(x, options)
  if (knitr::is_latex_output()) {
    return(paste0("\\begingroup\\scriptsize\n", out, "\n\\endgroup\n"))
  }
  out
})
library(glue)
options(knitr.kable.NA = "", width = 120)
if (knitr::is_html_output()) {
  # kableExtra's HTML table helper initializes tooltips with jQuery.
  knitr::knit_meta_add(list(jquerylib::jquery_core()))
}
```

::: {.content-visible when-format="pdf"}

> **Interactive HTML version (recommended):**
> The tutorial is best viewed as [interactive slides](https://pascal-kueng.github.io/05DyadicDataAnalysis/DyadicDataAnalysis.html).
> This PDF is a static archival version for indexing, citation, and long-term reference.

**DOI:** [10.5281/zenodo.17400655](https://doi.org/10.5281/zenodo.17400655)

**Suggested citation:** Küng, P. M. (2026). *Distinguishable and Exchangeable Dyads: Bayesian Multilevel Modelling* (Version 2.0.12) [Report]. Zenodo. [https://doi.org/10.5281/zenodo.17400655](https://doi.org/10.5281/zenodo.17400655)

The [GitHub repository](https://github.com/Pascal-Kueng/05DyadicDataAnalysis) contains the source, helpers, model code, and quick guides.

:::

## 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:

* [Distinguishable dyads quick guide](distinguishable-dyads-quick-guide.pdf)
* [Exchangeable dyads quick guide](exchangeable-dyads-quick-guide.pdf)

## Preparing Dyadic Data with `dyadMLM`

[`dyadMLM`](https://pascal-kueng.github.io/dyadMLM/) prepares dyadic data for multilevel models [@dyadMLM2026]. 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](https://pascal-kueng.r-universe.dev/dyadMLM),
using CRAN for dependencies:

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

[Getting started](https://pascal-kueng.github.io/dyadMLM/articles/getting-started.html) · [Function reference](https://pascal-kueng.github.io/dyadMLM/reference/index.html) · [Source code](https://github.com/Pascal-Kueng/dyadMLM)

## Loading Libraries and setting CmdStan Backend

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

```{r load-libraries, results='hide', message=FALSE, warning=FALSE}
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'
)
```


```{r simulate-cross-sectional-data, echo=FALSE, results='hide', cache=TRUE}
# --- Packages (explicit namespacing avoids conflicts) ---
set.seed(123)

# --- 1. Define Simulation Parameters ---
# Sample Size
n_het <- 550 # Number of heterosexual couples
n_ss  <- 400 # Number of same-sex couples
n_total <- n_het + n_ss

# Residual SDs by gender / dyad type (unchanged)
sd_e_het_female <- 1.1
sd_e_het_male   <- 1.4
sd_e_ss_female  <- 0.9   # (SS individuals coded female in this example)
# sd_e_ss_male <- ...    # if you add male–male couples

# --- NEW: Residual correlations by dyad type (compound symmetry) ---
rho_e_het <- 0.30
rho_e_ss  <- 0.40

# Fixed Effects (the "true" values we want to recover)
gamma_0_het_female <- 5.5
gamma_0_het_male   <- 4.5
gamma_0_ss_person  <- 5.35

gamma_a_het_female <- 1.6
gamma_a_het_male   <- 1.8
gamma_a_ss_person  <- 0.2

gamma_p_het_female <- 0.3
gamma_p_het_male   <- 0.2
gamma_p_ss_person  <- 0.25

# Variance Components
sd_x <- 1.5      # SD of the predictor (communication)
r_x  <- 0.4      # Correlation between partners' communication

# --- UPDATED: smaller, more realistic random slope SDs ---
# (Feel free to tweak; these are modest to reflect heterogeneity without dominating)
sd_u0j <- 0.8    # baseline random intercept SD (unused directly below; see group-specific)
sd_uaj <- 0.15   # (legacy overall; using group-specific below)
sd_upj <- 0.12
cor_u0a <- 0.20
cor_u0p <- 0.15
cor_uap <- -0.10
sd_e    <- 1.2   # legacy scalar residual SD (not used now that we do sex-specific)

# --- 2. Simulate Dyad-Level Variables ---
dyad_data <- tibble::tibble(coupleID = 1:n_total)

# -- 2a: SEPARATE variance components for each group (UPDATED slopes) --

# Heterosexual couples
sd_u0j_het <- 0.8
sd_uaj_het <- 0.15   # was 0.3
sd_upj_het <- 0.12   # was 0.2
cor_u0a_het <- 0.20  # was 0.25
cor_u0p_het <- 0.15  # same
cor_uap_het <- -0.10 # same

# Same-sex couples
sd_u0j_ss <- 1.0
sd_uaj_ss <- 0.12   # was 0.2
sd_upj_ss <- 0.15   # was 0.25
cor_u0a_ss <- 0.10
cor_u0p_ss <- -0.05
cor_uap_ss <- 0.00

# -- 2b: Variance-covariance matrices --
cov_matrix_u_het <- matrix(c(
  sd_u0j_het^2,                      cor_u0a_het*sd_u0j_het*sd_uaj_het, cor_u0p_het*sd_u0j_het*sd_upj_het,
  cor_u0a_het*sd_u0j_het*sd_uaj_het, sd_uaj_het^2,                      cor_uap_het*sd_uaj_het*sd_upj_het,
  cor_u0p_het*sd_u0j_het*sd_upj_het, cor_uap_het*sd_uaj_het*sd_upj_het, sd_upj_het^2
), nrow = 3, byrow = TRUE)

cov_matrix_u_ss <- matrix(c(
  sd_u0j_ss^2,                     cor_u0a_ss*sd_u0j_ss*sd_uaj_ss, cor_u0p_ss*sd_u0j_ss*sd_upj_ss,
  cor_u0a_ss*sd_u0j_ss*sd_uaj_ss,  sd_uaj_ss^2,                   cor_uap_ss*sd_uaj_ss*sd_upj_ss,
  cor_u0p_ss*sd_u0j_ss*sd_upj_ss,  cor_uap_ss*sd_uaj_ss*sd_upj_ss, sd_upj_ss^2
), nrow = 3, byrow = TRUE)

# -- 2c: Generate random effects separately for each group --
random_effects_het <- MASS::mvrnorm(n = n_het, mu = c(0, 0, 0), Sigma = cov_matrix_u_het)
random_effects_ss  <- MASS::mvrnorm(n = n_ss,  mu = c(0, 0, 0), Sigma = cov_matrix_u_ss)

# -- 2d: Combine into one dataframe --
het_df_re <- tibble::tibble(
  coupleID = 1:n_het,
  u_0j = random_effects_het[, 1],
  u_aj = random_effects_het[, 2],
  u_pj = random_effects_het[, 3]
)

ss_df_re <- tibble::tibble(
  coupleID = (n_het + 1):n_total,
  u_0j = random_effects_ss[, 1],
  u_aj = random_effects_ss[, 2],
  u_pj = random_effects_ss[, 3]
)

combined_re <- dplyr::bind_rows(het_df_re, ss_df_re)

# Generate correlated predictor scores for partners (same as before)
cov_matrix_x <- matrix(c(sd_x^2, r_x*sd_x^2, r_x*sd_x^2, sd_x^2), nrow = 2)
partner_scores <- MASS::mvrnorm(n = n_total, mu = c(5, 5), Sigma = cov_matrix_x)

# Join combined RE with main dyad data
dyad_data <- dyad_data %>%
  dplyr::mutate(
    communication_p1 = partner_scores[, 1],
    communication_p2 = partner_scores[, 2]
  ) %>%
  dplyr::left_join(combined_re, by = "coupleID")

# --- 3. Create Person-Level (Long) Dataframe ---
person1 <- dyad_data %>% dplyr::select(coupleID, communication = communication_p1) %>% dplyr::mutate(p_num = 1)
person2 <- dyad_data %>% dplyr::select(coupleID, communication = communication_p2) %>% dplyr::mutate(p_num = 2)

df <- dplyr::bind_rows(person1, person2) %>%
  dplyr::arrange(coupleID, p_num) %>%
  dplyr::mutate(userID = paste0(coupleID, "_", p_num)) %>%
  dplyr::select(userID, coupleID, communication, p_num)

# Assign gender (1 = female, 2 = male)
het_couples <- 1:n_het
ss_couples  <- (n_het + 1):n_total

df <- df %>%
  dplyr::mutate(
    gender = dplyr::case_when(
      coupleID %in% het_couples & p_num == 1 ~ 1,
      coupleID %in% het_couples & p_num == 2 ~ 2,
      coupleID %in% ss_couples               ~ 1
    )
  )

temp_df_for_outcome <- df %>%
  dplyr::group_by(coupleID) %>%
  dplyr::mutate(
    communication_actor   = communication,
    communication_partner = communication[2:1]
  ) %>%
  dplyr::ungroup() %>%
  dplyr::left_join(dyad_data, by = "coupleID") %>%
  dplyr::mutate(
    # Center predictors before generating outcome (grand-mean)
    c_comm_actor   = communication_actor   - mean(communication_actor),
    c_comm_partner = communication_partner - mean(communication_partner),

    # Deterministic part
    satisfaction_deterministic = dplyr::case_when(
      gender == 1 & coupleID <= n_het ~ gamma_0_het_female + gamma_a_het_female * c_comm_actor + gamma_p_het_female * c_comm_partner,
      gender == 2 & coupleID <= n_het ~ gamma_0_het_male   + gamma_a_het_male   * c_comm_actor + gamma_p_het_male   * c_comm_partner,
      gender == 1 & coupleID >  n_het ~ gamma_0_ss_person  + gamma_a_ss_person  * c_comm_actor + gamma_p_ss_person  * c_comm_partner
    )
  )

# --- NEW: Residuals with compound symmetry (per dyad) & sex-specific SDs -------
resid_pairs <- tibble::tibble(coupleID = 1:n_total) %>%
  dplyr::mutate(
    sd1 = ifelse(coupleID <= n_het, sd_e_het_female, sd_e_ss_female),  # p1: het=female, ss=female
    sd2 = ifelse(coupleID <= n_het, sd_e_het_male,   sd_e_ss_female),  # p2: het=male,   ss=female
    rho = ifelse(coupleID <= n_het, rho_e_het,       rho_e_ss)
  ) %>%
  dplyr::rowwise() %>%
  dplyr::mutate(
    e_pair = list(MASS::mvrnorm(
      n = 1, mu = c(0, 0),
      Sigma = matrix(c(sd1^2, rho*sd1*sd2,
                       rho*sd1*sd2, sd2^2), nrow = 2)
    ))
  ) %>%
  dplyr::ungroup() %>%
  tidyr::unnest_wider(e_pair, names_sep = "") %>%
  dplyr::rename(e1 = e_pair1, e2 = e_pair2)

# Attach residuals and finish outcome
temp_df_for_outcome <- temp_df_for_outcome %>%
  dplyr::left_join(resid_pairs, by = "coupleID") %>%
  dplyr::mutate(
    e_i = ifelse(p_num == 1, e1, e2),
    satisfaction = satisfaction_deterministic +
      u_0j + (u_aj * c_comm_actor) + (u_pj * c_comm_partner) + e_i
  ) %>%
  dplyr::select(-e1, -e2)

# Finalize: ORIGINAL OUTPUT FORMAT
df <- temp_df_for_outcome %>%
  dplyr::select(userID, coupleID, gender, communication, satisfaction) %>%
  dplyr::mutate(
    gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))
  )

```


# Distinguishable Dyads

# Distinguishable Dyads - Cross-Sectional APIM

## Distinguishable Dyads - Data: Simulated Dyads

```{r data-dist-apim-raw-preview}
print_df(head(df))
```


## Distinguishable Dyads - Cross-Sectional APIM

![Distinguishable APIM](Embed/APIM_dist.png)
[e.g., @kenny2006dyadic; @kennyPartnerEffectsRelationship1999; @kennyDyadicDataAnalysis2011]


## Distinguishable Dyads - Preparing Data

```{r data-dist-apim-prepare, results='hide', message=FALSE, warning=FALSE, cache=TRUE, dependson='simulate-cross-sectional-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

```{r data-dist-apim-preview, echo=FALSE}

print_apim_preview(df_apim)

```

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

```{r fit-dist-apim-glmmtmb, results='hide', message=FALSE, warning=FALSE, cache=TRUE, dependson='data-dist-apim-prepare'}
#| code-line-numbers: "1-13|15-16|18-21|23-27|29-30"
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

```{r results-dist-apim-glmmtmb, echo=FALSE, class.output='model-summary-output', cache=TRUE, dependson='fit-dist-apim-glmmtmb'}

summary(model_dist_apim_glmmtmb)

```

## glmmTMB: Diagnostics with Near-Zero Gaussian Sigma {.smaller}

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 [@R-glmmTMB]. 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
[@R-DHARMa]:

```r
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.

```{r fit-dist-apim-brms, results='hide', message=FALSE, warning=FALSE, cache=TRUE, dependson='data-dist-apim-prepare'}
#| code-line-numbers: "1-13|15-21|23-28|30-41"
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

```{r results-dist-apim-brms-raw, echo=FALSE, class.output='model-summary-output', cache=TRUE, dependson='fit-dist-apim-brms'}
summary(model_dist_apim_b)
```

## Distinguishable Dyads - Poisson Covariance in brms {.smaller}

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 [@R-brms]. We simulate an actual count outcome
rather than discretize continuous `satisfaction`.

## Distinguishable Dyads - Fitting the Poisson Model in brms

```{r simulate-dist-apim-poisson-data, include=FALSE, cache=TRUE, dependson='data-dist-apim-prepare'}
set.seed(124)

latent_by_dyad <- df_apim |>
  dplyr::distinct(coupleID) |>
  dplyr::arrange(coupleID)

latent_effects <- MASS::mvrnorm(
  n = nrow(latent_by_dyad),
  mu = c(0, 0),
  Sigma = matrix(
    c(0.25^2, 0.30 * 0.25 * 0.40,
      0.30 * 0.25 * 0.40, 0.40^2),
    nrow = 2
  )
)

latent_by_dyad <- latent_by_dyad |>
  dplyr::mutate(
    latent_female = latent_effects[, 1],
    latent_male = latent_effects[, 2]
  )

df_apim_count <- df_apim |>
  dplyr::left_join(latent_by_dyad, by = "coupleID") |>
  dplyr::mutate(
    count_eta =
      .is_female * (
        log(6) +
          0.18 * (.communication_actor - mean(.communication_actor)) +
          0.06 * (.communication_partner - mean(.communication_partner)) +
          latent_female
      ) +
      .is_male * (
        log(5) +
          0.20 * (.communication_actor - mean(.communication_actor)) +
          0.04 * (.communication_partner - mean(.communication_partner)) +
          latent_male
      ),
    satisfaction_count = stats::rpois(dplyr::n(), exp(count_eta))
  ) |>
  dplyr::select(-latent_female, -latent_male, -count_eta)
```

```{r fit-dist-apim-brms-poisson, results='hide', message=FALSE, warning=FALSE, cache=TRUE, dependson='simulate-dist-apim-poisson-data'}
#| code-line-numbers: "1-13|15-21|23-28|29-41"
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

```{r results-dist-apim-brms-poisson-raw, echo=FALSE, class.output='model-summary-output', cache=TRUE, dependson='fit-dist-apim-brms-poisson'}
summary(model_dist_apim_c)
```

## Distinguishable Dyads - Fitting the Model in glmmTMB

```{r fit-dist-apim-glmmtmb-poisson, results='hide', message=FALSE, warning=FALSE, cache=TRUE, dependson='simulate-dist-apim-poisson-data'}
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

```{r results-dist-apim-glmmtmb-poisson-raw, echo=FALSE, class.output='model-summary-output', cache=TRUE, dependson='fit-dist-apim-glmmtmb-poisson'}
summary(model_glmmTMB_dist_apim_c)
```

## 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):

```{r diag-dist-apim-hmc, cache = TRUE, message=TRUE, dependson='fit-dist-apim-brms'}
# In brms 2.23, draws from either Stan backend are stored as a stanfit here.
rstan::check_hmc_diagnostics(model_dist_apim_b$fit)
```

## Distinguishable Dyads - Check Model Convergence and Fit


```{r eval=FALSE}
plot(model_dist_apim_b, ask = FALSE)
```

```{r echo=FALSE}
n_plot_variables <- if (knitr::is_latex_output()) 1 else 3
p1 <- plot(
  model_dist_apim_b,
  ask = FALSE,
  nvariables = n_plot_variables,
  plot = FALSE
)
p1[[2]] # Inspect all parameters; show a compact subset in this tutorial.
```
## Distinguishable Dyads - Check Model Convergence and Fit

```{r pp-dist-apim-dens, cache = TRUE, dependson='fit-dist-apim-brms'}
pp_check(model_dist_apim_b, 'dens_overlay_grouped', group = '.is_male')
```

## Distinguishable Dyads - Check Model Convergence and Fit

```{r pp-dist-apim-ecdf, cache = TRUE, dependson='fit-dist-apim-brms'}
pp_check(model_dist_apim_b, 'ecdf_overlay_grouped', group = '.is_male')
```

## Distinguishable Dyads - Check Model Convergence and Fit

```{r diag-dist-apim-dharma, cache = TRUE, dependson='fit-dist-apim-brms'}
# 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 {.smaller}

- 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
  [@R-brms].


## 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

```{r results-dist-apim-brms, echo = FALSE, cache=TRUE, dependson='fit-dist-apim-brms'}

summary_dist_apim_b <- summarize_brms(
    model_dist_apim_b,
    exponentiate_sigma = TRUE,
    stats_to_report = c('CI', 'Rhat', 'ESS')
    )

print_df(
  summary_dist_apim_b,
  rows_to_pack = list(
    'Fixed Effects' = c(1,6),
    'Residual Structure' = c(7,9)
    )
  )

```

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

```{r simulate-longitudinal-data, echo=FALSE, results='hide', cache=TRUE}
# Simulate Data
set.seed(123)

# ===============================
# 1) Design
# ===============================
n_couples <- 100
days      <- 0:54
T_per     <- length(days)

# People per couple: p1 = female, p2 = male
persons <- tibble(
  coupleID  = rep(1:n_couples, each = 2),
  p_num     = rep(1:2, times = n_couples),
  gender    = ifelse(p_num == 1, 1L, 2L),  # 1=female, 2=male
  is_female = as.integer(gender == 1L),
  is_male   = as.integer(gender == 2L),
  userID    = paste0(coupleID, "_", p_num)
)

# ===============================
# 2) Fixed effects (edit as desired)
# ===============================
# Intercepts & time trends by gender
b0_female <- 5.5; b0_male <- 4.7 # depends on centering, not meaningful when recovered.
b_time_female <-  0.01; b_time_male <- -0.005 # must be small because over 55 days!

# Between-person APIM (effects of person-mean support)
b_bp_actor_f <- 1.5; b_bp_partner_f <- 0.5
b_bp_actor_m <- 1.2; b_bp_partner_m <- 0.3

# Within-person APIM (effects of daily deviations)
b_wp_actor_f <- 0.4; b_wp_partner_f <- 0.2
b_wp_actor_m <- 0.2; b_wp_partner_m <- 0.1 # hard to detect!

# ===============================
# 3) Couple-level random effects (full var-cov)
#    Vector: [b0_m, b0_f, time_m, time_f, wpA_m, wpP_m, wpA_f, wpP_f]
# ===============================
sd_re <- c(0.70, 0.80, 0.015, 0.015, 0.10, 0.08, 0.10, 0.08)
R     <- matrix(0.15, 8, 8); diag(R) <- 1
Sigma_re <- diag(sd_re) %*% R %*% diag(sd_re)

RE <- MASS::mvrnorm(n_couples, mu = rep(0, 8), Sigma = Sigma_re) |>
  as_tibble(.name_repair = ~c("re_b0_m","re_b0_f","re_time_m","re_time_f",
                              "re_wpA_m","re_wpP_m","re_wpA_f","re_wpP_f")) |>
  mutate(coupleID = 1:n_couples)

# ===============================
# 4) Provided support (predictor)
#    - Person means correlated within couple
#    - Daily AR(1) within-person; cross-partner correlated innovations
# ===============================
mu_sup_f <- 5.0; sd_mu_f <- 0.8
mu_sup_m <- 4.8; sd_mu_m <- 0.8
rho_couple_mean <- 0.50
Sigma_mu <- matrix(c(sd_mu_f^2, rho_couple_mean*sd_mu_f*sd_mu_m,
                     rho_couple_mean*sd_mu_f*sd_mu_m, sd_mu_m^2), 2)

mu_pairs <- MASS::mvrnorm(n_couples, mu = c(mu_sup_f, mu_sup_m), Sigma = Sigma_mu)

mu_df <- tibble(
  coupleID = rep(1:n_couples, each = 2),
  p_num    = rep(1:2, times = n_couples),
  mu_sup   = c(rbind(mu_pairs[,1], mu_pairs[,2]))
)

phi_sup <- 0.4
sd_sup_wp  <- 0.8
sd_eta_sup <- sqrt(sd_sup_wp^2 * (1 - phi_sup^2))
rho_sup_day <- 0.3
Sigma_eta_day <- matrix(c(sd_eta_sup^2, rho_sup_day*sd_eta_sup^2,
                          rho_sup_day*sd_eta_sup^2, sd_eta_sup^2), 2)

sim_support_one_couple <- function(cid) {
  eta <- MASS::mvrnorm(T_per, mu = c(0,0), Sigma = Sigma_eta_day)
  sup <- matrix(NA_real_, nrow = T_per, ncol = 2)
  sup[1,] <- eta[1,] / sqrt(1 - phi_sup^2)  # stationary init
  for (t in 2:T_per) sup[t,] <- phi_sup * sup[t-1,] + eta[t,]
  tibble(coupleID = cid, diaryday = days, sup_f = sup[,1], sup_m = sup[,2])
}

support_time <- purrr::map_dfr(1:n_couples, sim_support_one_couple)

support_long <- support_time |>
  tidyr::pivot_longer(c(sup_f, sup_m), names_to = "who", values_to = "sup_cwp") |>
  dplyr::mutate(p_num = ifelse(who == "sup_f", 1L, 2L)) |>
  dplyr::select(-who) |>
  dplyr::left_join(mu_df, by = c("coupleID","p_num")) |>
  dplyr::mutate(provided_support = mu_sup + sup_cwp) |>
  dplyr::select(coupleID, diaryday, p_num, provided_support, sup_cwp)

# ===============================
# 5) Assemble panel & partner lookups (safe keys)
# ===============================
panel <- tidyr::expand_grid(coupleID = 1:n_couples, diaryday = days, p_num = 1:2) |>
  dplyr::left_join(persons, by = c("coupleID","p_num")) |>
  dplyr::left_join(support_long, by = c("coupleID","diaryday","p_num"))

panel <- panel |>
  dplyr::group_by(coupleID, diaryday) |>
  dplyr::mutate(
    provided_support_partner = provided_support[3 - p_num],
    sup_cwp_partner          = sup_cwp[3 - p_num]
  ) |>
  dplyr::ungroup()

mu_actor_tbl   <- mu_df |> dplyr::rename(mu_actor = mu_sup)
mu_partner_tbl <- mu_df |> dplyr::mutate(p_num = 3 - p_num) |> dplyr::rename(mu_partner = mu_sup)

panel <- panel |>
  dplyr::left_join(mu_actor_tbl,   by = c("coupleID","p_num")) |>
  dplyr::left_join(mu_partner_tbl, by = c("coupleID","p_num"))

# ===============================
# 6) Same-day couple×day random effect (common shocks)
# ===============================
use_day_RE <- TRUE     # set FALSE to remove same-day shock
sd_day <- 0.17
if (use_day_RE) {
  day_shock <- tidyr::expand_grid(coupleID = 1:n_couples, diaryday = days) |>
    dplyr::mutate(z_day = rnorm(n(), 0, sd_day))
  panel <- panel |> dplyr::left_join(day_shock, by = c("coupleID","diaryday"))
} else {
  panel <- panel |> dplyr::mutate(z_day = 0)
}

# ===============================
# 7) Deterministic mean of daily closeness
# ===============================
grand_mu <- mean(mu_df$mu_sup)  # single grand mean for CBP

panel <- panel |>
  dplyr::left_join(RE, by = "coupleID") |>
  dplyr::mutate(
    # fixed parts by gender
    fix_b0   = ifelse(is_male==1, b0_male, b0_female),
    fix_time = ifelse(is_male==1, b_time_male, b_time_female),
    fix_bp_actor   = ifelse(is_male==1, b_bp_actor_m,  b_bp_actor_f),
    fix_bp_partner = ifelse(is_male==1, b_bp_partner_m,b_bp_partner_f),
    fix_wp_actor   = ifelse(is_male==1, b_wp_actor_m,  b_wp_actor_f),
    fix_wp_partner = ifelse(is_male==1, b_wp_partner_m,b_wp_partner_f),

    # couple-level REs by gender
    re_b0      = ifelse(is_male==1, re_b0_m,   re_b0_f),
    re_time    = ifelse(is_male==1, re_time_m, re_time_f),
    re_wp_actor= ifelse(is_male==1, re_wpA_m,  re_wpA_f),
    re_wp_part = ifelse(is_male==1, re_wpP_m,  re_wpP_f),

    # latent components
    cwp_actor   = sup_cwp,
    cwp_partner = sup_cwp_partner,
    cbp_actor   = mu_actor   - grand_mu,
    cbp_partner = mu_partner - grand_mu,

    # deterministic mean
    mu_it =
      fix_b0 + fix_time * diaryday +
      fix_bp_actor   * cbp_actor   + fix_bp_partner * cbp_partner +
      fix_wp_actor   * cwp_actor   + fix_wp_partner * cwp_partner +
      re_b0 + re_time * diaryday +
      re_wp_actor * cwp_actor + re_wp_part * cwp_partner +
      z_day
  )

# ===============================
# 8) Residuals: AR(1) within person, gender-specific σ
# ===============================
phi_eps <- 0.10
sd_eps_f <- 0.70; sd_eps_m <- 1.10
sd_eta_f <- sqrt(sd_eps_f^2 * (1 - phi_eps^2))
sd_eta_m <- sqrt(sd_eps_m^2 * (1 - phi_eps^2))

simulate_ar1_resid <- function(n, sd_eta, phi=phi_eps) {
  e <- numeric(n)
  e[1] <- stats::rnorm(1, 0, sd = sd_eta / sqrt(1 - phi^2))
  if (n >= 2) for (t in 2:n) e[t] <- phi * e[t-1] + stats::rnorm(1, 0, sd_eta)
  e
}

# Vectorized generation per user
idx <- panel |>
  dplyr::arrange(userID, diaryday) |>
  dplyr::distinct(userID, is_male)

eps_list <- mapply(
  FUN = simulate_ar1_resid,
  n = T_per,
  sd_eta = ifelse(idx$is_male==1, sd_eta_m, sd_eta_f),
  SIMPLIFY = FALSE
)
eps_vec <- unlist(eps_list, use.names = FALSE)

panel <- panel |>
  dplyr::arrange(userID, diaryday) |>
  dplyr::mutate(e_it = eps_vec) |>
  dplyr::arrange(coupleID, diaryday, p_num)

# ===============================
# 9) Outcomes & exports
# ===============================
df_long <- panel |>
  dplyr::mutate(closeness = mu_it + e_it) |>
  dplyr::select(userID, coupleID, diaryday, gender,
                closeness, provided_support) |>
  dplyr::arrange(coupleID, userID, diaryday) |>
  dplyr::mutate(
    gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))
  )

```

```{r data-dist-long-raw-preview}
print_df(head(df_long))
```

## Distinguishable Dyads - L-APIM - Data Preparation

```{r data-dist-long-apim-reshape, results='hide', warning=FALSE, message=FALSE, cache=TRUE, dependson='simulate-longitudinal-data'}

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


```


```{r data-long-preview-helper, echo = FALSE, results='hide', warning=FALSE, message=FALSE}

# Function for preview (not needed for analysis!)
print_couple_preview <- function(df, couple_id,
                                 days = 0:3) {
  # Filter/sort/select with original types
  dfc0 <- df %>%
    dplyr::filter(coupleID == couple_id) %>%
    dplyr::arrange(userID, diaryday)

  # Identify partners
  partners <- dfc0 %>%
    dplyr::distinct(userID) %>%
    dplyr::arrange(userID) %>%
    dplyr::pull(userID)

  if (length(partners) < 2)
    stop("This couple doesn't have two partners in the data.")

  # Slice by partner & day BEFORE coercion
  p1 <- dfc0 %>% dplyr::filter(userID == partners[1], diaryday %in% days)
  p2 <- dfc0 %>% dplyr::filter(userID == partners[2], diaryday %in% days)

  # Coerce everything to character explicitly (belt-and-suspenders on coupleID)
  to_chr <- function(x) x %>%
    dplyr::mutate(
      dplyr::across(where(is.numeric), ~round(., 2)),
      dplyr::across(dplyr::everything(), as.character),
      coupleID = as.character(coupleID)
      )

  p1c <- to_chr(p1)
  p2c <- to_chr(p2)

  # 1-row "..." tibble with same columns (all character)
  colset <- names(dfc0)
  ellipsis_row <- tibble::as_tibble(setNames(rep(list("..."), length(colset)), colset))

  dplyr::bind_rows(p1c, ellipsis_row, p2c)
}

print_couple_preview_for_output <- function(df, couple_id, days = 0:3,
                                            pdf_cols = NULL) {
  out <- print_couple_preview(df, couple_id = couple_id, days = days)
  if (knitr::is_latex_output() && !is.null(pdf_cols)) {
    out <- out |>
      dplyr::select(dplyr::any_of(pdf_cols))
  }
  out
}


```

## Distinguishable Dyads - L-APIM - Prepared Data

```{r data-dist-long-apim-preview, echo=FALSE}

print_df(
  print_couple_preview_for_output(
    df_long_apim,
    couple_id = "31",
    pdf_cols = c(
      "userID", "diaryday", "diaryday_c", ".is_male", ".is_female",
      "closeness", ".provided_support_cwp_actor",
      ".provided_support_cwp_partner"
    )
  ),
  font_size = if (knitr::is_latex_output()) 6 else 20
)

```


## Distinguishable Dyads - L-APIM - Centering

With repeated measures we need to center variables in order to not conflate levels of analysis [@bolgerIntensiveLongitudinalMethods2013]:

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

```{r fit-dist-apim-long-nlme, cache = TRUE, eval = FALSE}
#| classes: code-45vh

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

```{r fit-dist-apim-long-glmmtmb, results='hide', message=FALSE, warning=FALSE, cache=TRUE, dependson='data-dist-long-apim-reshape'}
#| code-line-numbers: "1-6|8-13|15-20|22-26|28-34"
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

```{r results-dist-apim-long-glmmtmb, echo=FALSE, class.output='model-summary-output', cache=TRUE, dependson='fit-dist-apim-long-glmmtmb'}
summary(model_dist_apim_long_glmmtmb)
```

## Comparing outputs between nlme and glmmTMB {.model-comparison}

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

::: {.model-compare-note}
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.
:::

::: {.model-compare-grid}
::: {.model-compare-panel}
**nlme:**
<pre class="model-compare-output"><code>                     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</code></pre>
:::
::: {.model-compare-panel}
**glmmTMB:**
<pre class="model-compare-output"><code>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</code></pre>
:::
:::

## 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!

```{r fit-dist-apim-long-simple-brms, results='hide', message=FALSE, warning=FALSE, cache=TRUE, dependson='data-dist-long-apim-reshape'}
#| code-line-numbers: "1-6|8-13|15-20|22-26|28-34|36-48|50-61"
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)

```{r results-dist-apim-long-simple-brms, echo=FALSE, cache=TRUE, dependson='fit-dist-apim-long-simple-brms'}

summary_dist_apim_long_simple <- summarize_brms(
    model_dist_apim_long_simple,
    exponentiate_sigma = TRUE,
    stats_to_report = c('CI', 'Rhat', 'ESS')
    )

print_df(
  summary_dist_apim_long_simple,
  rows_to_pack = list(
    'Fixed Effects' = c(1,12),
    'Random Effects' = c(13, 23),
    'Residual Structure' = c(24, 25)
    ),
  scroll_height = '55vh'
  )

```

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

## Distinguishable Dyads - L-APIM - AR(1)? {.smaller}

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 [@R-brms].

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

```r
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 [@R-glmmTMB].

## 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 [@nickellBiasesDynamicModels1981], 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`.


```{r fit-dist-apim-long-complex-brms, results='hide', message=FALSE, warning=FALSE, cache=TRUE, dependson='data-dist-long-apim-reshape'}
#| classes: code-50vh
#| code-line-numbers: "27-35|37-43"

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

```{r results-dist-apim-long-complex-brms, echo=FALSE, cache=TRUE, dependson='fit-dist-apim-long-complex-brms'}

summary_dist_apim_long_complex <- summarize_brms(
    model_dist_apim_long_complex,
    exponentiate_sigma = TRUE,
    stats_to_report = c('CI', 'Rhat', 'ESS')
    )

print_df(
  summary_dist_apim_long_complex,
  rows_to_pack = list(
    'Fixed Effects' = c(1,12),
    'Random Effects' = c(13, 49),
    'Residual Structure' = c(50, 51)
    ),
  scroll_height = '55vh'
  )

```


## Maximal Model: Diagnostics and Prediction Target {.smaller}

- 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
  [@buerknerEfficientLeaveOneOut2021].
- New-dyad prediction would instead require a grouped validation design
  [@loo2017a].

## Maximal Model: PSIS-LOO Diagnostics {.smaller}

```{r diag-dist-long-model-loo, eval=FALSE}
loo1 <- loo(model_dist_apim_long_simple)
loo2 <- loo(model_dist_apim_long_complex)

loo::pareto_k_table(loo1)
loo::pareto_k_table(loo2)
```

```{=html}
<pre class="model-compare-output"><code>All Pareto k estimates are good (k &lt; 0.7).

All Pareto k estimates are good (k &lt; 0.7).</code></pre>
```

Pareto $k$ diagnoses the importance-sampling approximation; high values are
not themselves evidence of overfitting [@loo2017a].

## Maximal Model: LOO Comparison {.smaller}

```{r diag-dist-long-model-loo-compare, eval=FALSE}
a <- loo_compare(loo1, loo2)
print(a)
```

```{=html}
<pre class="model-compare-output"><code>                             elpd_diff se_diff
model_dist_apim_long_complex   0.0       0.0
model_dist_apim_long_simple  -15.4       6.7</code></pre>
```

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
[@loo2025a]. If performance were similar, theory, identifiability, and
parsimony should guide simplification
[@barrRandomEffectsStructure2013; @delrosarioPracticalGuideSpecifying2025].

```{=latex}
\begin{landscape}
```

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

```{r results-dist-long-ground-truth, echo = FALSE, message=FALSE, warning=FALSE, cache=TRUE, dependson=c('fit-dist-apim-long-complex-brms', 'simulate-longitudinal-data', 'data-dist-long-apim-reshape')}

# HOW WELL DO THE RECOVERED PARAMETERS FIT THE SIMULATED VALUES?

library(posterior)

fit <- model_dist_apim_long_complex  # your model

# ---------------------------
# 1) Empirical "truth" from this simulation draw
# ---------------------------
mean_day <- mean(days)  # e.g., 27 for 0:54

# Fixed-effect truth on the same centered-predictor scale used by the model.
fixed_truth_data <- panel |>
  dplyr::mutate(
    fixed_only =
      fix_b0 + fix_time * diaryday +
      fix_bp_actor   * cbp_actor   + fix_bp_partner * cbp_partner +
      fix_wp_actor   * cwp_actor   + fix_wp_partner * cwp_partner
  ) |>
  dplyr::select(userID, coupleID, diaryday, fixed_only) |>
  dplyr::inner_join(df_long_apim, by = c("userID", "coupleID", "diaryday"))

fixed_truth_fit <- stats::lm(
  fixed_only ~
    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,
  data = fixed_truth_data
)
fixed_truth <- coef(fixed_truth_fit)

# Random effects on the same centered-time scale used by the fitted model.
true_re_tbl <- tibble::tibble(
  .is_male = RE$re_b0_m + RE$re_time_m * mean_day,
  .is_female = RE$re_b0_f + RE$re_time_f * mean_day,
  `.is_male:diaryday_c` = RE$re_time_m,
  `.is_female:diaryday_c` = RE$re_time_f,
  `.is_male:.provided_support_cwp_actor` = RE$re_wpA_m,
  `.is_male:.provided_support_cwp_partner` = RE$re_wpP_m,
  `.is_female:.provided_support_cwp_actor` = RE$re_wpA_f,
  `.is_female:.provided_support_cwp_partner` = RE$re_wpP_f
)

re_sd_truth <- purrr::map_dbl(true_re_tbl, sd)
re_cor_truth <- stats::cor(true_re_tbl)

re_sd_labels <- c(
  ".is_male" = "Couple intercept SD (male)",
  ".is_female" = "Couple intercept SD (female)",
  ".is_male:diaryday_c" = "Time slope SD (male)",
  ".is_female:diaryday_c" = "Time slope SD (female)",
  ".is_male:.provided_support_cwp_actor" = "Actor WP slope SD (male)",
  ".is_male:.provided_support_cwp_partner" = "Partner WP slope SD (male)",
  ".is_female:.provided_support_cwp_actor" = "Actor WP slope SD (female)",
  ".is_female:.provided_support_cwp_partner" = "Partner WP slope SD (female)"
)

re_cor_labels <- c(
  ".is_male" = "Intercept (male)",
  ".is_female" = "Intercept (female)",
  ".is_male:diaryday_c" = "Time slope (male)",
  ".is_female:diaryday_c" = "Time slope (female)",
  ".is_male:.provided_support_cwp_actor" = "Actor WP slope (male)",
  ".is_male:.provided_support_cwp_partner" = "Partner WP slope (male)",
  ".is_female:.provided_support_cwp_actor" = "Actor WP slope (female)",
  ".is_female:.provided_support_cwp_partner" = "Partner WP slope (female)"
)

true_resid_tbl <- panel |>
  dplyr::mutate(resid_total = z_day + e_it)

true_sigma_m <- sd(true_resid_tbl$resid_total[true_resid_tbl$is_male == 1])
true_sigma_f <- sd(true_resid_tbl$resid_total[true_resid_tbl$is_female == 1])

true_resid_cor <- true_resid_tbl |>
  dplyr::select(coupleID, diaryday, gender, resid_total) |>
  tidyr::pivot_wider(names_from = gender, values_from = resid_total, names_prefix = "gender_") |>
  dplyr::summarise(r = stats::cor(gender_1, gender_2)) |>
  dplyr::pull(r)

fixed_truth_tbl <- tibble::tribble(
  ~group, ~parameter, ~coef_name,
  "Fixed", "Intercept (male)", ".is_male",
  "Fixed", "Intercept (female)", ".is_female",
  "Fixed", "Time slope (male)", ".is_male:diaryday_c",
  "Fixed", "Time slope (female)", ".is_female:diaryday_c",
  "Between-person", "Actor BP (male)", ".is_male:.provided_support_cbp_actor",
  "Between-person", "Partner BP (male)", ".is_male:.provided_support_cbp_partner",
  "Between-person", "Actor BP (female)", ".is_female:.provided_support_cbp_actor",
  "Between-person", "Partner BP (female)", ".is_female:.provided_support_cbp_partner",
  "Within-person", "Actor WP (male)", ".is_male:.provided_support_cwp_actor",
  "Within-person", "Partner WP (male)", ".is_male:.provided_support_cwp_partner",
  "Within-person", "Actor WP (female)", ".is_female:.provided_support_cwp_actor",
  "Within-person", "Partner WP (female)", ".is_female:.provided_support_cwp_partner"
) |>
  dplyr::mutate(true = unname(fixed_truth[coef_name]), scale = "identity") |>
  dplyr::select(group, parameter, true, scale)

re_sd_truth_tbl <- tibble::tibble(
  group = "Random effects (SD)",
  parameter = unname(re_sd_labels[names(re_sd_truth)]),
  true = unname(re_sd_truth),
  scale = "identity"
)

re_cor_truth_tbl <- t(combn(names(re_cor_labels), 2)) |>
  as.data.frame() |>
  dplyr::rename(param1 = V1, param2 = V2) |>
  dplyr::mutate(
    group = "Random effects (cor)",
    parameter = paste0("cor(", re_cor_labels[param1], ", ", re_cor_labels[param2], ")"),
    true = purrr::map2_dbl(param1, param2, ~ re_cor_truth[.x, .y]),
    scale = "identity"
  ) |>
  dplyr::select(group, parameter, true, scale)

resid_truth_tbl <- tibble::tribble(
  ~group, ~parameter, ~true, ~scale,
  "Residuals", "Same-day residual correlation", true_resid_cor, "identity",
  "Residuals", "Residual SD (male) (exp)", true_sigma_m, "exp",
  "Residuals", "Residual SD (female) (exp)", true_sigma_f, "exp"
)

truth_tbl <- dplyr::bind_rows(
  fixed_truth_tbl,
  re_sd_truth_tbl,
  re_cor_truth_tbl,
  resid_truth_tbl
)


# ---------------------------
# 2) Extract estimates from the model (means + 95% CI)
# ---------------------------

# Helper to get posterior summaries for specific parameter names (regex)
get_summaries <- function(fit, pattern) {
  sm <- brms::posterior_summary(fit, pars = pattern, robust = FALSE, probs = c(.025,.975))
  if (is.null(dim(sm))) return(tibble())
  as_tibble(sm, rownames = "param") |>
    rename(Estimate = Estimate, Lower = Q2.5, Upper = Q97.5)
}

# (A) Population-level (fixed) effects — use posterior_summary so names have b_
fix_all <- get_summaries(fit, pattern = "^b_")

# Map model param names -> display labels (same as before)
map_fixed <- tibble::tribble(
  ~param_regex,                              ~label,
  "^b_\\.is_male$",                                      "Intercept (male)",
  "^b_\\.is_female$",                                    "Intercept (female)",
  "^b_\\.is_male:diaryday_c$",                           "Time slope (male)",
  "^b_\\.is_female:diaryday_c$",                         "Time slope (female)",
  "^b_\\.is_male:\\.provided_support_cbp_actor$",       "Actor BP (male)",
  "^b_\\.is_male:\\.provided_support_cbp_partner$",     "Partner BP (male)",
  "^b_\\.is_female:\\.provided_support_cbp_actor$",     "Actor BP (female)",
  "^b_\\.is_female:\\.provided_support_cbp_partner$",   "Partner BP (female)",
  "^b_\\.is_male:\\.provided_support_cwp_actor$",       "Actor WP (male)",
  "^b_\\.is_male:\\.provided_support_cwp_partner$",     "Partner WP (male)",
  "^b_\\.is_female:\\.provided_support_cwp_actor$",     "Actor WP (female)",
  "^b_\\.is_female:\\.provided_support_cwp_partner$",   "Partner WP (female)"
)

fix_labeled <- map_fixed |>
  rowwise() |>
  mutate(
    out = list(
      fix_all |>
        dplyr::filter(stringr::str_detect(param, param_regex)) |>
        dplyr::mutate(parameter = label, scale = "identity") |>
        dplyr::select(parameter, Estimate, Lower, Upper, scale)
    )
  ) |>
  ungroup() |>
  dplyr::select(out) |>
  tidyr::unnest(out)


# (B) Group-level SDs and correlations (random effects)
sd_all <- get_summaries(fit, pattern = "^sd_")

sd_labeled <- sd_all |>
  dplyr::mutate(
    component = stringr::str_remove(param, "^sd_coupleID__"),
    parameter = unname(re_sd_labels[component]),
    scale = "identity"
  ) |>
  dplyr::filter(!is.na(parameter)) |>
  dplyr::select(parameter, Estimate, Lower, Upper, scale)

cor_all <- get_summaries(fit, pattern = "^cor_coupleID__")

cor_labeled <- cor_all |>
  dplyr::mutate(
    components = stringr::str_remove(param, "^cor_coupleID__"),
    param1 = stringr::str_split_fixed(components, "__", 2)[, 1],
    param2 = stringr::str_split_fixed(components, "__", 2)[, 2],
    parameter = paste0("cor(", re_cor_labels[param1], ", ", re_cor_labels[param2], ")"),
    scale = "identity"
  ) |>
  dplyr::filter(!is.na(re_cor_labels[param1]), !is.na(re_cor_labels[param2])) |>
  dplyr::select(parameter, Estimate, Lower, Upper, scale)

# (C) Residual sigmas (distributional parameters on log scale)
# With sigma ~ 0 + .is_male + .is_female, brms names population-level params:
#   b_sigma_.is_male, b_sigma_.is_female  (on log scale)
draws <- as_draws_df(fit)
resid_map <- tribble(
  ~param_draw,            ~label,
  "b_sigma_.is_male",      "Residual SD (male) (exp)",
  "b_sigma_.is_female",    "Residual SD (female) (exp)"
)

resid_summ <- resid_map |>
  rowwise() |>
  mutate(
    vals = list(exp(draws[[param_draw]])),  # exp to go from log-sigma to sigma
    Estimate = mean(vals),
    Lower = quantile(vals, 0.025),
    Upper = quantile(vals, 0.975),
    scale = "exp"
  ) |>
  ungroup() |>
  dplyr::select(parameter = label, Estimate, Lower, Upper, scale)

cortime_all <- get_summaries(fit, pattern = "^cortime")
cortime_summ <- if (nrow(cortime_all) > 0) {
  cortime_all |>
    dplyr::slice(1) |>
    dplyr::mutate(
      parameter = "Same-day residual correlation",
      scale = "identity"
    ) |>
    dplyr::select(parameter, Estimate, Lower, Upper, scale)
} else {
  tibble::tibble(parameter = character(), Estimate = numeric(), Lower = numeric(), Upper = numeric(), scale = character())
}

# Combine estimates
est_tbl <- bind_rows(fix_labeled, sd_labeled, cor_labeled, cortime_summ, resid_summ)

# ---------------------------
# 3) Join truth + estimates and compute diagnostics
# ---------------------------
full_tbl <- truth_tbl |>
  dplyr::select(group, parameter, true, scale) |>
  left_join(est_tbl, by = c("parameter","scale")) |>
  mutate(
    covered  = !is.na(Lower) & true >= pmin(Lower, Upper) & true <= pmax(Lower, Upper),
    abs_error = abs(Estimate - true),
    Recovery = dplyr::case_when(
      is.na(Estimate)                ~ "—",
      covered                        ~ "Covered",
      TRUE                           ~ "Miss"
    ),
    `95% CI` = ifelse(is.na(Lower), "", sprintf("[%.3f, %.3f]", Lower, Upper))
  ) |>
  dplyr::select(
    Group = group,
    Parameter = parameter,
    `True (Simulated)` = true,
    Estimate, `95% CI`,
    Coverage = covered,
    `Abs. Error` = abs_error,
    Recovery
  ) |>
  mutate(
    Coverage = ifelse(is.na(Estimate), "", ifelse(Coverage, "Yes", "No")),
    `Abs. Error` = ifelse(is.na(Estimate), "", sprintf("%.3f", `Abs. Error`))
  )

# ---------------------------
# 4) Render (scrollable)
# ---------------------------
print_df(
  full_tbl,
  scroll_height = '50vh'
)



```

```{=latex}
\end{landscape}
```










# Exchangeable Dyads

# Exchangeable Dyads - Cross-Sectional

## Exchangeable Dyads - Cross-Sectional APIM

![Exchangeable APIM](Embed/APIM_ind_1.png)
[e.g., @delrosarioPracticalGuideSpecifying2025; @kenny_ackerman_2023]

## Exchangeable Dyads - Cross-Sectional APIM

![Exchangeable APIM](Embed/APIM_ind_3.png){width=72%}
[e.g., @delrosarioPracticalGuideSpecifying2025; @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., @delrosarioPracticalGuideSpecifying2025; @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

```{r data-exch-prepare, results='hide', message=FALSE, warning=FALSE, cache=TRUE, dependson='simulate-cross-sectional-data'}
# Reuse the distinguishable preparation, including its arbitrary
# -1/+1 member contrast.
df_exchangeable_apim <- df_apim
```

## Exchangeable Dyads - Cross-Sectional APIM: Prepared Data

```{r data-exch-apim-preview}

print_apim_preview(df_exchangeable_apim)

```

## 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:

```{r fit-exch-apim-glmmtmb, results='hide', message=FALSE, warning=FALSE, cache=TRUE, dependson='data-exch-prepare'}
#| code-line-numbers: "1-5|7-8|10-16"
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 {.smaller}

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

```{r results-exch-apim-glmmtmb, echo=FALSE, class.output='model-summary-output', cache=TRUE, dependson='fit-exch-apim-glmmtmb'}
dyadMLM::recover_exchangeable_covariance(
  model_exch_apim_glmmtmb,
  block_pairings = list(
    shared_block = "(1 | coupleID)",
    difference_block = "(0 + .member_contrast_arbitrary | coupleID)",
    difference_indicator = ".member_contrast_arbitrary"
  )
)
```

## Exchangeable Dyads - Cross-Sectional APIM: Interpreting the Covariance {.smaller}

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$ [@delrosarioPracticalGuideSpecifying2025]. 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:

```{r fit-exch-apim-brms, results='hide', message=FALSE, warning=FALSE, cache=TRUE, dependson='data-exch-prepare'}
#| code-line-numbers: "1-13|15-20|22-33"
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

```{r results-exch-apim-brms, echo=FALSE, cache=TRUE, dependson='fit-exch-apim-brms'}

summary_exch_apim <- summarize_brms(
    model_exch_apim,
    stats_to_report = c('CI', 'Rhat', 'ESS'),
    exponentiate_sigma = TRUE
    )

communication_grand_mean <- mean(df_exchangeable_apim$communication)
apim_centered_intercept <-
  fixef(model_exch_apim)["Intercept", "Estimate"] +
  communication_grand_mean * (
    fixef(model_exch_apim)[".communication_actor", "Estimate"] +
    fixef(model_exch_apim)[".communication_partner", "Estimate"]
  )

print_df(
  summary_exch_apim,
  rows_to_pack = list(
    'Fixed Effects' = c(1,3),
    'Residual Structure' = c(4,5)
    )
  )

```


## Constraints vs Prediction {.smaller}

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

```{r compare-exch-dist-apim-glmmtmb, echo=FALSE, cache=TRUE, dependson=c('fit-dist-apim-glmmtmb', 'fit-exch-apim-glmmtmb')}
comparison_exch_dist <- dyadMLM::compare_nested_models(
  model_dist_apim_glmmtmb,
  model_exch_apim_glmmtmb
)

if (knitr::is_latex_output()) {
  comparison_table <- comparison_exch_dist
  class(comparison_table) <- c("anova", "data.frame")
  print(comparison_table)
  cat(
    "\n",
    paste(strwrap(attr(comparison_exch_dist, "conclusion"), width = 82),
          collapse = "\n"),
    "\n",
    sep = ""
  )
} else {
  print(comparison_exch_dist)
}
```

The following observation-level PSIS-LOO comparison instead evaluates
conditional predictive performance [@buerknerEfficientLeaveOneOut2021]. It is
not a formal test of distinguishability or prediction of a new dyad.

## Exchangeable Dyads - APIM: LOO Comparison {.smaller}

| 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 [@loo2017a; @loo2025a].

## Exchangeable Dyads - Cross-Sectional DIM

```{r diagram-exch-dim, echo=FALSE, cache=TRUE}

if (knitr::is_html_output()) {
  DiagrammeR::grViz(glue("
    digraph DIM {{
      graph [
        size='11,11!'
        layout=dot
        rankdir='LR'
        label='Dyad–Individual Model (DIM)'
        labelloc=t
      ]

      node [
        shape=box
        fixedsize=false
        width=1.9
        margin='0.14,0.10'
        fontsize=10
      ]
      pre1 [label='Communication\\nBetween-couple\\n(cbc)']
      out1 [label='Satisfaction\\n(Self)']
      pre2 [label='Communication\\nWithin-couple\\n(cwc)']

      edge [
        minlen=2.5
        labelfloat=true
        fontsize=12
      ]
      pre1 -> out1 [label='between']
      pre2 -> out1 [taillabel='within' labeldistance=4 labelangle=-15 constraint=false]
    }}
  "))
} else {
  # Native R graphics keep the PDF build independent of HTML-widget screenshot tools.
  plot.new()
  plot.window(xlim = c(0, 10), ylim = c(0, 10))
  title("Dyad–Individual Model (DIM)")

  draw_node <- function(x, y, label) {
    rect(x - 1.8, y - 0.9, x + 1.8, y + 0.9)
    text(x, y, label, cex = 0.78)
  }

  draw_node(2.2, 7, "Communication\nBetween-couple\n(cbc)")
  draw_node(2.2, 3, "Communication\nWithin-couple\n(cwc)")
  draw_node(7.8, 5, "Satisfaction\n(Self)")
  arrows(4.1, 6.7, 5.9, 5.6, length = 0.1)
  arrows(4.1, 3.3, 5.9, 4.4, length = 0.1)
  text(5, 6.5, "between", pos = 3)
  text(5, 3.5, "within", pos = 1)
}


```

---

## Exchangeable Dyads - Cross-Sectional DIM: Data

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

```{r data-exch-dim-source-preview, echo=FALSE, cache=TRUE, dependson='simulate-cross-sectional-data'}

df_exchangeable_dim <- prepare_dyad_data(
  data = df,
  dyad = coupleID,
  member = userID,
  role = gender,
  predictors = communication,
  model_types = "dim",
  keep_compositions = "female-male",
  set_exchangeable_compositions = "female-male",
  seed = 123
)

print_df(
  head(
    dplyr::select(
      df_exchangeable_dim,
      userID, coupleID, satisfaction, communication,
      .communication_dyad_mean_gmc, .communication_within_dyad_dev,
      .member_contrast_arbitrary
    )
  )
)

```

## 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

```{r data-exch-dim-preview, echo=FALSE}
print_df(
  head(
    dplyr::select(
      df_exchangeable_dim,
      userID, coupleID, satisfaction,
      .communication_dyad_mean_gmc, .communication_within_dyad_dev,
      .member_contrast_arbitrary
    )
  )
)
```


## Exchangeable Dyads - Cross-Sectional DIM: Model

```{r fit-exch-dim-brms, results='hide', message=FALSE, warning=FALSE, cache=TRUE, dependson='data-exch-dim-source-preview'}
#| code-line-numbers: "1-7|9-14|16-27"
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

```{r results-exch-dim-brms, echo=FALSE, cache=TRUE, dependson='fit-exch-dim-brms'}

summary_exch_dim <- summarize_brms(
    model_exch_dim,
    stats_to_report = c('CI', 'Rhat', 'ESS')
    )

print_df(
  summary_exch_dim,
  rows_to_pack = list(
    'Fixed Effects' = c(1,3),
    'Residual Structure' = c(4,5)
    )
)
```


## Equivalence APIM and DIM {.smaller}

`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 [@Bolger2025Unified]:

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

::: columns
::: column

```{r compare-exch-apim-dim-residual-apim, echo=FALSE}
print_df(
  summary_exch_apim[1:5, 1:2],
  rows_to_pack = list(
    'Fixed Effects' = c(1,3),
    'Residual Structure' = c(4,5)
    )
  )

```

:::

::: column

```{r compare-exch-apim-dim-residual-dim, echo=FALSE}
print_df(
  summary_exch_dim[1:5, 1:2],
  rows_to_pack = list(
    'Fixed Effects' = c(1,3),
    'Residual Structure' = c(4,5)
    )
  )

```

:::
:::


## Equivalence APIM and DIM
:::columns
:::column

```{r compare-exch-apim-dim-fixed-apim, echo=FALSE}
print_df(
  summary_exch_apim[1:3, 1:2],
  rows_to_pack = list(
    'Fixed Effects' = c(1,3)
    )
  )

```

:::

:::column

```{r compare-exch-apim-dim-fixed-dim, echo=FALSE}
print_df(
  summary_exch_dim[1:3, 1:2],
  rows_to_pack = list(
    'Fixed Effects' = c(1,3)
    )
  )

```

:::
:::

```{=html}

  <p style="margin-top:0;font-size:1.2rem;max-width:1100px;">
    Top sliders: grand-mean-centered <em>Communication</em> for <strong>Actor</strong> and <strong>Partner</strong> (APIM).<br>
    Bottom sliders: reparameterized communication — <strong>Centered Between-Couple</strong> (<code>x<sub>cbc</sub></code>) and <strong>Centered Within-Couple</strong> (<code>x<sub>cwc</sub></code>) (DIM).<br>
  </p>

  <style>
    .slider-grid{
      display:grid;
      grid-template-columns:250px 1fr 80px auto;
      gap:14px;
      align-items:center;
      font-size:1.2rem;
      margin-bottom:1.2rem;
      width:100%;
    }
    .slider-grid>label{grid-column:1;font-weight:600}
    .slider-grid>input[type="range"]{grid-column:2;width:100%;accent-color:steelblue;height:10px;cursor:pointer}
    .slider-grid>input[data-group="lower"]{accent-color:coral}
    .slider-grid>output{grid-column:3;font-variant-numeric:tabular-nums;text-align:right;width:3.8ch;font-size:1.15rem}
    .reset-cell{grid-column:4;grid-row:1/-1;align-self:center;justify-self:start}
    #reset-btn{background:#4B5563;color:#FAFAFA;border:none;padding:.6rem 1.2rem;border-radius:8px;font-size:1.05rem;cursor:pointer;transition:background .2s}
    #reset-btn:hover{background:#1F2937}

    .equations{margin-top:1.1rem}
    .math{
      font-family:'STIX Two Text','Times New Roman',serif;
      font-size:1.30rem;
      line-height:1.6;
      font-variant-numeric:lining-nums tabular-nums;
      letter-spacing:.05px;
    }
    .eq-table{
      display:grid;
      grid-template-columns:
        1.5ch auto 1.4ch minmax(7.8ch,max-content) 1.4ch minmax(19ch,max-content)
        1.4ch minmax(19ch,max-content) 1.4ch minmax(8.6ch,max-content);
      column-gap:.7rem; row-gap:.45rem; align-items:baseline; justify-content:start; width:100%;
    }
    .cell.label{font-weight:750;color:#334155;font-size:1.32rem}
    .cell.op{opacity:.9;text-align:center}

    .num{display:inline-flex;align-items:baseline}
    .sign{display:inline-block;width:1.1ch;text-align:right;padding-right:.18ch}
    .minus{visibility:visible}.phantom{visibility:hidden}
    .digits{display:inline-block;min-width:5.6ch;text-align:left}
    .times{padding:0 .24rem;font-weight:700;font-size:1.10em}

    .stack{display:inline-flex;flex-direction:column;align-items:center;gap:.08rem}
    .vlabel{font-size:.82rem;line-height:1;letter-spacing:.2px}

    .apx .digits,.apx .sign{color:steelblue}
    .dimx .digits,.dimx .sign{color:coral}
    .apx .vlabel{color:steelblue;font-weight:650}
    .dimx .vlabel{color:coral;font-weight:650}
    .coefwrap .vlabel{color:#6B7280;font-weight:500}
    .result .digits,.result .sign{color:#0f172a;font-weight:800}
  </style>

  <div class="slider-grid">
    <!-- Compact symbolic labels -->
    <label style="color:steelblue;">x<sub>actor_gmc</sub></label>
    <input id="slider-actor" type="range" min="-5" max="5" step="0.1"
      value="0"
      data-group="upper">
    <output id="out-actor">0.00</output>

    <label style="color:steelblue;">x<sub>partner_gmc</sub></label>
    <input id="slider-partner" type="range" min="-5" max="5" step="0.1"
      value="0"
      data-group="upper">
    <output id="out-partner">0.00</output>

    <label style="color:#E26A5C;">x<sub>cbc</sub></label>
    <input id="slider-mean" type="range" min="-5" max="5" step="0.1" value="0" data-group="lower">
    <output id="out-mean">0.00</output>

    <label style="color:#E26A5C;">x<sub>cwc</sub></label>
    <input id="slider-within" type="range" min="-5" max="5" step="0.1" value="0" data-group="lower">
    <output id="out-within">0.00</output>

    <div class="reset-cell"><button id="reset-btn">Reset to 0</button></div>
  </div>

  <!-- Model coefficients -->
  <div id="model-coefs"
       data-apim-b0="`r if (knitr::is_html_output()) round(apim_centered_intercept, 2) else ''`"
       data-apim-ba="`r if (knitr::is_html_output()) round(fixef(model_exch_apim)['.communication_actor','Estimate'], 2) else ''`"
       data-apim-bp="`r if (knitr::is_html_output()) round(fixef(model_exch_apim)['.communication_partner','Estimate'], 2) else ''`"
       data-dim-b0 ="`r if (knitr::is_html_output()) round(fixef(model_exch_dim)['Intercept','Estimate'], 2) else ''`"
       data-dim-cb ="`r if (knitr::is_html_output()) round(fixef(model_exch_dim)['.communication_dyad_mean_gmc','Estimate'], 2) else ''`"
       data-dim-cw ="`r if (knitr::is_html_output()) round(fixef(model_exch_dim)['.communication_within_dyad_dev','Estimate'], 2) else ''`"></div>

  <div class="equations">
    <div id="eq-table" class="eq-table math"></div>
  </div>

  <script>
(function(){
  const MIN=-5, MAX=5;
  const clamp=(x,lo,hi)=>Math.min(hi,Math.max(lo,x));
  const r1 = x => Math.round(x*10)/10;
  const r2 = x => { const y=Math.round(Number(x)*100)/100; return Math.abs(y)<1e-12?0:y; };

  // Elements
  const actor   = document.getElementById('slider-actor');   // x_actor_gmc
  const partner = document.getElementById('slider-partner'); // x_partner_gmc
  const mean    = document.getElementById('slider-mean');    // x_cbc
  const within  = document.getElementById('slider-within');  // x_cwc

  const outA = document.getElementById('out-actor');
  const outP = document.getElementById('out-partner');
  const outM = document.getElementById('out-mean');
  const outW = document.getElementById('out-within');

  const TIMES='×', MINUS='−';
  const coefEl=document.getElementById('model-coefs').dataset;
  const C = {
    apim:{ b0:+coefEl.apimB0, b_actor_gmc:+coefEl.apimBa, b_partner_gmc:+coefEl.apimBp },
    dim: { b0:+coefEl.dimB0,  b_cbc:+coefEl.dimCb,        b_cwc:+coefEl.dimCw }
  };

  const table=document.getElementById('eq-table');

  function SignedNumHTML(v, extra=''){
    const y=r2(v), neg=y<0, abs=Math.abs(y).toFixed(2);
    return `<span class="num ${extra}"><span class="sign ${neg?'minus':'phantom'}">${MINUS}</span><span class="digits">${abs}</span></span>`;
  }
  function VarNumber(label, value, colorClass){
    return `<span class="stack ${colorClass}"><span>${SignedNumHTML(value)}</span><span class="vlabel">${label}</span></span>`;
  }
  function CoefNumber(label, value){
    return `<span class="stack coefwrap"><span>${SignedNumHTML(value,'coef')}</span><span class="vlabel">${label}</span></span>`;
  }
  function row({modelLabel,b0Label,b0,b1Label,b1,x1Label,x1Val,x1Class,b2Label,b2,x2Label,x2Val,x2Class,yVal}){
    return `
      <div></div>
      <div class="cell label">${modelLabel}:&nbsp;Y</div>
      <div class="cell op">=</div>
      <div class="cell">${CoefNumber(b0Label,b0)}</div>
      <div class="cell op">+</div>
      <div class="cell">${CoefNumber(b1Label,b1)}<span class="times">${TIMES}</span>${VarNumber(x1Label,x1Val,x1Class)}</div>
      <div class="cell op">+</div>
      <div class="cell">${CoefNumber(b2Label,b2)}<span class="times">${TIMES}</span>${VarNumber(x2Label,x2Val,x2Class)}</div>
      <div class="cell op">=</div>
      <div class="cell result">${SignedNumHTML(yVal,'result')}</div>`;
  }

  function render(){
    const A=+actor.value, P=+partner.value, M=+mean.value, W=+within.value;
    const yAPIM = C.apim.b0 + C.apim.b_actor_gmc*A + C.apim.b_partner_gmc*P;
    const yDIM  = C.dim.b0  + C.dim.b_cbc*M        + C.dim.b_cwc*W;
    table.innerHTML =
      row({modelLabel:'APIM',
        b0Label:'b<sub>0</sub>', b0:C.apim.b0,
        b1Label:'b<sub>actor_gmc</sub>', b1:C.apim.b_actor_gmc,
        x1Label:'x<sub>actor_gmc</sub>', x1Val:A, x1Class:'apx',
        b2Label:'b<sub>partner_gmc</sub>', b2:C.apim.b_partner_gmc,
        x2Label:'x<sub>partner_gmc</sub>', x2Val:P, x2Class:'apx',
        yVal:yAPIM
      }) +
      row({modelLabel:'DIM',
        b0Label:'b<sub>0</sub>', b0:C.dim.b0,
        b1Label:'b<sub>cbc</sub>', b1:C.dim.b_cbc,
        x1Label:'x<sub>cbc</sub>', x1Val:M, x1Class:'dimx',
        b2Label:'b<sub>cwc</sub>', b2:C.dim.b_cwc,
        x2Label:'x<sub>cwc</sub>', x2Val:W, x2Class:'dimx',
        yVal:yDIM
      });
  }

  function updateOutputs(){
    outA.value=(+actor.value).toFixed(2);
    outP.value=(+partner.value).toFixed(2);
    outM.value=(+mean.value).toFixed(2);
    outW.value=(+within.value).toFixed(2);
  }

  /* =======================
     Feasible-interval helpers
     ======================= */
  // For moving A with P fixed:
  function feasibleA(p){
    // m=(a+p)/2 in [MIN,MAX]  -> a in [2*MIN - p, 2*MAX - p]
    // w=(a-p)/2 in [MIN,MAX]  -> a in [2*MIN + p, 2*MAX + p]
    const lo = Math.max(MIN, 2*MIN - p, 2*MIN + p);
    const hi = Math.min(MAX, 2*MAX - p, 2*MAX + p);
    return [lo, hi];
  }
  // For moving P with A fixed:
  function feasibleP(a){
    const lo = Math.max(MIN, 2*MIN - a, -2*MAX + a); // from w=(a-p)/2 -> p in [a-2*MAX, a-2*MIN]
    const hi = Math.min(MAX, 2*MAX - a, -2*MIN + a);
    // cleaner symmetrical form:
    return [Math.max(MIN, 2*MIN - a, a - 2*MAX),
            Math.min(MAX, 2*MAX - a, a - 2*MIN)];
  }
  // For moving M with W fixed:
  function feasibleM(w){
    // a=m+w ∈ [MIN,MAX] -> m ∈ [MIN - w, MAX - w]
    // p=m-w ∈ [MIN,MAX] -> m ∈ [MIN + w, MAX + w]
    const lo = Math.max(MIN, MIN - w, MIN + w);
    const hi = Math.min(MAX, MAX - w, MAX + w);
    return [lo, hi];
  }
  // For moving W with M fixed:
  function feasibleW(m){
    // a=m+w ∈ [MIN,MAX] -> w ∈ [MIN - m, MAX - m]
    // p=m-w ∈ [MIN,MAX] -> w ∈ [m - MAX, m - MIN]
    const lo = Math.max(MIN, MIN - m, m - MAX);
    const hi = Math.min(MAX, MAX - m, m - MIN);
    return [lo, hi];
  }

  /* =======================
     Move handlers (clamp the moved slider FIRST)
     ======================= */
  function moveActor(){
    const p = r1(+partner.value);
    const [lo,hi] = feasibleA(p);
    const a = r1(clamp(+actor.value, lo, hi)); // hard-stop the thumb
    actor.value = a;
    const m = r1((a + p)/2), w = r1((a - p)/2);
    mean.value = m; within.value = w;
    updateOutputs(); render();
  }
  function movePartner(){
    const a = r1(+actor.value);
    const [lo,hi] = feasibleP(a);
    const p = r1(clamp(+partner.value, lo, hi));
    partner.value = p;
    const m = r1((a + p)/2), w = r1((a - p)/2);
    mean.value = m; within.value = w;
    updateOutputs(); render();
  }
  function moveMean(){
    const w = r1(+within.value);
    const [lo,hi] = feasibleM(w);
    const m = r1(clamp(+mean.value, lo, hi));
    mean.value = m;
    const a = r1(m + w), p = r1(m - w);
    actor.value = a; partner.value = p;
    updateOutputs(); render();
  }
  function moveWithin(){
    const m = r1(+mean.value);
    const [lo,hi] = feasibleW(m);
    const w = r1(clamp(+within.value, lo, hi));
    within.value = w;
    const a = r1(m + w), p = r1(m - w);
    actor.value = a; partner.value = p;
    updateOutputs(); render();
  }

  // Listeners
  actor  .addEventListener('input', moveActor);
  partner.addEventListener('input', movePartner);
  mean   .addEventListener('input', moveMean);
  within .addEventListener('input', moveWithin);

  document.getElementById('reset-btn').addEventListener('click', ()=>{
    actor.value=partner.value=mean.value=within.value=0;
    updateOutputs(); render();
  });

updateOutputs(); render();
})();
</script>
```

## Equivalence APIM and DIM {data-background-color="#101820" .sumdiff-slide}


```{=html}

  <style>
    .reveal .slides section.sumdiff-slide{
      background:#101820;
      padding:22px 28px;
      border-radius:0;
    }
    .reveal .sumdiff-slide h2{
      color:#F8FAFC;
      margin-bottom:12px;
    }
    .sumdiff-v2{
      --actor:#4E94F8;
      --partner:#A78BFA;
      --cbc:#E26A5C;
      --cwc:#FBBF24;
      --ink:#101820;
      --panel:#17212f;
      --line:#314154;
      --text:#F8FAFC;
      --muted:#AAB6C6;
      color:var(--text);
      display:grid;
      grid-template-columns:minmax(260px,320px) minmax(0,1fr);
      gap:18px;
      align-items:stretch;
      width:min(100%,1000px);
      height:560px;
      margin:0 auto;
      font-size:15px;
    }
    .sumdiff-v2 *{box-sizing:border-box}
    .sumdiff-v2 h3{
      margin:0;
      color:#D6E4F0;
      font-size:.78rem;
      letter-spacing:.12em;
      text-transform:uppercase;
    }
    .sd2-panel{
      background:rgba(23,33,47,.92);
      border:1px solid rgba(170,182,198,.22);
      border-radius:12px;
      padding:14px;
      display:flex;
      flex-direction:column;
      gap:14px;
      min-width:0;
    }
    .sd2-head{
      display:flex;
      justify-content:space-between;
      align-items:center;
      gap:12px;
      border-bottom:1px solid rgba(170,182,198,.16);
      padding-bottom:10px;
    }
    #sd2-reset{
      border:1px solid rgba(170,182,198,.25);
      background:#263447;
      color:#F8FAFC;
      padding:7px 11px;
      border-radius:6px;
      font-size:.72rem;
      text-transform:uppercase;
      letter-spacing:.08em;
      cursor:pointer;
    }
    .sd2-group{
      display:flex;
      flex-direction:column;
      gap:10px;
    }
    .sd2-group-title{
      display:flex;
      align-items:center;
      gap:8px;
      color:var(--muted);
      font-weight:700;
      text-transform:uppercase;
      letter-spacing:.08em;
      font-size:.72rem;
      border-bottom:1px solid rgba(170,182,198,.14);
      padding-bottom:7px;
    }
    .sd2-dot{width:9px;height:9px;border-radius:999px;display:inline-block}
    .sd2-control{display:grid;grid-template-columns:1fr 58px;gap:8px 12px;align-items:center}
    .sd2-control label{font-size:.78rem;font-weight:700;color:#D6E4F0}
    .sd2-control output{
      justify-self:end;
      font-family:ui-monospace,SFMono-Regular,Menlo,monospace;
      font-size:.88rem;
      font-weight:800;
      min-width:52px;
      text-align:right;
    }
    .sd2-control input{grid-column:1 / 3;width:100%;cursor:pointer}
    .sd2-control-actor input{accent-color:var(--actor)}
    .sd2-control-partner input{accent-color:var(--partner)}
    .sd2-control-cbc input{accent-color:var(--cbc)}
    .sd2-control-cwc input{accent-color:var(--cwc)}
    .sd2-relations{
      border-top:1px solid rgba(170,182,198,.14);
      padding-top:12px;
      color:var(--muted);
      font-size:.74rem;
      line-height:1.45;
    }
    .sd2-main{
      display:flex;
      flex-direction:column;
      gap:14px;
      min-width:0;
    }
    .sd2-equations{
      background:rgba(23,33,47,.74);
      border:1px solid rgba(170,182,198,.18);
      border-radius:12px;
      padding:11px 14px;
      display:grid;
      gap:7px;
      font-family:'STIX Two Text','Times New Roman',serif;
      font-size:.98rem;
      font-variant-numeric:tabular-nums;
    }
    .sd2-equations > div{white-space:nowrap}
    .sd2-eq-label{
      font-family:system-ui,-apple-system,BlinkMacSystemFont,'Segoe UI',sans-serif;
      font-size:.64rem;
      letter-spacing:.12em;
      text-transform:uppercase;
      color:var(--muted);
      font-weight:800;
      margin-right:8px;
    }
    .sd2-var{font-style:italic}
    .sd2-var sub,
    .sd2-beta sub{
      font-style:normal;
      font-size:.68em;
    }
    .sd2-num{
      display:inline-block;
      min-width:5.4ch;
      text-align:right;
      font-variant-numeric:tabular-nums;
    }
    .sd2-result{
      display:inline-block;
      min-width:5.7ch;
      text-align:right;
      font-weight:800;
      font-variant-numeric:tabular-nums;
    }
    .sd2-sign-pad{
      display:inline-block;
      width:1ch;
    }
    .sd2-times{padding:0 .08em}
    .sd2-grid-wrap{
      display:grid;
      grid-template-columns:minmax(320px,380px) minmax(210px,260px);
      gap:16px;
      align-items:start;
    }
    #sd2-svg{
      width:100%;
      aspect-ratio:1 / 1;
      background:rgba(15,23,42,.46);
      border:1px solid rgba(170,182,198,.18);
      border-radius:14px;
      cursor:crosshair;
      touch-action:none;
    }
    .sd2-explain{
      background:rgba(23,33,47,.74);
      border:1px solid rgba(170,182,198,.18);
      border-radius:12px;
      padding:13px;
      color:#D6E4F0;
      font-size:.76rem;
      line-height:1.36;
      display:flex;
      flex-direction:column;
      gap:10px;
    }
    .sd2-explain strong{color:#fff}
    .sd2-demo{
      border:1px solid rgba(214,228,240,.22);
      background:rgba(214,228,240,.08);
      color:#F8FAFC;
      border-radius:6px;
      padding:4px 8px;
      font-size:.65rem;
      text-transform:uppercase;
      letter-spacing:.08em;
      font-weight:800;
      cursor:pointer;
      margin-left:6px;
    }
    .sd2-mini-eq{
      color:#F8FAFC;
      background:rgba(15,23,42,.45);
      border:1px solid rgba(170,182,198,.12);
      border-radius:8px;
      padding:8px 10px;
      font-family:'STIX Two Text','Times New Roman',serif;
      text-align:center;
    }
  </style>

  <div class="sumdiff-v2">
    <div class="sd2-panel">
      <div class="sd2-head">
        <h3>Coordinate Inputs</h3>
          <button id="sd2-reset" type="button">Reset</button>
      </div>

      <div class="sd2-group">
        <div class="sd2-group-title">
          <span class="sd2-dot" style="background:var(--actor)"></span>
          <span class="sd2-dot" style="background:var(--partner)"></span>
          APIM: individual coordinates
        </div>
        <div class="sd2-control sd2-control-actor">
          <label for="sd2-actor">Actor predictor, <span class="sd2-var">x<sub>actor</sub></span></label>
          <output id="sd2-out-actor" style="color:var(--actor)">0.0</output>
          <input id="sd2-actor" type="range" min="-5" max="5" step="0.1" value="0">
        </div>
        <div class="sd2-control sd2-control-partner">
          <label for="sd2-partner">Partner predictor, <span class="sd2-var">x<sub>partner</sub></span></label>
          <output id="sd2-out-partner" style="color:var(--partner)">0.0</output>
          <input id="sd2-partner" type="range" min="-5" max="5" step="0.1" value="0">
        </div>
      </div>

      <div class="sd2-group">
        <div class="sd2-group-title">
          <span class="sd2-dot" style="background:var(--cbc)"></span>
          <span class="sd2-dot" style="background:var(--cwc)"></span>
          DIM: mean and deviation coordinates
        </div>
        <div class="sd2-control sd2-control-cbc">
          <label for="sd2-cbc">Couple mean, <span class="sd2-var">x<sub>cbc</sub></span></label>
          <output id="sd2-out-cbc" style="color:var(--cbc)">0.0</output>
          <input id="sd2-cbc" type="range" min="-5" max="5" step="0.1" value="0">
        </div>
        <div class="sd2-control sd2-control-cwc">
          <label for="sd2-cwc">Couple deviation, <span class="sd2-var">x<sub>cwc</sub></span></label>
          <output id="sd2-out-cwc" style="color:var(--cwc)">0.0</output>
          <input id="sd2-cwc" type="range" min="-5" max="5" step="0.1" value="0">
        </div>
      </div>

      <div class="sd2-relations">
        <div><span class="sd2-var">x<sub>cbc</sub></span> = (<span class="sd2-var">x<sub>actor</sub></span> + <span class="sd2-var">x<sub>partner</sub></span>) / 2</div>
        <div><span class="sd2-var">x<sub>cwc</sub></span> = (<span class="sd2-var">x<sub>actor</sub></span> − <span class="sd2-var">x<sub>partner</sub></span>) / 2</div>
      </div>
    </div>

    <main class="sd2-main">
      <div class="sd2-equations">
        <div><span class="sd2-eq-label">APIM</span><span id="sd2-eq-apim"></span></div>
        <div><span class="sd2-eq-label">DIM</span><span id="sd2-eq-dim"></span></div>
        <div><span class="sd2-eq-label">Coefficients</span><span id="sd2-eq-coef"></span></div>
      </div>

      <div class="sd2-grid-wrap">
        <svg id="sd2-svg" viewBox="0 0 320 320" role="img" aria-label="APIM and DIM coordinate grid">
          <defs>
            <clipPath id="sd2-clip"><rect x="30" y="30" width="260" height="260" rx="10"></rect></clipPath>
          </defs>
          <rect x="30" y="30" width="260" height="260" rx="10" fill="#0F172A" fill-opacity=".52" stroke="#64748B" stroke-width="1.5"></rect>
          <g id="sd2-grid-lines" clip-path="url(#sd2-clip)"></g>
          <g id="sd2-axis-labels"></g>
          <circle id="sd2-halo" cx="160" cy="160" r="13" fill="white" fill-opacity=".12"></circle>
          <circle id="sd2-point" cx="160" cy="160" r="6" fill="#4E94F8" stroke="#fff" stroke-width="2"></circle>
        </svg>

        <div class="sd2-explain">
          <div>
            <strong>Between-couple effect: the mean.</strong>
            If both partners increase by 1, the couple mean increases by 1 and the deviation stays 0.
            <button class="sd2-demo" id="sd2-demo-mean" type="button">Demonstrate</button>
            <div class="sd2-mini-eq"><span class="sd2-beta">b<sub>cbc</sub></span> = <span class="sd2-beta">b<sub>actor</sub></span> + <span class="sd2-beta">b<sub>partner</sub></span></div>
          </div>
          <div>
            <strong>Within-couple effect: the deviation.</strong>
            If the actor increases by 1 and the partner decreases by 1, the couple mean stays 0 and the deviation increases by 1.
            <button class="sd2-demo" id="sd2-demo-deviation" type="button">Demonstrate</button>
            <div class="sd2-mini-eq"><span class="sd2-beta">b<sub>cwc</sub></span> = <span class="sd2-beta">b<sub>actor</sub></span> − <span class="sd2-beta">b<sub>partner</sub></span></div>
          </div>
          <div>
            The dot lives in APIM space: horizontal = actor predictor, vertical = partner predictor. The diagonal axes show the same point in DIM space.
          </div>
        </div>
      </div>
    </main>
  </div>

  <div id="sd2-coefs"
       data-b0="`r if (knitr::is_html_output()) round(apim_centered_intercept, 2) else ''`"
       data-ba="`r if (knitr::is_html_output()) round(fixef(model_exch_apim)['.communication_actor','Estimate'], 2) else ''`"
       data-bp="`r if (knitr::is_html_output()) round(fixef(model_exch_apim)['.communication_partner','Estimate'], 2) else ''`"></div>

  <script>
  (function(){
    const root = document.currentScript.closest('section');
    const MIN = -5, MAX = 5, PLOT_DOMAIN = 5;
    const c = {actor:'#4E94F8', partner:'#A78BFA', cbc:'#E26A5C', cwc:'#FBBF24'};
    const coefData = root.querySelector('#sd2-coefs').dataset;
    const b0 = Number(coefData.b0), bA = Number(coefData.ba), bP = Number(coefData.bp);
    const bCbc = +(bA + bP).toFixed(2);
    const bCwc = +(bA - bP).toFixed(2);

    const el = id => root.querySelector('#' + id);
    const actor = el('sd2-actor'), partner = el('sd2-partner'), cbc = el('sd2-cbc'), cwc = el('sd2-cwc');
    const outA = el('sd2-out-actor'), outP = el('sd2-out-partner'), outCbc = el('sd2-out-cbc'), outCwc = el('sd2-out-cwc');
    const eqA = el('sd2-eq-apim'), eqD = el('sd2-eq-dim'), eqC = el('sd2-eq-coef');
    const point = el('sd2-point'), halo = el('sd2-halo'), svg = el('sd2-svg');

    const r1 = x => Math.round(Number(x) * 10) / 10;
    const fmt = x => (Math.abs(x) < 1e-10 ? 0 : x).toFixed(2);
    const mathNum = (x, color = null, result = false) => {
      const value = Math.abs(x) < 1e-10 ? 0 : x;
      const sign = value < 0 ? '−' : '<span class="sd2-sign-pad" aria-hidden="true"></span>';
      const style = color ? ` style="color:${color}"` : '';
      const cls = result ? 'sd2-result' : 'sd2-num';
      return `<span class="${cls}"${style}>${sign}${Math.abs(value).toFixed(2)}</span>`;
    };
    const beta = name => `<span class="sd2-beta">b<sub>${name}</sub></span>`;
    const xPix = x => 160 + x * (130 / PLOT_DOMAIN);
    const yPix = y => 160 - y * (130 / PLOT_DOMAIN);
    const clamp = (x, lo, hi) => Math.min(hi, Math.max(lo, Number(x)));

    function feasibleA(p){
      return [Math.max(MIN, 2 * MIN - p, 2 * MIN + p),
              Math.min(MAX, 2 * MAX - p, 2 * MAX + p)];
    }
    function feasibleP(a){
      return [Math.max(MIN, 2 * MIN - a, a - 2 * MAX),
              Math.min(MAX, 2 * MAX - a, a - 2 * MIN)];
    }
    function feasibleM(w){
      return [Math.max(MIN, MIN - w, MIN + w),
              Math.min(MAX, MAX - w, MAX + w)];
    }
    function feasibleW(m){
      return [Math.max(MIN, MIN - m, m - MAX),
              Math.min(MAX, MAX - m, m - MIN)];
    }

    function setFromAP(a, p){
      const aClamped = r1(clamp(a, MIN, MAX));
      const pClamped = r1(clamp(p, MIN, MAX));
      actor.value = aClamped;
      partner.value = pClamped;
      cbc.value = r1((aClamped + pClamped) / 2);
      cwc.value = r1((aClamped - pClamped) / 2);
      render();
    }
    function setFromDIM(m, w){
      const mClamped = r1(clamp(m, ...feasibleM(w)));
      const wClamped = r1(clamp(w, ...feasibleW(mClamped)));
      cbc.value = mClamped;
      cwc.value = wClamped;
      actor.value = r1(mClamped + wClamped);
      partner.value = r1(mClamped - wClamped);
      render();
    }
    function moveActor(){
      const p = r1(partner.value);
      const bounds = feasibleA(p);
      const a = r1(clamp(actor.value, bounds[0], bounds[1]));
      setFromAP(a, p);
    }
    function movePartner(){
      const a = r1(actor.value);
      const bounds = feasibleP(a);
      const p = r1(clamp(partner.value, bounds[0], bounds[1]));
      setFromAP(a, p);
    }
    function moveCbc(){
      const w = r1(cwc.value);
      const bounds = feasibleM(w);
      const m = r1(clamp(cbc.value, bounds[0], bounds[1]));
      setFromDIM(m, w);
    }
    function moveCwc(){
      const m = r1(cbc.value);
      const bounds = feasibleW(m);
      const w = r1(clamp(cwc.value, bounds[0], bounds[1]));
      setFromDIM(m, w);
    }

    function render(){
      const a = Number(actor.value), p = Number(partner.value), m = Number(cbc.value), w = Number(cwc.value);
      outA.value = a.toFixed(1); outP.value = p.toFixed(1); outCbc.value = m.toFixed(2); outCwc.value = w.toFixed(2);
      const yA = b0 + bA * a + bP * p;
      const yD = b0 + bCbc * m + bCwc * w;
      eqA.innerHTML = `<span class="sd2-var">Ŷ</span> = ${mathNum(b0)} + ${mathNum(bA)} <span class="sd2-times">×</span> ${mathNum(a, c.actor)} + ${mathNum(bP)} <span class="sd2-times">×</span> ${mathNum(p, c.partner)} = ${mathNum(yA, null, true)}`;
      eqD.innerHTML = `<span class="sd2-var">Ŷ</span> = ${mathNum(b0)} + ${mathNum(bCbc)} <span class="sd2-times">×</span> ${mathNum(m, c.cbc)} + ${mathNum(bCwc)} <span class="sd2-times">×</span> ${mathNum(w, c.cwc)} = ${mathNum(yD, null, true)}`;
      eqC.innerHTML = `${beta('cbc')} = ${mathNum(bA)} + ${mathNum(bP)} = ${mathNum(bCbc)}; &nbsp; ${beta('cwc')} = ${mathNum(bA)} − ${mathNum(bP)} = ${mathNum(bCwc)}`;
      point.setAttribute('cx', xPix(a)); point.setAttribute('cy', yPix(p));
      halo.setAttribute('cx', xPix(a)); halo.setAttribute('cy', yPix(p));
    }

    function drawGrid(){
      const g = el('sd2-grid-lines');
      const ticks = [-4, -2, 0, 2, 4];
      let html = '';
      ticks.forEach(t => {
        const main = t === 0 ? 1.5 : 1;
        html += `<line x1="${xPix(t)}" y1="30" x2="${xPix(t)}" y2="290" stroke="${c.actor}" stroke-opacity="${t===0?.48:.20}" stroke-width="${main}" stroke-dasharray="${t===0?'':'2 2'}"></line>`;
        html += `<line x1="30" y1="${yPix(t)}" x2="290" y2="${yPix(t)}" stroke="${c.partner}" stroke-opacity="${t===0?.48:.20}" stroke-width="${main}" stroke-dasharray="${t===0?'':'2 2'}"></line>`;
        html += `<line x1="${xPix(-PLOT_DOMAIN)}" y1="${yPix(2*t + PLOT_DOMAIN)}" x2="${xPix(PLOT_DOMAIN)}" y2="${yPix(2*t - PLOT_DOMAIN)}" stroke="${c.cbc}" stroke-opacity="${t===0?.50:.22}" stroke-width="${main}" stroke-dasharray="${t===0?'':'3 3'}"></line>`;
        html += `<line x1="${xPix(-PLOT_DOMAIN)}" y1="${yPix(-PLOT_DOMAIN - 2*t)}" x2="${xPix(PLOT_DOMAIN)}" y2="${yPix(PLOT_DOMAIN - 2*t)}" stroke="${c.cwc}" stroke-opacity="${t===0?.50:.22}" stroke-width="${main}" stroke-dasharray="${t===0?'':'3 3'}"></line>`;
      });
      g.innerHTML = html;
      el('sd2-axis-labels').innerHTML = `
        <text x="286" y="148" fill="${c.actor}" font-size="12" font-style="italic" font-weight="800" text-anchor="end">x<tspan baseline-shift="sub" font-size="8.5" font-style="normal">actor</tspan></text>
        <text x="150" y="22" fill="${c.partner}" font-size="12" font-style="italic" font-weight="800" text-anchor="middle">x<tspan baseline-shift="sub" font-size="8.5" font-style="normal">partner</tspan></text>
        <g transform="translate(${xPix(3.4)},${yPix(3.4)}) rotate(-45)"><text y="-7" fill="${c.cbc}" font-size="11" font-style="italic" font-weight="800" text-anchor="middle">x<tspan baseline-shift="sub" font-size="8" font-style="normal">cbc</tspan></text></g>
        <g transform="translate(${xPix(3.4)},${yPix(-3.4)}) rotate(45)"><text y="-7" fill="${c.cwc}" font-size="11" font-style="italic" font-weight="800" text-anchor="middle">x<tspan baseline-shift="sub" font-size="8" font-style="normal">cwc</tspan></text></g>`;
    }

    function gridToAP(evt){
      const rect = svg.getBoundingClientRect();
      const clientX = evt.touches ? evt.touches[0].clientX : evt.clientX;
      const clientY = evt.touches ? evt.touches[0].clientY : evt.clientY;
      const sx = ((clientX - rect.left) / rect.width) * 320;
      const sy = ((clientY - rect.top) / rect.height) * 320;
      setFromAP(r1((sx - 160) / (130 / PLOT_DOMAIN)), r1((160 - sy) / (130 / PLOT_DOMAIN)));
    }
    let dragging = false;
    svg.addEventListener('mousedown', evt => { dragging = true; gridToAP(evt); });
    svg.addEventListener('mousemove', evt => { if (dragging) gridToAP(evt); });
    window.addEventListener('mouseup', () => { dragging = false; });
    svg.addEventListener('touchstart', evt => { dragging = true; gridToAP(evt); }, {passive:true});
    svg.addEventListener('touchmove', evt => { if (dragging) gridToAP(evt); }, {passive:true});
    window.addEventListener('touchend', () => { dragging = false; });

    actor.addEventListener('input', moveActor);
    partner.addEventListener('input', movePartner);
    cbc.addEventListener('input', moveCbc);
    cwc.addEventListener('input', moveCwc);
    function bindButton(id, fn){
      const button = el(id);
      button.addEventListener('pointerdown', evt => evt.stopPropagation());
      button.addEventListener('click', evt => {
        evt.preventDefault();
        evt.stopPropagation();
        fn();
      });
    }
    bindButton('sd2-reset', () => setFromAP(0, 0));
    bindButton('sd2-demo-mean', () => setFromAP(1, 1));
    bindButton('sd2-demo-deviation', () => setFromAP(1, -1));

    drawGrid();
    setFromAP(0, 0);
  })();
  </script>
```

::: {.content-visible when-format="pdf"}

In the interactive HTML slides, this section contains sliders that demonstrate
the APIM–DIM reparameterization. In this static version:

* actor and partner predictors are rewritten as a dyad-level mean and a
  within-dyad half-difference; and
* the DIM coefficients are a linear reparameterization of the APIM coefficients,
  so fitted values are unchanged.

```text
x_cbc = (x_actor_gmc + x_partner_gmc) / 2
x_cwc = (x_actor_gmc - x_partner_gmc) / 2
```

:::


## 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}$$

[@Bolger2025Unified]

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

```{r compare-exch-apim-dim-hypothesis, cache=TRUE, dependson='fit-exch-apim-brms'}
a <- hypothesis(
  model_exch_apim,
  ".communication_actor + .communication_partner = 0"
)

round(a$hypothesis[,c(2,3,4,5)], 2)
```

:::columns
:::column

```{r compare-exch-apim-hyp-apim, echo=FALSE}
print_df(
  summary_exch_apim[1:3, 1:2],
  rows_to_pack = list(
    'Fixed Effects' = c(1,3)
    )
  )

```

:::

:::column

```{r compare-exch-apim-hyp-dim, echo=FALSE}
print_df(
  summary_exch_dim[1:3, 1:2],
  rows_to_pack = list(
    'Fixed Effects' = c(1,3)
    )
  )

```

:::
:::

## 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., @stadlerHowCompanionshipRelated2023] is equivalent to the corresponding full APIM (aside from prior impact in Bayesian fits) [@Bolger2025Unified; @iidaModelsInterdependentIndividuals2018].



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

```{r data-exch-long-apim-source, results='hide', cache=TRUE, dependson='simulate-longitudinal-data'}
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)
  )

```


```{r data-exch-long-apim-preview, echo=FALSE}

print_df(
  print_couple_preview_for_output(
    df_long_dim,
    couple_id = "31",
    days = if (knitr::is_latex_output()) 0:2 else 0:3,
    pdf_cols = c(
      "userID", ".member_contrast_arbitrary", "diaryday", "diaryday_c",
      "closeness", ".provided_support_cwp_actor",
      ".provided_support_cwp_partner"
    )
  ),
  font_size = if (knitr::is_latex_output()) 6 else 20
)

```


# 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 [@bolgerIntensiveLongitudinalMethods2013].

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

::: {.content-visible when-format="pdf"}

```{r data-exch-long-dim-preview, echo=FALSE, cache=TRUE, dependson='data-exch-long-apim-source', eval=knitr::is_latex_output()}
print_df(
  print_couple_preview_for_output(
    df_long_dim,
    couple_id = "31",
    pdf_cols = c(
      "userID", ".member_contrast_arbitrary", "diaryday", "diaryday_c",
      ".provided_support_cbp_dyad_mean",
      ".provided_support_cbp_within_dyad_dev"
    )
  ),
  font_size = 6
)
```

:::


::: {.content-visible when-format="pdf"}

```{=latex}
\Needspace{20\baselineskip}
```

:::


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

```{r data-exch-long-member-contrast, echo=FALSE, cache=TRUE, dependson='data-exch-long-apim-source'}
print_df(
  print_couple_preview_for_output(
    df_long_dim,
    couple_id = "31",
    days = 0:4,
    pdf_cols = c(
      "userID", ".member_contrast_arbitrary", "diaryday", "diaryday_c",
      ".provided_support_cwp_dyad_mean",
      ".provided_support_cwp_within_dyad_dev"
    )
  ),
  font_size = if (knitr::is_latex_output()) 6 else 20
)

```


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

```{r fit-exch-long-apim-glmmtmb, eval=FALSE, results='hide', message=FALSE, warning=FALSE}
#| code-line-numbers: "1-4|6-11|13-18|20-21|23-28|30-34|36-39"

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

```{r results-exch-long-apim-glmmtmb, eval=FALSE, echo=FALSE, class.output='model-summary-output'}
summary(model_exch_long_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:

::: {.brms-prior-note}
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.
:::

```{r fit-exch-long-apim-brms, results='hide', message=FALSE, warning=FALSE, cache=TRUE, dependson='data-exch-long-member-contrast'}
#| code-line-numbers: "1-4|6-12|14-15|17-22|24-25|27-29"

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:

```{r fit-exch-long-dim-brms, results='hide', message=FALSE, warning=FALSE, cache=TRUE, dependson='data-exch-long-member-contrast'}
#| code-line-numbers: "7-11|13-17| 19-20"

# 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 [@delrosarioPracticalGuideSpecifying2025; @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.

(@delrosarioPracticalGuideSpecifying2025 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 [@delrosarioPracticalGuideSpecifying2025].

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

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

:::columns
:::column

```{r results-exch-long-apim-fixed, echo = FALSE, cache=TRUE, dependson='fit-exch-long-apim-brms'}
summary_exch_long_apim <- summarize_brms(
    model_exch_long_apim,
    stats_to_report = c('CI', 'Rhat', 'ESS')
    )

print_df(
  summary_exch_long_apim[1:6,1:2],
  rows_to_pack = list(
    'Fixed Effects' = c(1,6)
    )
  )

```

:::

:::column

```{r results-exch-long-dim-fixed, echo = FALSE, cache=TRUE, dependson='fit-exch-long-dim-brms'}
summary_exch_long_dim <- summarize_brms(
    model_exch_long_dim,
    stats_to_report = c('CI', 'Rhat', 'ESS')
    )

print_df(
  summary_exch_long_dim[1:6,1:2],
  rows_to_pack = list(
    'Fixed Effects' = c(1,6)
    )
  )

```

:::
:::

- Equivalence holds on both levels. Numeric deviations arise due to priors, sampling noise, and rounding.
  

::: {.content-visible when-format="pdf"}

```{=latex}
\begin{landscape}
```

:::

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

APIM Results in Detail:

```{r results-exch-long-apim-detail, echo=FALSE, cache=TRUE, dependson='fit-exch-long-apim-brms'}
summary_exch_long_apim <- summarize_brms(
    model_exch_long_apim,
    stats_to_report = c('CI', 'Rhat', 'ESS')
    )

print_df(
  summary_exch_long_apim,
  rows_to_pack = list(
    'Fixed Effects' = c(1,6),
    'Random Effects' = c(7, 27),
    'Residual Structure' = c(28, 28)
    ),
  scroll_height = '50vh'
)

```

::: {.content-visible when-format="pdf"}

```{=latex}
\end{landscape}
```

:::


## 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

```{r extract-exch-long-apim-covariance, cache=TRUE, dependson='fit-exch-long-apim-brms'}
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

```{r results-exch-long-apim-cov-within}
print_df(a$within_person_covariance_matrix)
```


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

```{r results-exch-long-apim-cov-cross}
print_df(a$cross_person_covariance_matrix)
```


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

```{r results-exch-long-apim-cov-full}
print_wide_matrix(a$full_covariance_matrix, font_size = 6)
```


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

```{r results-exch-long-apim-re-sds}
print_df(
  a$sd_summary %>%
    mutate(across(where(is.numeric), ~ round(.x, 3)))
)
```

## Random effect correlation matrix

```{r results-exch-long-apim-re-cor}
print_wide_matrix(a$full_correlation_matrix, font_size = 6)
```

## Full reporting summary

```{r results-exch-long-apim-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'
)
```


## References

```{r write-package-bibliography, echo=FALSE, results='hide'}

pkgs <- c(
  "brms", "bmlm", "easystats", "DHARMa", "tidyverse", "cmdstanr",
  "loo", "rstan", "MASS", "nlme", "glmmTMB"
)

knitr::write_bib(pkgs, "Embed/Rpackages.bib")

```