inst/doc/dim.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())

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

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

## ----internal-setup, include=FALSE--------------------------------------------
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"),
  add_apim_gmc_predictors = TRUE,
  # 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
    .provided_support_dyad_mean_gmc +

    # Within-dyad effect
    .provided_support_within_dyad_dev +

    # Residual Gaussian covariance structure
    us(1 | coupleID) +
    us(0 + .member_contrast_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
    .provided_support_gmc_actor + .provided_support_gmc_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 + .member_contrast_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

b0_apim <- apim_coef[["(Intercept)"]]
b_actor <- apim_coef[[".provided_support_gmc_actor"]]
b_partner <- apim_coef[[".provided_support_gmc_partner"]]

b0_dim <- dim_coef[["(Intercept)"]]
b_mean <- dim_coef[[".provided_support_dyad_mean_gmc"]]
b_dev <- dim_coef[[".provided_support_within_dyad_dev"]]


cat("From APIM model:\n",
     "  intercept:                     ", round(b0_apim, 3), "\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",
     "  intercept:                     ", round(b0_dim, 3), "\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
)

## ----check-cross-dim, dev.args = list(bg = "white"), eval = has_glmmTMB, fig.width = 8, fig.height = 6, out.width = "100%", fig.alt = "Four predictive-check histograms for exchangeable dyads: common member SD, partner correlation, dyad-average SD, and half-difference RMS. Red lines mark observed values."----
dim_simulations <- dyadMLM::simulate_dyad_responses(dim_1, seed = 123)
dyadMLM::check_partner_dependence(
  dim_simulations,
  dyad = coupleID,
  role = gender,
  # Supply the fitting data because gender is not in the model formula.
  data = cross_exchangeable_data
)

## ----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
    .provided_support_cwp_dyad_mean +
    .provided_support_cwp_within_dyad_dev +

    # Between-person DIM
    .provided_support_cbp_dyad_mean +
    .provided_support_cbp_within_dyad_dev +

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

    # Residual (same-day) exchangeable dyad-level covariance
    us(1 | coupleID:diaryday) +
    us(0 + .member_contrast_arbitrary | coupleID:diaryday)

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

summary(dim_ILD)


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

    diaryday +

    # Within-person DIM
    .provided_support_cwp_dyad_mean +
    .provided_support_cwp_within_dyad_dev +

    # Between-person DIM
    .provided_support_cbp_dyad_mean +
    .provided_support_cbp_within_dyad_dev +

    # 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_exchangeable_data,
  # Non-default settings help this example converge.
  control = glmmTMB::glmmTMBControl(profile = TRUE)
)

glmmTMB::VarCorr(dim_ILD_ar_model)