inst/doc/apim.R

In dyadMLM: Tools for Dyadic Multilevel Models

## ----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,
  dev.args = list(bg = "transparent")
)

## ----setup--------------------------------------------------------------------
library(dyadMLM)

## ----internal-setup, include=FALSE--------------------------------------------
has_glmmTMB <- requireNamespace("glmmTMB", quietly = TRUE)
apim_distinguishable_fitted_alt <-
  "Fitted distinguishable APIM diagram unavailable."
apim_exchangeable_fitted_alt <-
  "Fitted exchangeable APIM diagram unavailable."

## ----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, with both predictors centered using one pooled grand mean. 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. Grand-mean-centered female and male predictors each have an actor path to their own outcome and a partner path to the other member's outcome. Female and male outcomes have separate intercepts, and their residuals covary."----
draw_apim_diagram("distinguishable", predictors_centered = TRUE)

## ----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. Both predictors use one pooled grand-mean reference. 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 with grand-mean-centered predictors. 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 outcomes have separate intercepts and their residuals covary."----
draw_apim_member_diagram("distinguishable", predictors_centered = TRUE)

## ----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",
  add_apim_gmc_predictors = TRUE,
  # 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",
  include_arbitrary_member_contrast = TRUE,
  seed = 123
)

print(apim_distinguishable_data, n=4)

## ----fit-cross-distinguishable-gaussian, eval = has_glmmTMB-------------------
apim_distinguishable_model <- glmmTMB::glmmTMB(
  closeness ~

    # Gender-specific intercepts
    0 +
    .is_female +
    .is_male +

    # Gender-specific actor effects
    .is_female:.provided_support_gmc_actor +
    .is_male:.provided_support_gmc_actor +

    # Gender-specific partner effects
    .is_female:.provided_support_gmc_partner +
    .is_male:.provided_support_gmc_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 +
         .is_female +
         .is_male
       | coupleID)

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

summary(apim_distinguishable_model)

## ----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"]]
)

## ----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"),
  predictors_centered = TRUE
)

## ----exchangeable-apim-diagram, echo=FALSE, fig.width=9, fig.height=4.8, out.width="100%", fig.cap="Conceptual cross-sectional APIM for exchangeable dyads, with both predictors centered using one pooled grand mean. 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 with grand-mean-centered predictors. 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", predictors_centered = TRUE)

## ----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 predictors use one pooled grand-mean reference. The 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 with grand-mean-centered predictors. 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", predictors_centered = TRUE)

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

summary(apim_exchangeable_model)

## ----backtransform-cross-exchangeable, eval = has_glmmTMB---------------------
backtransformed <- dyadMLM::recover_exchangeable_covariance(apim_exchangeable_model)

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

## ----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"]]
)

## ----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"),
  predictors_centered = TRUE
)

## ----compare-cross-distinguishability, eval = has_glmmTMB---------------------
dyadMLM::compare_nested_models(
  apim_exchangeable_model,
  apim_distinguishable_model
)

## ----check-cross-apim, dev.args = list(bg = "white"), eval = has_glmmTMB, fig.width = 10, fig.height = 6.5, out.width = "100%", fig.alt = "Six predictive-check histograms. Member SDs and partner correlation are in the top row; dyad-average and half-difference summaries are in the bottom row. Red lines mark observed values."----
apim_simulations <- dyadMLM::simulate_dyad_responses(
  apim_exchangeable_model, seed = 123
)
dyadMLM::check_partner_dependence(
  apim_simulations,
  dyad = coupleID,
  role = gender,
  # Supply the fitting data because gender is not in the model formula.
  data = apim_distinguishable_data
)

## ----check-residuals-apim, dev.args = list(bg = "white"), eval = has_glmmTMB, dev = "png", fig.ext = "png", dpi = 96, fig.width = 10, fig.height = 8, out.width = "100%", fig.alt = c("Residual distribution and pattern for women and men: uniform QQ plots with blue envelopes in the top row and PIT quartiles across predicted outcomes in the bottom row. Red marks observed values.", "Residual spread for women and men: mean PIT distance histograms. Red marks observed values.", "Residual patterns for women and men: PIT quartiles across the partner's provided support with blue envelopes. Red curves show observed values.")----
dyadMLM::check_dyad_residuals(
  apim_simulations,
  dyad = coupleID,
  role = gender,
  predictors = ".provided_support_gmc_partner",
  data = apim_distinguishable_data
)

