---
title: "Actor-Partner Interdependence Model (APIM)"
author: "Pascal Küng"
bibliography: references.bib
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Actor-Partner Interdependence Model (APIM)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

diagram_helper <- if (file.exists("diagram-helpers.Rinc")) {
  "diagram-helpers.Rinc"
} else {
  file.path("vignettes", "diagram-helpers.Rinc")
}
sys.source(diagram_helper, envir = knitr::knit_global())

diagram_device <- .vignette_diagram_device()
knitr::opts_chunk$set(dev = diagram_device, fig.ext = diagram_device)
```

```{r setup}
library(dyadMLM)
has_glmmTMB <- requireNamespace("glmmTMB", quietly = TRUE)
apim_distinguishable_fitted_alt <-
  "Fitted distinguishable APIM diagram unavailable."
apim_exchangeable_fitted_alt <-
  "Fitted exchangeable APIM diagram unavailable."
```

This vignette focuses on Gaussian cross-sectional and intensive longitudinal
Actor-Partner Interdependence models for distinguishable and exchangeable
dyads. The intensive longitudinal examples cover concurrent associations and a
simple dynamic model of stability and partner influence.

For the broader package workflow and an overview of the available
model-specific vignettes, including the
[Dyad-Individual Model](dim.html) and
[Dyadic Score Model](dsm.html), see the
[online package overview](https://pascal-kueng.github.io/dyadMLM/).

A vignette for non-Gaussian generalized models is planned.


# Cross-sectional APIMs

## The distinguishable APIM

A conceptual example for distinguishable female-male dyads:

```{r distinguishable-apim-diagram, echo=FALSE, fig.width=9, fig.height=4.8, out.width="100%", fig.cap="Conceptual cross-sectional APIM for distinguishable female-male dyads. Intercepts $b_\\mathrm{0}$, actor effects $a$, and partner effects $p$ can differ by the role of the outcome member (F and M), and the two outcome residuals covary within dyads.", fig.alt="Path diagram for a distinguishable cross-sectional APIM. Female and male outcomes have separate intercepts. Female and male predictors each have an actor path to their own outcome and a partner path to the other member's outcome. The female and male outcome residuals covary."}
draw_apim_diagram("distinguishable")
```

For univariate MLM software like `glmmTMB`, this model is fitted in long format 
with one outcome row per member, which can be visualized as:

```{r distinguishable-apim-member-diagram, echo=FALSE, fig.width=9, fig.height=5.2, out.width="100%", fig.cap="Individual-level representation of the distinguishable cross-sectional APIM used for the long-format multilevel model. For the female outcome, the female predictor is the actor predictor and the male predictor is the partner predictor. These roles reverse for the male outcome. Intercepts, actor coefficients, and partner coefficients may differ by outcome role, and the two member residuals may have different variances and covary.", fig.alt="Two-panel path diagram for a distinguishable female-male APIM. The female and male outcomes have separate intercepts. In the female outcome panel, female X is the actor predictor and male X is the partner predictor of female Y, with coefficients a F and p F. In the male outcome panel, male X is the actor predictor and female X is the partner predictor of male Y, with coefficients a M and p M. The female and male outcome residuals covary."}
draw_apim_member_diagram("distinguishable")
```


### Residual random-effects structure {#distinguishable-residual-structure}

For a distinguishable female-male dyad, the two members can have different
residual variances. The within-dyad residual covariance block (shared across dyads) is:

$$
\operatorname{Cov}
\begin{pmatrix}
\epsilon_{Fi} \\
\epsilon_{Mi}
\end{pmatrix}
= \boldsymbol{\Sigma}_{\epsilon}
= \begin{bmatrix}
\sigma_{\epsilon_F}^{2}
& \rho_{\epsilon_F\epsilon_M}\sigma_{\epsilon_F}\sigma_{\epsilon_M} \\
\rho_{\epsilon_F\epsilon_M}\sigma_{\epsilon_F}\sigma_{\epsilon_M}
& \sigma_{\epsilon_M}^{2}
\end{bmatrix}
$$

And the full residual covariance matrix for all dyads (first three shown) is then block-diagonal:

$$
\boldsymbol{\Sigma}_{\mathrm{model}}
= \begin{bmatrix}
\boldsymbol{\Sigma}_{\epsilon} & \boldsymbol{0} & \boldsymbol{0} & \cdots \\
\boldsymbol{0} & \boldsymbol{\Sigma}_{\epsilon} & \boldsymbol{0} & \cdots \\
\boldsymbol{0} & \boldsymbol{0} & \boldsymbol{\Sigma}_{\epsilon} & \cdots \\
\vdots & \vdots & \vdots & \ddots
\end{bmatrix}
$$


This structure is estimated with an unstructured random-effects block such as
`us(0 + .dy_is_female_x_male_female + .dy_is_female_x_male_male | coupleID)`
and `dispformula = ~ 0`. 

### Fitting the distinguishable APIM with glmmTMB

We first prepare the example data with `dyadMLM::prepare_dyad_data()`:

```{r prepare-cross-distinguishable-apim}
apim_distinguishable_data <- dyadMLM::prepare_dyad_data(
  data = dyads_cross,
  dyad = coupleID,
  member = personID,
  role = gender,
  predictors = provided_support,
  model_types = "apim",
  # All three observed compositions in `dyads_cross` are detected and retained by
  # default. This example focuses on `female-male` dyads, so we restrict the
  # analysis here.
  keep_compositions = "female-male"
)

