## ----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)
dsm_fitted_alt <- "Fitted DSM diagram unavailable."

## ----prepare-cross-dsm--------------------------------------------------------
cross_dsm_data <- dyadMLM::prepare_dyad_data(
  dyads_cross,
  dyad = coupleID,
  member = personID,
  role = gender,
  predictors = provided_support,
  model_types = "dsm",
  # 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",
  dsm_role_order = c("female", "male")
)

print(cross_dsm_data, n = 4)

## ----conceptual-dsm-diagram, echo=FALSE, fig.width=9, fig.height=4.8, out.width="100%", fig.cap="Conceptual cross-sectional DSM. Predictor mean and predictor difference each predict both outcome scores.", fig.alt="Path diagram for a cross-sectional dyadic score model. The centered female-male predictor mean and female-minus-male predictor difference each predict the female-male outcome mean and female-minus-male outcome difference. Paths are labelled a11, a12, a21, and a22, and outcome intercepts are labelled a10 and a20."----
draw_dsm_diagram()

## ----conceptual-dsm-member-diagram, echo=FALSE, fig.width=9, fig.height=5.6, out.width="100%", fig.cap="Individual-level representation of the cross-sectional DSM used for the long-format multilevel model. The centered predictor mean and female-minus-male predictor difference appear on both member rows. For the female outcome, the intercept and slopes add one-half of the corresponding outcome-difference parameters. For the male outcome, they subtract one-half. The female and male residuals may have different variances and covary.", fig.alt="Two-panel path diagram for a female-minus-male dyadic score model. Both panels contain the centered predictor mean and female-minus-male predictor difference. For the female outcome, the intercept is a10 plus half a20, the predictor-mean coefficient is a11 plus half a21, and the predictor-difference coefficient is a12 plus half a22. For the male outcome, the same combinations use minus signs. The two member residuals covary."----
draw_dsm_member_diagram()

## ----fit-cross-dsm, eval=has_glmmTMB------------------------------------------
dsm_model <- glmmTMB::glmmTMB(
  closeness ~

    # Outcome-level intercept
    1 +

    # Predictor level -> outcome level (a11)
    .dy_provided_support_dyad_mean_gmc +

    # Predictor difference -> outcome level (a12)
    .dy_provided_support_within_dyad_diff +

    # Outcome-difference intercept (a20)
    .dy_dsm_role_contrast +

    # Predictor level -> outcome difference (a21)
    .dy_provided_support_dyad_mean_gmc:.dy_dsm_role_contrast +

    # Predictor difference -> outcome difference (a22)
    .dy_provided_support_within_dyad_diff:.dy_dsm_role_contrast +

    # Outcome-level and outcome-difference residual variances and their covariance
    us(1 + .dy_dsm_role_contrast | coupleID),
  dispformula = ~ 0,
  family = gaussian(),
  data = cross_dsm_data
)

summary(dsm_model)

## ----prepare-fitted-dsm-diagram, include=FALSE, eval=has_glmmTMB--------------
dsm_diagram_values <- .extract_dsm_diagram_values(dsm_model)
dsm_diagram_estimates <- dsm_diagram_values$estimates
dsm_diagram_residuals <- dsm_diagram_values$residuals

dsm_fitted_alt <- sprintf(
  paste(
    "Fitted DSM. Intercepts a10 %.2f and a20 %.2f; paths a11 %.2f, a12 %.2f,",
    "a21 %.2f, and a22 %.2f; residual SDs %.2f and %.2f, with correlation %.2f."
  ),
  dsm_diagram_estimates[["a10"]],
  dsm_diagram_estimates[["a20"]],
  dsm_diagram_estimates[["a11"]],
  dsm_diagram_estimates[["a12"]],
  dsm_diagram_estimates[["a21"]],
  dsm_diagram_estimates[["a22"]],
  dsm_diagram_residuals[["sd_mean"]],
  dsm_diagram_residuals[["sd_difference"]],
  dsm_diagram_residuals[["correlation"]]
)

## ----fitted-dsm-diagram, echo=FALSE, eval=has_glmmTMB, fig.width=9, fig.height=4.8, out.width="100%", fig.cap="Fitted cross-sectional DSM for the example data. The nodes identify the mean and difference scores; edge and intercept labels show the estimated DSM coefficients, and the residual labels show the estimated score-component standard deviations and correlation.", fig.alt=dsm_fitted_alt----
draw_dsm_diagram(
  model = dsm_model,
  labels = c(predictor = "provided support", outcome = "closeness")
)