## ----check-outcomes-apim, dev.args = list(bg = "white"), eval = has_glmmTMB, dev = "png", fig.ext = "png", dpi = 96, fig.width = 10, fig.height = 12, out.width = "100%", fig.alt = "Outcome checks for women and men: outcome ECDFs in the top row, then histograms of the response SD and the largest absolute deviation. Red marks observed values."----
dyadMLM::check_dyad_outcomes(
  apim_simulations,
  dyad = coupleID,
  role = gender,
  data = apim_distinguishable_data
)

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

## ----fit-ild-distinguishable-apim, eval = has_glmmTMB-------------------------
ild_distinguishable_model <- glmmTMB::glmmTMB(
  closeness ~
    0 +

    # Role-specific intercepts
    .is_female +
    .is_male +

    # Role-specific time trends
    .is_female:diaryday_gmc +
    .is_male:diaryday_gmc +

    # Role-specific within-person actor effects
    .is_female:.provided_support_cwp_actor +
    .is_male:.provided_support_cwp_actor +

    # Role-specific within-person partner effects
    .is_female:.provided_support_cwp_partner +
    .is_male:.provided_support_cwp_partner +

    # Role-specific between-person actor effects
    .is_female:.provided_support_cbp_actor +
    .is_male:.provided_support_cbp_actor +

    # Role-specific between-person partner effects
    .is_female:.provided_support_cbp_partner +
    .is_male:.provided_support_cbp_partner +

    # Stable dyad-level covariance
    us(0 + .is_female + .is_male | coupleID) +

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

## ----fit-ild-distinguishable-ar1, eval = has_glmmTMB--------------------------
ild_distinguishable_ar_model <- glmmTMB::glmmTMB(
  closeness ~
    0 +

    # Role-specific intercepts
    .is_female +
    .is_male +

    # Role-specific time trends
    .is_female:diaryday_gmc +
    .is_male:diaryday_gmc +

    # Role-specific within-person actor effects
    .is_female:.provided_support_cwp_actor +
    .is_male:.provided_support_cwp_actor +

    # Role-specific within-person partner effects
    .is_female:.provided_support_cwp_partner +
    .is_male:.provided_support_cwp_partner +

    # Role-specific between-person actor effects
    .is_female:.provided_support_cbp_actor +
    .is_male:.provided_support_cbp_actor +

    # Role-specific between-person partner effects
    .is_female:.provided_support_cbp_partner +
    .is_male:.provided_support_cbp_partner +

    # Stable dyad-level covariance
    us(0 + .is_female + .is_male | coupleID) +

    # Same-occasion covariance
    us(0 + .is_female + .is_male |
         coupleID:diaryday) +

    # Role-specific residual persistence
    ar1(0 + .is_female:factor(diaryday) | coupleID) +
    ar1(0 + .is_male:factor(diaryday) | coupleID),
  dispformula = ~ 0,
  family = gaussian(),
  data = ild_distinguishable_data,
  # Non-default settings help this example converge.
  control = glmmTMB::glmmTMBControl(profile = TRUE)
)

glmmTMB::VarCorr(ild_distinguishable_ar_model)

## ----ild-distinguishable-random-slopes, eval = FALSE--------------------------
# us(
#   0 +
#     .is_female +
#     .is_male +
#     .is_female:.provided_support_cwp_actor +
#     .is_male:.provided_support_cwp_actor
#   | coupleID
# )

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

## ----fit-ild-exchangeable-apim, eval = has_glmmTMB----------------------------
ild_apim_model <- glmmTMB::glmmTMB(
  closeness ~
    1 +

    diaryday +

    # Within-person actor and partner effects
    .provided_support_cwp_actor +
    .provided_support_cwp_partner +

    # Between-person actor and partner effects
    .provided_support_cbp_actor +
    .provided_support_cbp_partner +

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

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

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

summary(ild_apim_model)

## ----fit-ild-exchangeable-ar1, eval = has_glmmTMB-----------------------------
ild_exchangeable_ar_model <- glmmTMB::glmmTMB(
  closeness ~
    1 +

    diaryday +

    # Within-person actor and partner effects
    .provided_support_cwp_actor +
    .provided_support_cwp_partner +

    # Between-person actor and partner effects
    .provided_support_cbp_actor +
    .provided_support_cbp_partner +

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

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

    # Pooled residual persistence across member series
    ar1(0 + factor(diaryday) | coupleID:personID),
  dispformula = ~ 0,
  family = gaussian(),
  data = ild_apim_data,
  # Non-default settings help this example converge.
  control = glmmTMB::glmmTMBControl(profile = TRUE)
)

glmmTMB::VarCorr(ild_exchangeable_ar_model)

## ----backtransform-ild-exchangeable-apim, eval = has_glmmTMB------------------

recovered_covariance <- dyadMLM::recover_exchangeable_covariance(
  ild_exchangeable_ar_model
)

print(recovered_covariance)


## ----fit-ild-exchangeable-random-slopes, eval = has_glmmTMB-------------------
ild_apim_random_slope_model <- update(
  ild_apim_model,
  formula = . ~ . -
    us(1 | coupleID) -
    us(0 + .member_contrast_arbitrary | coupleID) +
    us(1 + .provided_support_cwp_actor | coupleID) +
    us(0 +
         .member_contrast_arbitrary +
         .member_contrast_arbitrary:
           .provided_support_cwp_actor
       | coupleID),
  # Non-default settings help this example converge.
  control = glmmTMB::glmmTMBControl(
    profile = TRUE,
    optimizer = stats::optim,
    optArgs = list(method = "BFGS")
  )
)

random_slope_covariance <- dyadMLM::recover_exchangeable_covariance(
  ild_apim_random_slope_model
)

print(random_slope_covariance)

## ----fit-ild-no-contrast-slope, eval = has_glmmTMB----------------------------
ild_apim_no_contrast_slope <- update(
  ild_apim_random_slope_model,
  formula = . ~ . -
    us(0 +
         .member_contrast_arbitrary +
         .member_contrast_arbitrary:
           .provided_support_cwp_actor
       | coupleID) +
    us(0 + .member_contrast_arbitrary | coupleID)
)

dyadMLM::compare_nested_models(
  ild_apim_no_contrast_slope,
  ild_apim_random_slope_model
)

## ----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 + .provided_support_cwp_actor | coupleID)",
      difference_block =
        "us(0 + .member_contrast_arbitrary | coupleID)",
      difference_indicator =
        ".member_contrast_arbitrary"
    )
  )
)