print(apim_distinguishable_data, n=4)

```

The generated `.dy_*` columns can be used directly in the model formula. Here is
a simple example:

```{r fit-cross-distinguishable-gaussian, eval = has_glmmTMB}
apim_distinguishable_model <- glmmTMB::glmmTMB(
  closeness ~

    # Gender-specific intercepts
    0 +
    .dy_is_female_x_male_female +
    .dy_is_female_x_male_male +

    # Gender-specific actor effects
    .dy_is_female_x_male_female:.dy_provided_support_actor +
    .dy_is_female_x_male_male:.dy_provided_support_actor +

    # Gender-specific partner effects
    .dy_is_female_x_male_female:.dy_provided_support_partner +
    .dy_is_female_x_male_male:.dy_provided_support_partner +

    # Dyad-level unstructured random effects represent the two partner
    # residual variances and their covariance when dispformula = ~ 0.
    # This is glmmTMB-specific syntax! `brms` uses different syntax.
    us(0 +
         .dy_is_female_x_male_female +
         .dy_is_female_x_male_male
       | coupleID)

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

summary(apim_distinguishable_model)
```



```{r prepare-fitted-distinguishable-apim-diagram, include=FALSE, eval=has_glmmTMB}
apim_distinguishable_diagram_values <- .extract_apim_diagram_values(
  apim_distinguishable_model,
  type = "distinguishable"
)
apim_distinguishable_diagram_estimates <-
  apim_distinguishable_diagram_values$estimates
apim_distinguishable_diagram_residuals <-
  apim_distinguishable_diagram_values$residuals

apim_distinguishable_fitted_alt <- sprintf(
  paste(
    "Fitted distinguishable APIM. Female and male intercepts %.2f and %.2f;",
    "actor effects %.2f and %.2f; partner effects %.2f and %.2f;",
    "residual SDs %.2f and %.2f, with correlation %.2f."
  ),
  apim_distinguishable_diagram_estimates[["intercept_female"]],
  apim_distinguishable_diagram_estimates[["intercept_male"]],
  apim_distinguishable_diagram_estimates[["actor_female"]],
  apim_distinguishable_diagram_estimates[["actor_male"]],
  apim_distinguishable_diagram_estimates[["partner_female"]],
  apim_distinguishable_diagram_estimates[["partner_male"]],
  apim_distinguishable_diagram_residuals[["sd_female"]],
  apim_distinguishable_diagram_residuals[["sd_male"]],
  apim_distinguishable_diagram_residuals[["correlation"]]
)
```

The estimated coefficients map as follows: 

```{r fitted-distinguishable-apim-diagram, echo=FALSE, eval=has_glmmTMB, fig.width=9, fig.height=4.8, out.width="100%", fig.cap="Fitted cross-sectional distinguishable APIM for the example data. Fixed effects, residual standard deviations, and the residual correlation are extracted from the fitted model.", fig.alt=apim_distinguishable_fitted_alt}
draw_apim_diagram(
  "distinguishable",
  model = apim_distinguishable_model,
  labels = c(predictor = "Provided support", outcome = "Closeness")
)
```

## The exchangeable APIM

Conceptually, the exchangeable APIM constrains several of the effects to be equal: 

```{r exchangeable-apim-diagram, echo=FALSE, fig.width=9, fig.height=4.8, out.width="100%", fig.cap="Conceptual cross-sectional APIM for exchangeable dyads. The two members share one intercept, one actor effect, and one partner effect. Their outcome residuals have equal variances, yet still covary within dyads.", fig.alt="Path diagram for an exchangeable cross-sectional APIM. Both outcomes have the same intercept. Each member's predictor has the same actor effect on their own outcome and the same partner effect on the other member's outcome. The two outcome residuals have equal variances and covary."}
draw_apim_diagram("exchangeable")
```

Because the member labels are arbitrary, swapping members 1 and 2 does **not**
change the model. 

To estimate this model in a univariate MLM framework, we can draw a conceptual
diagram as such: 

```{r exchangeable-apim-member-diagram, echo=FALSE, fig.width=9, fig.height=5.2, out.width="100%", fig.cap="Individual-level representation of the exchangeable cross-sectional APIM used for the long-format multilevel model. Both members share the same intercept. Each member's own predictor has the shared actor effect, and the other member's predictor has the shared partner effect. The two residual variances are equal and the residuals may covary.", fig.alt="Two-panel path diagram for an exchangeable APIM. Both outcomes have the same intercept. For arbitrary member 1, X 1 is the actor predictor and X 2 is the partner predictor of Y 1. For arbitrary member 2, X 2 is the actor predictor and X 1 is the partner predictor of Y 2. Both panels use the same actor coefficient a and partner coefficient p, and their outcome residuals covary."}
draw_apim_member_diagram("exchangeable")
```


### Modeling the residual random-effects structure {#exchangeable-residual-structure}

For a simple random-intercept structure, equal member variances and their
covariance can be specified directly. The specification becomes more difficult
when the same exchangeability constraints must cover several random slopes. A
single homogeneous structure would impose one variance across intercepts and
slopes, whereas separate structures would omit their correlations.

The shared/difference representation works for residuals and other random-effect
terms, including random slopes. Following
@delrosarioPracticalGuideSpecifying2025, `dyadMLM::prepare_dyad_data()`
generates an arbitrary member-difference column, named `.dy_member_contrast_*`.
This contrast is `+1` for one member and `-1` for the other. The exchangeable
residual structure is represented by two separate random-effects terms: a
shared dyad random intercept and a random coefficient for this difference
column. Additional random slopes can be included in both blocks without
changing this logic.

We will now fit a simple exchangeable APIM and then use the `dyadMLM::recover_exchangeable_covariance()`
function that back-transforms the structure to the often more interpretable
member-level residual covariance matrix.


### Fitting the exchangeable APIM with glmmTMB

We use the same dataset as before, but do not distinguish males and females.
We can test distinguishability later by comparing this model with the prior
model.

We use `set_exchangeable_compositions` for the exchangeability constraints.
Another option would be to omit roles.

```{r prepare-cross-distinguishability}
apim_exchangeable_data <- dyadMLM::prepare_dyad_data(
  dyads_cross,
  dyad = coupleID,
  member = personID,
  role = gender,
  predictors = provided_support,
  keep_compositions = "female-male",
  set_exchangeable_compositions = "female-male",
  seed = 123
)

print(apim_exchangeable_data, n = 4)
```

We then use the columns to fit the model as follows: 

```{r fit-cross-distinguishability, eval = has_glmmTMB}
apim_exchangeable_model <- glmmTMB::glmmTMB(
  closeness ~
    
    # Pooled single intercept
    1 +
    
    # Pooled single actor and partner effects
    .dy_provided_support_actor +
    .dy_provided_support_partner +
    
    # Residual variance covariance matrix via the shared/difference
    # specification in two uncorrelated blocks
    us(1 | coupleID) +
    us(0 + .dy_member_contrast_female_x_male_arbitrary | coupleID),
  dispformula = ~ 0,
  family = gaussian(),
  data = apim_exchangeable_data
)

summary(apim_exchangeable_model)
```

We use the `dyadMLM::recover_exchangeable_covariance()` to recover the interpretable
variance-covariance matrix:

```{r backtransform-cross-exchangeable, eval = has_glmmTMB}
backtransformed <- dyadMLM::recover_exchangeable_covariance(apim_exchangeable_model)

# residual variance-covariance and SD-correlation representations
print(backtransformed)
```

The back-transformation follows directly from the shared and member-difference
random effects. If $u_j$ is the shared effect for dyad $j$ and
$\widetilde{u}_j$ its member-difference effect, the two member effects are

$$
u_{1j} = u_j + \widetilde{u}_j,
\qquad
u_{2j} = u_j - \widetilde{u}_j.
$$

Because the two fitted blocks are independent,

$$
\operatorname{Var}(u_{1j}) = \operatorname{Var}(u_{2j})
= \operatorname{Var}(u_j) + \operatorname{Var}(\widetilde{u}_j),
\qquad
\operatorname{Cov}(u_{1j}, u_{2j})
= \operatorname{Var}(u_j) - \operatorname{Var}(\widetilde{u}_j).
$$


```{r prepare-fitted-exchangeable-apim-diagram, include=FALSE, eval=has_glmmTMB}
apim_exchangeable_diagram_values <- .extract_apim_diagram_values(
  apim_exchangeable_model,
  type = "exchangeable"
)
apim_exchangeable_diagram_estimates <-
  apim_exchangeable_diagram_values$estimates
apim_exchangeable_diagram_residuals <-
  apim_exchangeable_diagram_values$residuals

apim_exchangeable_fitted_alt <- sprintf(
  paste(
    "Fitted exchangeable APIM. Intercept %.2f, actor effect %.2f, partner",
    "effect %.2f, common residual SD %.2f, and residual correlation %.2f."
  ),
  apim_exchangeable_diagram_estimates[["intercept"]],
  apim_exchangeable_diagram_estimates[["actor"]],
  apim_exchangeable_diagram_estimates[["partner"]],
  apim_exchangeable_diagram_residuals[["sd"]],
  apim_exchangeable_diagram_residuals[["correlation"]]
)
```



The output can now be mapped as follows:

```{r fitted-exchangeable-apim-diagram, echo=FALSE, eval=has_glmmTMB, fig.width=9, fig.height=4.8, out.width="100%", fig.cap="Fitted cross-sectional exchangeable APIM for the example data. The common member residual standard deviation and residual correlation are back-transformed from the fitted mean and difference components.", fig.alt=apim_exchangeable_fitted_alt}
draw_apim_diagram(
  "exchangeable",
  model = apim_exchangeable_model,
  labels = c(predictor = "Provided support", outcome = "Closeness")
)
```

## Testing distinguishability {#testing-distinguishability}

Distinguishability can be evaluated by comparing a full model in which the two
roles may differ with a restricted exchangeable model. This comparison tests
the imposed equality constraints jointly. Here, they concern the fixed
intercepts, actor effects, partner effects, and residual variances.

The two parameterizations require different generated columns, but both models
above use the same original observations.

[`dyadMLM::compare_nested_glmmTMB_models()`](https://pascal-kueng.github.io/dyadMLM/reference/compare_nested_glmmTMB_models.html)
verifies that both models use equivalent original
observations before performing the likelihood-ratio test:

```{r compare-cross-distinguishability, eval = has_glmmTMB}
dyadMLM::compare_nested_glmmTMB_models(
  apim_exchangeable_model,
  apim_distinguishable_model
)
```

The test provides evidence against all restrictions jointly, but it does not
show which parameter differs. The helper also cannot determine whether the
models are mathematically nested; that remains a modeling requirement. The
usual chi-squared reference distribution may be unreliable when a tested
variance parameter lies on its boundary.


## Intensive longitudinal APIMs

For longitudinal APIMs, time-varying predictors are decomposed into
within-person and between-person components before actor and partner variables
are constructed [@bolgerIntensiveLongitudinalMethods2013;
@gistelinckMultilevelAutoregressiveDyadic2020]. The default `"auto"` selects
`"2l"` when both `time` and `predictors` are supplied. The within-person (`cwp`)
component captures occasion-specific deviations from each member's observed
mean, whereas the between-person (`cbp`) component captures each member's
observed mean relative to the sample grand mean.

Note that observed person means used to construct the between-person (`cbp`)
predictors can be unreliable when each member contributes **few occasions**,
which can bias between-person estimates
[@gottfredsonStraightforwardApproachCoping2019].


### Concurrent ILD Gaussian APIM for distinguishable dyads

The decomposition above can be combined with role-specific effects. We first
retain the female-male distinction when preparing the data:

```{r prepare-ild-distinguishable-apim}
ild_distinguishable_data <- dyadMLM::prepare_dyad_data(
  dyads_ild,
  dyad = coupleID,
  member = personID,
  role = gender,
  time = diaryday,
  predictors = provided_support,
  model_types = "apim",
  keep_compositions = "female-male"
) |>
  dplyr::mutate(
    # we grand-mean center in this example to help convergence and
    # interpretation
    diaryday_gmc = diaryday - mean(diaryday)
  )
```

The following model allows the intercepts, time trends, and within- and
between-person actor and partner effects to differ between female and male
members. Centering `diaryday` makes the role-specific intercepts refer to the
average study day:

```{r fit-ild-distinguishable-apim, eval = has_glmmTMB}
ild_distinguishable_model <- glmmTMB::glmmTMB(
  closeness ~
    0 +

    # Role-specific intercepts
    .dy_is_female_x_male_female +
    .dy_is_female_x_male_male +

    # Role-specific time trends
    .dy_is_female_x_male_female:diaryday_gmc +
    .dy_is_female_x_male_male:diaryday_gmc +

    # Role-specific within-person actor effects
    .dy_is_female_x_male_female:.dy_provided_support_cwp_actor +
    .dy_is_female_x_male_male:.dy_provided_support_cwp_actor +

    # Role-specific within-person partner effects
    .dy_is_female_x_male_female:.dy_provided_support_cwp_partner +
    .dy_is_female_x_male_male:.dy_provided_support_cwp_partner +

    # Role-specific between-person actor effects
    .dy_is_female_x_male_female:.dy_provided_support_cbp_actor +
    .dy_is_female_x_male_male:.dy_provided_support_cbp_actor +

    # Role-specific between-person partner effects
    .dy_is_female_x_male_female:.dy_provided_support_cbp_partner +
    .dy_is_female_x_male_male:.dy_provided_support_cbp_partner +

    # Stable dyad-level covariance
    us(0 + .dy_is_female_x_male_female + .dy_is_female_x_male_male | coupleID) +

    # Same-occasion covariance
    us(0 + .dy_is_female_x_male_female + .dy_is_female_x_male_male |
         coupleID:diaryday),
  dispformula = ~ 0,
  family = gaussian(),
  data = ild_distinguishable_data
)
```

#### Random slopes

For example, role-specific within-person actor random slopes can be added by
replacing the stable dyad-level block above with:

```{r ild-distinguishable-random-slopes, eval = FALSE}
us(
  0 +
    .dy_is_female_x_male_female +
    .dy_is_female_x_male_male +
    .dy_is_female_x_male_female:.dy_provided_support_cwp_actor +
    .dy_is_female_x_male_male:.dy_provided_support_cwp_actor
  | coupleID
)
```

This block estimates the covariance among the two role-specific random
intercepts and actor slopes.

### Concurrent ILD Gaussian APIM for exchangeable dyads

In longitudinal Gaussian exchangeable APIMs, the sum-and-difference
parametrization from
@delrosarioPracticalGuideSpecifying2025 can be extended to the dyad-occasion
level to represent same-occasion residual dependence.

We first prepare within-person (`cwp`) and between-person (`cbp`) actor and
partner predictors. We retain the female-female dyads as one substantively
exchangeable composition:

```{r prepare-ild-exchangeable-apim}
ild_apim_data <- dyadMLM::prepare_dyad_data(
  dyads_ild,
  dyad = coupleID,
  member = personID,
  role = gender,
  time = diaryday,
  predictors = provided_support,
  model_types = "apim",
  keep_compositions = "female-female",
  seed = 123
)

print(ild_apim_data, n = 4)
```

The example below estimates same-day associations between support and closeness
and includes `diaryday` to adjust for a linear time trend. We also include an
actor random slope in both stable dyad-level blocks so that we can demonstrate
the random-slope back-transformation below.

```{r fit-ild-exchangeable-apim, eval = has_glmmTMB}
# Profile the fixed effects and use BFGS for stable cross-platform optimization.
ild_apim_control <- glmmTMB::glmmTMBControl(
  profile = TRUE,
  optimizer = stats::optim,
  optArgs = list(method = "BFGS")
)

ild_apim_model <- glmmTMB::glmmTMB(
  closeness ~
    1 +

    diaryday +

    # Within-person actor and partner effects
    .dy_provided_support_cwp_actor +
    .dy_provided_support_cwp_partner +

    # Between-person actor and partner effects
    .dy_provided_support_cbp_actor +
    .dy_provided_support_cbp_partner +

    # Stable exchangeable dyad-level covariance with actor random slopes
    us(1 + .dy_provided_support_cwp_actor | coupleID) + # shared intercept and slope
    us(0 + .dy_member_contrast_female_x_female_arbitrary + # difference intercept
         .dy_member_contrast_female_x_female_arbitrary:
           .dy_provided_support_cwp_actor                       # difference slope
       | coupleID) +

    # Same-occasion exchangeable covariance
    us(1 | coupleID:diaryday) +
    us(0 + .dy_member_contrast_female_x_female_arbitrary | coupleID:diaryday)

  , dispformula = ~ 0
  , family = gaussian()
  , data = ild_apim_data
  , control = ild_apim_control
)

summary(ild_apim_model)
```

We can now recover the member-level covariance matrices for both the stable
dyad effects and the same-occasion residual dependence. The two matched block
pairs are returned separately:

```{r backtransform-ild-exchangeable-apim, eval = has_glmmTMB}

recovered_covariance <- dyadMLM::recover_exchangeable_covariance(ild_apim_model)

print(recovered_covariance)

```

The `cwp` terms estimate actor and partner associations for occasion-specific
deviations from each member's usual support. The `cbp` terms estimate actor and
partner associations involving members' usual support levels. These are
concurrent associations; they do not by themselves represent temporal
carryover.

#### Extension to exchangeable random slopes {#exchangeable-random-slope-back-transformation}

The same shared/difference back-transformation described above applies
separately to every random slope. It also applies to covariances among the
random intercept, actor slope, and partner slope.

#### Testing random-effect constraints {#fitted-constraints-and-omitted-blocks}

The full model above estimates both a shared and a member-contrast actor random
slope. We can first test a smaller model that omits only the actor random slope
from the member-contrast block. The member-contrast random intercept and both
same-occasion blocks remain in the model:

```{r fit-ild-no-contrast-slope, eval = has_glmmTMB}
ild_apim_no_contrast_slope <- update(
  ild_apim_model,
  formula = . ~ . -
    us(0 +
         .dy_member_contrast_female_x_female_arbitrary +
         .dy_member_contrast_female_x_female_arbitrary:
           .dy_provided_support_cwp_actor
       | coupleID) +
    us(0 + .dy_member_contrast_female_x_female_arbitrary | coupleID)
)

dyadMLM::compare_nested_glmmTMB_models(
  ild_apim_no_contrast_slope,
  ild_apim_model
)
```

Without the member-contrast slope, the two members have identical actor random
slopes at the stable dyad level. We tell the back-transformation which fitted
blocks represent this constraint.

Since we omitted terms, automatic matching is no longer possible and we
need to tell `dyadMLM::recover_exchangeable_covariance()` what blocks
belong together.

```{r backtransform-ild-no-contrast-slope, eval = has_glmmTMB}
no_contrast_slope_covariance <- dyadMLM::recover_exchangeable_covariance(
  ild_apim_no_contrast_slope,
  block_pairings = list(
    dyad = list(
      shared_block =
        "us(1 + .dy_provided_support_cwp_actor | coupleID)",
      difference_block =
        "us(0 + .dy_member_contrast_female_x_female_arbitrary | coupleID)",
      difference_indicator =
        ".dy_member_contrast_female_x_female_arbitrary"
    )
  )
)

print(no_contrast_slope_covariance, representation = "sdcor")
```

We can impose the stronger constraint by omitting the full member-contrast
block at the stable dyad level. The same-occasion member-contrast block again
remains in the model:

```{r fit-ild-no-contrast-block, eval = has_glmmTMB}
ild_apim_no_contrast_block <- update(
  ild_apim_model,
  formula = . ~ . -
    us(0 +
         .dy_member_contrast_female_x_female_arbitrary +
         .dy_member_contrast_female_x_female_arbitrary:
           .dy_provided_support_cwp_actor
       | coupleID)
)

dyadMLM::compare_nested_glmmTMB_models(
  ild_apim_no_contrast_block,
  ild_apim_model
)
```

Here, both members have identical random intercepts and actor random slopes at
the stable dyad level. Because the full member-contrast block is absent, we
specify it as `NULL`:

```{r backtransform-ild-no-contrast-block, eval = has_glmmTMB}
no_contrast_block_covariance <- dyadMLM::recover_exchangeable_covariance(
  ild_apim_no_contrast_block,
  block_pairings = list(
    dyad = list(
      shared_block =
        "us(1 + .dy_provided_support_cwp_actor | coupleID)",
      difference_block = NULL,
      difference_indicator =
        ".dy_member_contrast_female_x_female_arbitrary"
    )
  )
)

print(no_contrast_block_covariance, representation = "sdcor")

```

These are constraints on the stable dyad-level random effects, not on the
same-occasion residual structure. Because variance constraints lie on the
boundary of the parameter space, the usual chi-squared reference distribution
for the likelihood-ratio tests should be interpreted cautiously.


### Current limitations of dyadic ILD designs in R

Concurrent dyadic ILD models can adjust for a time trend and account for stable
dyadic dependence and same-occasion partner dependence. They do not model
residual serial dependence, however, and therefore estimate concurrent
associations under the assumption that residuals from different days are
independent.

If the goal is to retain these concurrent APIM associations while accounting
for serial dependence, the closest extension is a dyadic residual dynamic
structural equation model (RDSEM). It keeps the concurrent regression separate
from a VAR model for its residuals [@asparouhovComparisonModelsAnalysis2020;
@mcneishPrimerTwoLevel2020]. This residual-VAR structure is not directly
available through the `glmmTMB` interface used here, and open-source support for
dyadic dynamic models remains very limited
[@delrosarioPracticalGuideSpecifying2025]. Within open-source R, such a model
generally requires custom TMB or Stan code.

### Dynamic ILD APIM example {#dynamic-models}

A **practical alternative** is a model with lagged outcomes, especially when
carryover or temporal dynamics are part of the research question
[@gistelinckMultilevelAutoregressiveDyadic2020]. In such a model, the
member's own lagged outcome represents stability and the partner's lagged
outcome represents influence. All other APIM predictor effects then describe
associations conditional on the members' prior outcomes.

By adding the outcome to `predictors` and selecting it with `lag1_predictors`,
`dyadMLM::prepare_dyad_data()` returns lag-1 raw and within-person scores alongside
the contemporaneous scores. Between-person scores are not lagged because they
describe stable differences between members.

For this example, we obtain the lagged actor and partner outcome columns through
the `lag1_predictors` argument:

```{r prepare-ild-apim-dynamic-model}
ild_apim_data_dynamic <- dyadMLM::prepare_dyad_data(
  dyads_ild,
  dyad = coupleID,
  member = personID,
  role = gender,
  time = diaryday,
  predictors = closeness,
  lag1_predictors = closeness,
  model_types = "apim",
  keep_compositions = "female-female",
  seed = 123
)
```

This returns all necessary variables for either choice, including lag-1 raw
and within-person closeness scores. Lags are matched at exactly
`diaryday - 1`, so omitted diary days are not bridged.

A simple fixed-slope dyadic stability and influence model
[@delrosarioPracticalGuideSpecifying2025]:

```{r fit-ild-stability-influence, eval = has_glmmTMB}
stability_influence <- glmmTMB::glmmTMB(
  closeness ~ 1 +

    # Stability (actor effect across time)
    .dy_closeness_actor_lag1 +

    # Influence (partner effect across time)
    .dy_closeness_partner_lag1 +

    # Linear time trend
    diaryday +

    # Stable exchangeable dyad-level covariance
    us(1 | coupleID) +
    us(0 + .dy_member_contrast_female_x_female_arbitrary | coupleID) +

    # Same-day exchangeable dyad-level covariance
    us(1 | coupleID:diaryday) +
    us(0 + .dy_member_contrast_female_x_female_arbitrary | coupleID:diaryday)

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

summary(stability_influence)
```

This model can be extended with contemporaneous actor and partner predictor
associations and other lagged predictors. Time-varying predictors should usually
be separated into within-person and between-person components. Their
contemporaneous coefficients are then conditional on both partners' prior
outcomes.

The model above uses raw lagged outcomes. The corresponding
within-person-centered lagged terms are `.dy_closeness_cwp_actor_lag1` and
`.dy_closeness_cwp_partner_lag1`.

Whether to use a raw or within-person-centered lagged outcome depends on the
research question and the data. Person-mean centering the outcome lag can bias
the average carryover estimate downward [@hamakerCenterNotCenter2015]. This
downward bias is known as Nickell bias [@nickellBiasesDynamicModels1981]. A raw
lag avoids this centering bias, but a standard random-intercept model can still
be biased because it assumes that the lag is unrelated to stable member levels
[@gistelinckMultilevelAutoregressiveSmall2021].

Both problems matter most when there are few occasions. In one simple
simulation, most approaches other than person-mean centering performed
acceptably from about ten occasions, but this is not a universal cutoff
[@gistelinckMultilevelAutoregressiveSmall2021]. For shorter panels, consider an
LD-APIM, which lets the first outcomes relate to the members' stable levels in a
wide-format SEM [@gistelinckMultilevelAutoregressiveDyadic2020].

These manifest-lag models are not equivalent to Mplus DSEM, which uses latent
person-mean centering by default and can estimate multivariate or residual
dynamics jointly [@mcneishPrimerTwoLevel2020]. For very short panels, however,
default DSEM may still be biased or unstable, so the LD-APIM recommendation
above may be preferable [@gistelinckMultilevelAutoregressiveSmall2021].

These lagged-outcome issues are separate from the earlier concern about
unreliable person means used to construct `cbp` predictors with few occasions.


---

From here, choose the model-specific vignette that matches the research
question:

- [Dyad-Individual Model vignette](dim.html) for the exchangeable DIM
  parameterization; or
- [Dyadic Score Model vignette](dsm.html) for the distinguishable DSM
  parameterization.

Or return to the
[online package overview](https://pascal-kueng.github.io/dyadMLM/).

A vignette with non-Gaussian generalized APIM examples is planned. 


## References

::: {#refs}
:::