## ----prepare-cross-dsm-inverted-----------------------------------------------
cross_dsm_data_inverted <- dyadMLM::prepare_dyad_data(
  dyads_cross,
  dyad = coupleID,
  member = personID,
  role = gender,
  predictors = provided_support,
  # Request APIM columns too for comparison below.
  model_types = c("dsm", "apim"),
  keep_compositions = "female-male",
  dsm_role_order = c("male", "female")
)

## ----fit-cross-dsm-inverted, eval=has_glmmTMB---------------------------------
dsm_model_inverted <- glmmTMB::glmmTMB(
  closeness ~
    .dy_provided_support_dyad_mean_gmc +
    .dy_provided_support_within_dyad_diff +
    .dy_dsm_role_contrast +
    .dy_provided_support_dyad_mean_gmc:.dy_dsm_role_contrast +
    .dy_provided_support_within_dyad_diff:.dy_dsm_role_contrast +
    us(1 + .dy_dsm_role_contrast | coupleID),
  dispformula = ~ 0,
  family = gaussian(),
  data = cross_dsm_data_inverted
)

female_minus_male <- glmmTMB::fixef(dsm_model)$cond
male_minus_female <- glmmTMB::fixef(dsm_model_inverted)$cond

knitr::kable(
  data.frame(
    `model term` = names(female_minus_male),
    `female - male` = unname(female_minus_male),
    `male - female` = unname(male_minus_female),
    check.names = FALSE
  ),
  digits = 3,
  align = c("l", "r", "r")
)

## ----fit-cross-apim, eval=has_glmmTMB-----------------------------------------
apim_model <- glmmTMB::glmmTMB(
  closeness ~
    # Role-specific intercepts
    0 +
    .dy_is_female_x_male_female +
    .dy_is_female_x_male_male +

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

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

    # Role-specific Gaussian residual covariance structure
    us(0 +
         .dy_is_female_x_male_female +
         .dy_is_female_x_male_male
       | coupleID),
  dispformula = ~ 0,
  family = gaussian(),
  data = cross_dsm_data_inverted
)

## ----compare-cross-fit, eval=has_glmmTMB--------------------------------------
data.frame(
  model = c("DSM: female - male", "DSM: male - female", "APIM"),
  AIC = round(c(AIC(dsm_model), AIC(dsm_model_inverted), AIC(apim_model)), 3),
  BIC = round(c(BIC(dsm_model), BIC(dsm_model_inverted), BIC(apim_model)), 3),
  logLik = round(c(
    as.numeric(logLik(dsm_model)),
    as.numeric(logLik(dsm_model_inverted)),
    as.numeric(logLik(apim_model))
  ), 3)
)

## ----compare-cross-coefficients, eval=has_glmmTMB, echo=FALSE-----------------
apim_coef <- glmmTMB::fixef(apim_model)$cond
dsm_coef <- glmmTMB::fixef(dsm_model)$cond

b0_female <- apim_coef[[".dy_is_female_x_male_female"]]
b0_male <- apim_coef[[".dy_is_female_x_male_male"]]
b_actor_female <- apim_coef[[
  ".dy_is_female_x_male_female:.dy_provided_support_actor"
]]
b_actor_male <- apim_coef[[
  ".dy_is_female_x_male_male:.dy_provided_support_actor"
]]
b_partner_female <- apim_coef[[
  ".dy_is_female_x_male_female:.dy_provided_support_partner"
]]
b_partner_male <- apim_coef[[
  ".dy_is_female_x_male_male:.dy_provided_support_partner"
]]

mu_X <- mean(
  (cross_dsm_data_inverted$.dy_provided_support_actor +
     cross_dsm_data_inverted$.dy_provided_support_partner) / 2 -
    cross_dsm_data_inverted$.dy_provided_support_dyad_mean_gmc,
  na.rm = TRUE
)

b_sum_female <- b_actor_female + b_partner_female
b_sum_male <- b_actor_male + b_partner_male
b_difference_female <- b_actor_female - b_partner_female
b_difference_male <- b_actor_male - b_partner_male

a11_from_apim <- (b_sum_female + b_sum_male) / 2
a21_from_apim <- b_sum_female - b_sum_male
a12_from_apim <- (b_difference_female - b_difference_male) / 4
a22_from_apim <- (b_difference_female + b_difference_male) / 2

