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

dim_grid_helper <- if (file.exists("dim-grid.Rinc")) {
  "dim-grid.Rinc"
} else {
  file.path("vignettes", "dim-grid.Rinc")
}
sys.source(dim_grid_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)
has_htmltools <- requireNamespace("htmltools", quietly = TRUE)
dim_fitted_alt <- "Fitted DIM diagram unavailable."

## ----prepare-cross-exchangeable-----------------------------------------------
cross_exchangeable_data <- dyadMLM::prepare_dyad_data(
  dyads_cross,
  dyad = coupleID,
  member = personID,
  role = gender,
  predictors = provided_support,
  # Create both APIM and DIM columns for comparison.
  model_types = c("apim", "dim"),
  # All three observed compositions in `dyads_cross` are detected and retained by
  # default. This example focuses on `female-female` dyads, so we restrict the
  # analysis here.
  keep_compositions = "female-female",
  seed = 123
)

# Print the first two dyads.
print(cross_exchangeable_data, n = 4)

## ----conceptual-dim-diagram, echo=FALSE, fig.width=9, fig.height=4.8, out.width="100%", fig.cap="Cross-sectional DIM. The dyad mean has a between-dyad effect, and the within-dyad member deviation has a within-dyad effect. Mean and deviation residuals are uncorrelated under exchangeability. Both members' residuals and their correlation can be obtained from these residuals.", fig.alt="Path diagram for a cross-sectional Dyad-Individual Model. The grand-mean-centered dyad mean predicts the outcome mean through the between-dyad effect b mean. Member i's within-dyad member deviation predicts the same member's outcome deviation through the within-dyad effect b dev. There are no cross-paths, and only the outcome-mean equation has intercept b zero."----
draw_dim_diagram()

## ----conceptual-dim-member-diagram, echo=FALSE, fig.width=9, fig.height=5.2, out.width="100%", fig.cap="Individual-level representation of the cross-sectional DIM used for the long-format multilevel model.", fig.alt="Path diagram for two arbitrarily labelled members of an exchangeable dyad. For each member, the grand-mean-centered dyad mean and the member's own within-dyad member deviation predict the individual outcome. Both members share the same between-dyad and within-dyad coefficients and residual standard deviation, while their outcome residuals are correlated."----
draw_dim_member_diagram()

## ----fit-cross-dim, eval=has_glmmTMB------------------------------------------

dim_1 <- glmmTMB::glmmTMB(
  closeness ~

    # Pooled fixed intercept
    1 +

    # Between-dyad effect
    .dy_provided_support_dyad_mean_gmc +

    # Within-dyad effect
    .dy_provided_support_within_dyad_dev +

    # Residual Gaussian covariance structure
    us(1 | coupleID) +
    us(0 + .dy_member_contrast_female_x_female_arbitrary | coupleID)
  , dispformula = ~ 0
  , family = gaussian()
  , data = cross_exchangeable_data
)

summary(dim_1)

## ----prepare-fitted-dim-diagram, include=FALSE, eval=has_glmmTMB--------------
dim_diagram_values <- .extract_dim_diagram_values(dim_1)
dim_diagram_estimates <- dim_diagram_values$estimates
dim_diagram_residuals <- dim_diagram_values$residuals

dim_fitted_alt <- sprintf(
  paste(
    "Fitted DIM. Intercept %.2f, between-dyad effect %.2f, within-dyad",
    "effect %.2f, and mean/deviation residual SDs %.2f and %.2f; their",
    "correlation is fixed at zero."
  ),
  dim_diagram_estimates[["b0"]],
  dim_diagram_estimates[["b_mean"]],
  dim_diagram_estimates[["b_dev"]],
  dim_diagram_residuals[["sd_mean"]],
  dim_diagram_residuals[["sd_difference"]]
)

## ----fitted-dim-diagram, echo=FALSE, eval=has_glmmTMB, fig.width=9, fig.height=4.8, out.width="100%", fig.cap="Estimated fixed effects and residual-component standard deviations from the cross-sectional Gaussian DIM in its mean-and-deviation representation. The intercept belongs to the outcome-mean equation.", fig.alt=dim_fitted_alt----
draw_dim_diagram(
  model = dim_1,
  effect_unit = "dyad",
  labels = c(predictor = "provided support", outcome = "closeness")
)

## ----fit-cross-apim, eval = has_glmmTMB---------------------------------------

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

    # Fixed effects APIM
    .dy_provided_support_actor + .dy_provided_support_partner +

    # Since both models are equivalent, the same random-effects structure
    # can be used. See the APIM vignette to learn how to back-transform
    # these blocks to a full actor-partner covariance matrix.
    us(1 | coupleID) +
    us(0 + .dy_member_contrast_female_x_female_arbitrary | coupleID)
  , dispformula = ~ 0
  , family = gaussian()
  , data = cross_exchangeable_data
)


## ----compare-cross-fit, eval=has_glmmTMB--------------------------------------
data.frame(
  model = c("DIM", "APIM"),
  AIC = c(AIC(dim_1), AIC(apim_1)),
  BIC = c(BIC(dim_1), BIC(apim_1)),
  logLik = c(as.numeric(logLik(dim_1)), as.numeric(logLik(apim_1)))
)

## ----compare-cross-coefficients, echo = TRUE, eval = has_glmmTMB--------------
apim_coef <- glmmTMB::fixef(apim_1)$cond
dim_coef <- glmmTMB::fixef(dim_1)$cond

b_actor <- apim_coef[[".dy_provided_support_actor"]]
b_partner <- apim_coef[[".dy_provided_support_partner"]]

b_mean <- dim_coef[[".dy_provided_support_dyad_mean_gmc"]]
b_dev <- dim_coef[[".dy_provided_support_within_dyad_dev"]]


cat("From APIM model:\n",
     "  actor effect:                  ", round(b_actor, 3), "\n",
     "  partner effect:                ", round(b_partner, 3), "\n\n",

     "DIM transformation:\n",
     "  b_mean = b_actor + b_partner:  ", round(b_actor + b_partner, 3), "\n",
     "  b_dev = b_actor - b_partner:   ", round(b_actor - b_partner, 3), "\n\n",

     "From DIM model:\n",
     "  between-dyad effect:           ", round(b_mean, 3), "\n",
     "  within-dyad effect:            ", round(b_dev, 3), "\n"
)


## ----interactive-dim-grid, echo=FALSE, eval=has_htmltools---------------------
dim_reparameterization_grid(
  if (has_glmmTMB && exists("apim_1")) apim_1 else NULL,
  predictor_label = "Provided support",
  outcome_label = "the linear predictor",
  limit = 5
)

## ----prepare-ild-exchangeable-------------------------------------------------
ild_exchangeable_data <- dyadMLM::prepare_dyad_data(
  dyads_ild,
  dyad = coupleID,
  member = personID,
  role = gender,
  time = diaryday,
  predictors = provided_support,
  model_types = c("apim", "dim"),
  keep_compositions = "female-female",
  seed = 123
)

print(ild_exchangeable_data)

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

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

    diaryday +

    # Within-person DIM
    .dy_provided_support_cwp_dyad_mean +
    .dy_provided_support_cwp_within_dyad_dev +

    # Between-person DIM
    .dy_provided_support_cbp_dyad_mean +
    .dy_provided_support_cbp_within_dyad_dev +

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

    # Residual (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_exchangeable_data
)

summary(dim_ILD)