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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

recovered_covariance <- dyadMLM::recover_exchangeable_covariance(ild_apim_model)

print(recovered_covariance)


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

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

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

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


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

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