apim_to_dsm <- c(
  a10 = (b0_female + b0_male) / 2 + mu_X * a11_from_apim,
  a11 = a11_from_apim,
  a12 = a12_from_apim,
  a20 = b0_female - b0_male + mu_X * a21_from_apim,
  a21 = a21_from_apim,
  a22 = a22_from_apim
)

a10_dsm <- dsm_coef[["(Intercept)"]]
a11_dsm <- dsm_coef[[".dy_provided_support_dyad_mean_gmc"]]
a12_dsm <- dsm_coef[[".dy_provided_support_within_dyad_diff"]]
a20_dsm <- dsm_coef[[".dy_dsm_role_contrast"]]
a21_dsm <- dsm_coef[[
  ".dy_provided_support_dyad_mean_gmc:.dy_dsm_role_contrast"
]]
a22_dsm <- dsm_coef[[
  ".dy_provided_support_within_dyad_diff:.dy_dsm_role_contrast"
]]

knitr::kable(
  data.frame(
    `DSM path` = names(apim_to_dsm),
    `From APIM transformation` = round(unname(apim_to_dsm), 3),
    `From DSM model` = round(c(
      a10_dsm,
      a11_dsm,
      a12_dsm,
      a20_dsm,
      a21_dsm,
      a22_dsm
    ), 3),
    check.names = FALSE
  ),
  align = c("l", "r", "r"),
  caption = paste0(
    "APIM-to-DSM fixed-effect transformation (centering constant = ",
    round(mu_X, 3),
    ")."
  )
)

## ----compare-cross-random-effects, eval=has_glmmTMB---------------------------
apim_vcov <- as.matrix(glmmTMB::VarCorr(apim_model)$cond$coupleID)
dsm_vcov <- as.matrix(glmmTMB::VarCorr(dsm_model)$cond$coupleID)

rotation <- rbind(
  outcome_level = c(0.5, 0.5),
  outcome_difference = c(1, -1)
)
apim_to_dsm_vcov <- rotation %*% apim_vcov %*% t(rotation)

data.frame(
  parameter = c(
    "Var(outcome mean)",
    "Cov(outcome mean, outcome diff)",
    "Var(outcome diff)"
  ),
  from_DSM = round(c(
    dsm_vcov[1, 1],
    dsm_vcov[1, 2],
    dsm_vcov[2, 2]
  ), 3),
  from_APIM_transformation = round(c(
    apim_to_dsm_vcov[1, 1],
    apim_to_dsm_vcov[1, 2],
    apim_to_dsm_vcov[2, 2]
  ), 3)
)

## ----prepare-ild-dsm----------------------------------------------------------
ild_dsm_data <- dyadMLM::prepare_dyad_data(
  dyads_ild,
  dyad = coupleID,
  member = personID,
  role = gender,
  time = diaryday,
  predictors = provided_support,
  model_types = "dsm",
  keep_compositions = "female-male",
  dsm_role_order = c("female", "male")
)

print(ild_dsm_data, n = 4)

## ----fit-ild-dsm, eval = has_glmmTMB------------------------------------------

dsm_ILD <- glmmTMB::glmmTMB(
  closeness ~

    # Outcome-level intercept and linear time trend
    1 +
    diaryday +

    # Within-person predictor level -> outcome level
    .dy_provided_support_cwp_dyad_mean +

    # Within-person predictor difference -> outcome level
    .dy_provided_support_cwp_within_dyad_diff +

    # Between-person predictor level -> outcome level
    .dy_provided_support_cbp_dyad_mean +

    # Between-person predictor difference -> outcome level
    .dy_provided_support_cbp_within_dyad_diff +

    # Outcome-difference intercept and linear time trend
    .dy_dsm_role_contrast +
    diaryday:.dy_dsm_role_contrast +

    # Within-person predictor level and difference -> outcome difference
    .dy_provided_support_cwp_dyad_mean:.dy_dsm_role_contrast +
    .dy_provided_support_cwp_within_dyad_diff:.dy_dsm_role_contrast +

    # Between-person predictor level and difference -> outcome difference
    .dy_provided_support_cbp_dyad_mean:.dy_dsm_role_contrast +
    .dy_provided_support_cbp_within_dyad_diff:.dy_dsm_role_contrast +

    # Stable outcome-level and outcome-difference covariance
    us(1 + .dy_dsm_role_contrast | coupleID) +

    # Same-day outcome-level and outcome-difference covariance
    us(1 + .dy_dsm_role_contrast | coupleID:diaryday),
  dispformula = ~ 0,
  family = gaussian(),
  data = ild_dsm_data
)

summary(dsm_ILD)