---
title: "Dyadic Data Analysis with R"
subtitle: "Intensive Longitudinal Workshop Exercises"
author:
  - "Pascal Küng"
  - "Turu Stadler"
format:
  html:
    embed-resources: true
    toc: true
    toc-title: "ILD Exercises"
editor: source
---

# Exercises

These exercises extend the cross-sectional APIM workflow.

Download the complete workshop materials at any time from the
[workshop page](https://pascal-kueng.github.io/dyadMLM/workshop/). Package
documentation and vignettes are available on the
[`dyadMLM` homepage](https://pascal-kueng.github.io/dyadMLM/).


# Getting started

Keep `04_exercises_ild.Rmd`, `dyadic-ild-exercise.rds`, `00_setup.R`,
`hints.R`, and `helpers.R` in the same folder.

If you have not already run the workshop setup, run the next chunk once, in a
fresh R session. It installs or updates the workshop packages and any
dependencies that need updating, but does not update unrelated packages. This
may take several minutes.

```{r eval=FALSE}
source("00_setup.R")
```

Load the required packages, optional hints, and ILD diagnostic helpers:

```{r message=FALSE, warning=FALSE}
library(glmmTMB)
library(dyadMLM)
library(wbCorr)
library(easystats, quietly = TRUE)
library(DHARMa)
library(dplyr)
library(ggplot2)

source("hints.R")
source("helpers.R")
```

Optional hints appear as commented `hint()` calls. Uncomment and run them one
at a time.


# The ILD workflow

The same six-step workflow still applies:

1. Define dyads, roles, and research questions
2. Prepare the data and validate the dyadic structure
3. Describe and visualize measures and observed relationships
4. Estimate the models
5. Verify model adequacy, role differences, and robustness
6. Visualize, interpret, and report results

ILD additionally requires attention to matched occasions, time gaps,
within- and between-person variation, time trends, and residual temporal
dependence.


# 1. Define dyads, roles, and research questions

The cross-sectional exercise let you choose the predictor and outcome. Here,
everyone uses one prescribed question and a focused simulated dataset so that
the focus remains on the new ILD decisions and results can be compared across
groups. This is a teaching choice, not a claim that there is only one sensible
question or causal direction.

Load and inspect the focused workshop dataset:

```{r}
raw_ild_data <- readRDS("dyadic-ild-exercise.rds")

raw_ild_data
```

The data contain **simulated daily scores for female--male dyads**.
Gender is the distinguishing variable. Higher scores indicate more provided
support and greater closeness.

Our concurrent research question is:

> **How are daily deviations and stable between-person differences in provided
> support associated with one's own and one's partner's closeness on the same
> day?**

Unlike the prospective workshop example, this question concerns the same day.
Therefore, no lagged predictors are needed.

Translate this into four questions or hypotheses:

> [Insert the within-person actor question]

> [Insert the within-person partner question]

> [Insert the between-person actor question]

> [Insert the between-person partner question]

Briefly state what the concurrent design can and cannot establish:

> [Insert your answer]


# 2. Prepare the data and validate the dyadic structure

Prepare a **distinguishable** APIM dataset with
`dyadMLM::prepare_dyad_data()`.

You need actor and partner variables for both the within-person and
between-person components. For help with most functions, run `?function_name`.

Save the result as `my_ild_data`.

```{r eval=FALSE}
# hint("ild_preparation", "solution")
```

```{r eval=FALSE}
my_ild_data <- dyadMLM::prepare_dyad_data(
  data = raw_ild_data,
  # Complete the required arguments.
)
```

Add two time variables:

- `diaryday_c`: centered at the study midpoint and scaled in full-study units
- `diaryday_f`: a factor containing every scheduled day from 0 through 13

These names are used below.

```{r eval=FALSE}
# hint("ild_preparation", "time_variables")
```

```{r eval=FALSE}
my_ild_data <- my_ild_data |>
  dplyr::mutate(
    # Create diaryday_c and diaryday_f.
  )
```

Inspect the prepared data and its generated columns:

```{r eval=FALSE}
my_ild_data
```

```{r eval=FALSE}
names(my_ild_data)
```


## Validate the dyadic and temporal structure

This teaching dataset is deliberately complete and balanced. Use separate
summaries for women and men to check that both roles have the expected number
of observations and day range, and inspect variable ranges and missingness.

```{r eval=FALSE}
my_ild_data |>
  dplyr::filter(gender == "male") |>
  dplyr::select(diaryday, provided_support, closeness) |>
  report::report_table()

my_ild_data |>
  dplyr::filter(gender == "female") |>
  dplyr::select(diaryday, provided_support, closeness) |>
  report::report_table()
```

Summarize what you found:

> [Insert your data-structure and missingness summary]


# 3. Describe and visualize measures and observed relationships


## Describe the focal variables

Reuse the simple tools from the cross-sectional exercise: summarize
`provided_support` and `closeness` overall and by role, then inspect
role-specific histograms. Means and standard deviations are sufficient here.

```{r}
# hint("ild_descriptives", "summary")
```

```{r}
# Add your overall and role-specific summaries here.
```

```{r}
# hint("ild_descriptives", "histograms")
```

```{r}
# Add your role-specific histograms here.
```


## Plot paired trajectories

Plot both variables over time for a manageable selection of dyads. Make the two
members identifiable and avoid drawing every trajectory in one panel.

```{r eval=FALSE}
trajectory_dyads <- sort(unique(my_ild_data$coupleID))[1:4]

trajectory_data <- my_ild_data |>
  dplyr::filter(coupleID %in% trajectory_dyads) |>
  dplyr::select(
    coupleID, personID, diaryday, gender,
    provided_support, closeness
  ) |>
  tidyr::pivot_longer(
    cols = c(provided_support, closeness),
    names_to = "measure",
    values_to = "score"
  )
```

```{r}
# hint("ild_trajectories", "structure")
# hint("ild_trajectories", "solution")
```

```{r}
# Complete the mappings and facet specification, then run the plot:
# ggplot2::ggplot(
#   trajectory_data,
#   ggplot2::aes(
#     x = ...,
#     y = ...,
#     color = ...,
#     group = ...
#   )
# ) +
#   ggplot2::geom_line(na.rm = TRUE) +
#   ggplot2::geom_point(na.rm = TRUE) +
#   ggplot2::facet_grid(... ~ ..., scales = "free_y") +
#   ggplot2::labs(x = "Diary day", y = NULL, color = "Gender") +
#   ggplot2::theme_minimal()
```


## Optional extension: separate within- and between-person information

If time permits, adapt the wide-data correlation workflow from the
cross-sectional exercise. Create one row per couple and day, with separate
female and male columns for each focal variable. Then use `wbCorr()` with
`coupleID` as the cluster to inspect:

- ICCs for each role-specific series
- correlations among daily deviations within couples
- correlations among role-specific person means between couples

Map the relevant entries to the actor and partner questions. Treat these
correlations as descriptive rather than as replacements for the APIM.
Multilevel reliability would also matter for multi-item measures, but the
item-level data needed to estimate it are not included here.

```{r}
# hint("ild_level_associations", "structure")
# hint("ild_level_associations", "solution")
```

```{r}
# Add your optional level-specific association code here.
```

Describe the main patterns:

> [Insert your descriptive interpretation]


# 4. Estimate the role-specific AR(1) model

Fit a Gaussian distinguishable APIM with:

- role-specific intercepts and linear time trends
- role-specific within-person actor and partner effects
- role-specific between-person actor and partner effects
- stable role-specific covariance at the dyad level
- separate female and male AR(1) residual processes
- same-occasion role-specific covariance at the dyad--day level

Use `dispformula = ~ 0` so the modeled dyad--day block represents the Gaussian
observation-level residual structure. Adapt the distinguishable ILD model from
the applied tutorial and inspect `names(my_ild_data)` when selecting columns.

The applied tutorial shows how residual temporal dependence motivates AR(1).
For this focused extension, fit the planned role-specific AR(1) specification
directly. This avoids spending workshop time estimating a deliberately
misspecified no-AR intermediate model.

To keep this exercise focused and fitting times manageable, use random
intercepts but no random slopes. In substantive analyses, consider
theoretically plausible within-person random slopes when the data can support
them; their omission here is a teaching and computational simplification.

The covariance model has several components. Use this tested optimizer control
for every candidate so that the comparison does not also change the estimation
strategy:

```{r}
ild_control <- glmmTMB::glmmTMBControl(
  profile = TRUE,
  optimizer = stats::optim,
  optArgs = list(method = "BFGS"),
  optCtrl = list(maxit = 2000, reltol = 1e-10)
)
```

```{r eval=FALSE}
# hint("ild_ar1", "role_specific_syntax")
# hint("ild_models", "distinguishable_ar1")
```

```{r eval=FALSE}
ild_ar1_start <- make_ild_ar1_start(my_ild_data)

distinguishable_ar1_model <- glmmTMB::glmmTMB(
  formula = YOUR_ROLE_SPECIFIC_AR1_FORMULA,
  dispformula = ~ 0,
  family = gaussian(),
  data = my_ild_data,
  control = ild_control,
  start = ild_ar1_start
)

summary(distinguishable_ar1_model)
```

The start helper fits only an ordinary fixed-effects regression, then supplies
generic interior values for the covariance parameters. It does not fit a
mixed model or use the simulation truth.

Inspect the model output for convergence or Hessian warnings, boundary
estimates, and the plausibility of all estimates before continuing. A zero
convergence code is not enough when the Hessian is non-positive-definite or a
covariance estimate is on its boundary.


# 5. Verify model adequacy, role differences, and robustness


## Does the model reproduce the temporal dependence?

Use the workshop helper to simulate block-whitened DHARMa residuals. The helper
accounts for fitted dependence within each dyad rather than treating every
person-day as independent. Pass `role` and request `by_role = TRUE` so the
female and male lag-1 checks are not averaged together.

```{r eval=FALSE}
distinguishable_dharma_residuals <- simulate_ild_dharma(
  model = distinguishable_ar1_model,
  dyad = my_ild_data$coupleID,
  member = my_ild_data$personID,
  time = my_ild_data$diaryday,
  role = my_ild_data$gender,
  n = 250,
  seed = 123
)

plot(
  distinguishable_dharma_residuals,
  title = "Block-whitened DHARMa residuals"
)

test_ild_lag1(distinguishable_dharma_residuals, by_role = TRUE)
```

Consider:

- overall calibration and visible residual patterns
- the magnitude of the observed lag-1 correlation for each role
- whether each correlation is typical under simulations from the fitted model

The plot uses block-whitened scaled residuals. The lag-1 helper instead uses a
raw response-residual correlation and compares it with the same statistic in
simulations from the fitted model.

The exercise uses 250 simulations to keep this extension manageable during the
workshop. Increase that number for a final substantive analysis.

Do not base the decision on a significance label alone.

> [Insert your diagnostic interpretation]

```{r eval=FALSE}
# hint("ild_diagnostics", "solution")
```

Use one stable name for the model comparison:

```{r eval=FALSE}
distinguishable_model <- distinguishable_ar1_model
```


## Do the focal associations differ by role?

The cross-sectional exercise covered a global exchangeability test. Here, use a
more focused comparison: pool only the four within- and between-person
actor/partner associations. Keep role-specific intercepts, time trends, stable
covariance, same-occasion covariance, and both AR(1) processes unchanged.

This creates a regular four-degree-of-freedom test of one clear question:

> **Do the four focal APIM associations differ between women and men?**

Because the temporal and covariance terms are identical in both candidates,
the comparison does not mix that question with a change in residual structure.

```{r eval=FALSE}
# hint("ild_models", "pooled_apim_effects")
```

```{r eval=FALSE}
pooled_apim_effects_start <- make_pooled_apim_effects_start(
  distinguishable_model
)

pooled_apim_effects_model <- update(
  distinguishable_model,
  formula = YOUR_POOLED_EFFECTS_FORMULA_UPDATE,
  start = pooled_apim_effects_start
)

summary(pooled_apim_effects_model)
```

Check estimation, then compare the matched models:

```{r eval=FALSE}
dyadMLM::compare_nested_models(
  distinguishable_model,
  pooled_apim_effects_model
)
```

Choose the model you will retain. Consider statistical evidence, parsimony, and
whether common actor and partner associations are substantively defensible.
Do not call the restricted model exchangeable: its baseline, covariance, and
AR(1) parameters remain role-specific.

> [Insert your role-difference conclusion and decision]

Use one stable name for all further processing and reporting:

After making your role-difference decision, you can reveal the optional
checkpoint and compare whether you reached the same conclusion.

```{r eval=FALSE}
# hint("ild_model_choices", "role_solution")
```

```{r eval=FALSE}
retained_model <- YOUR_RETAINED_MODEL
```


## If pooled: recheck adequacy

If you retained `pooled_apim_effects_model`, repeat the DHARMa and lag-1 checks
with `retained_model`. If you retained the already-checked distinguishable
model, skip this step.

```{r eval=FALSE}
retained_model_dharma_residuals <- simulate_ild_dharma(
  model = retained_model,
  dyad = my_ild_data$coupleID,
  member = my_ild_data$personID,
  time = my_ild_data$diaryday,
  role = my_ild_data$gender,
  n = 250,
  seed = 123
)

plot(
  retained_model_dharma_residuals,
  title = "Retained-model block-whitened residuals"
)

test_ild_lag1(retained_model_dharma_residuals, by_role = TRUE)
```

Briefly assess whether the pooled-effects model remains adequate:

> [Insert your conclusion, or state that you skipped this step]


## Robustness

If time permits, assess a different plausible time trend and compare the focal
estimates. In a substantive analysis, you could also assess theoretically
plausible random slopes if the design supports them. Do not select a
specification solely because it changes a significance label.

> [Insert your robustness conclusion]

For a future nonnegative, right-skewed outcome, `sqrt()` or `log1p()` may
improve the scale; inspect its histogram, refit, and rerun all diagnostics. An
appropriate non-Gaussian family may instead be preferable, but generalized
families have different assumptions and require case-specific choices.


# 6. Visualize, interpret, and report results

Use `retained_model` for every table, plot, and interpretation below.

Print the fixed and random effects, and inspect the fitted covariance
parameters:

```{r}
# hint("ild_reporting", "solution")
```

```{r}
# Add your reporting code here.
```

Interpret each AR(1) SD as the stationary SD of that role's persistent latent
residual process. The dyad--day correlation describes only the remaining
nonpersistent same-occasion component; it is not the correlation of the
complete residual after the AR(1) components are added.

Plot the focal fixed effects or model-implied associations:

```{r}
# hint("ild_reporting", "plot")
```

```{r}
# Add your plotting code here.
```

Write a brief results paragraph that:

- describes the analysis sample and temporal structure
- explains the selected temporal specification
- explains the role-difference decision
- reports within-person and between-person actor and partner effects with
  estimates and 95% confidence intervals
- summarizes stable dyad covariance, same-occasion partner covariance, and any
  fitted AR(1) parameters
- reports the adequacy and robustness checks

Interpret the estimates as associations, not causal effects. If effects are
role-specific, a significant estimate for one role and a nonsignificant
estimate for the other do not by themselves establish that the roles differ.

```{r}
# hint("ild_reporting", "paragraph_template")
```

> [Insert your brief results paragraph]


You have now completed the ILD APIM workflow!