print(no_contrast_slope_covariance, representation = "sdcor")

## ----fit-ild-no-contrast-block, eval = has_glmmTMB----------------------------
ild_apim_no_contrast_block <- update(
  ild_apim_random_slope_model,
  formula = . ~ . -
    us(0 +
         .member_contrast_arbitrary +
         .member_contrast_arbitrary:
           .provided_support_cwp_actor
       | coupleID)
)

dyadMLM::compare_nested_models(
  ild_apim_no_contrast_block,
  ild_apim_random_slope_model
)

## ----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 + .provided_support_cwp_actor | coupleID)",
      difference_block = NULL,
      difference_indicator =
        ".member_contrast_arbitrary"
    )
  )
)

print(no_contrast_block_covariance, representation = "sdcor")


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

## ----specify-ild-lagged-outcomes, eval = FALSE--------------------------------
# lagged_outcome_example <- glmmTMB::glmmTMB(
#   closeness ~ 1 +
# 
#     # Own-outcome lag (stability association)
#     .closeness_actor_lag1 +
# 
#     # Partner-outcome lag association
#     .closeness_partner_lag1 +
# 
#     # Linear time trend
#     diaryday +
# 
#     # Stable exchangeable dyad-level covariance
#     us(1 | coupleID) +
#     us(0 + .member_contrast_arbitrary | coupleID) +
# 
#     # Same-day exchangeable dyad-level covariance
#     us(1 | coupleID:diaryday) +
#     us(0 + .member_contrast_arbitrary | coupleID:diaryday)
# 
#   , dispformula = ~ 0
#   , family = gaussian()
#   , data = ild_apim_data_dynamic
#   # Non-default settings help this example converge.
#   , control = glmmTMB::glmmTMBControl(
#       profile = TRUE,
#       optimizer = stats::optim,
#       optArgs = list(method = "BFGS")
#     )
# )