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

# Exercises

These exercises follow the six-step cross-sectional APIM workflow from the
applied tutorial.

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

# Getting started

An R Markdown (`.Rmd`) file combines *Markdown* text -- such as this paragraph --
with runnable R code chunks.

Run code chunks individually via the small green play button in the
top-right of each chunk. Try it by running the next chunk. It installs or updates
the workshop packages and any dependencies that need updating. This may take
several minutes; you can read ahead while it runs.

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



In RStudio, insert a new R chunk with
Ctrl+Alt+I on Windows or Linux, or Cmd+Option+I on macOS.

Try it and insert a new chunk between the green lines:

---



---


Load all required packages:

```{r}
library(glmmTMB)
library(dyadMLM)
library(easystats, quietly = TRUE) # This package likes to be especially noisy
library(DHARMa)
library(dplyr)

source("hints.R")

```


While going through the exercises you can sometimes spot optional hints:

```{r}
# hint("example_hint", "run_chunk_without_comment")
```

To reveal a hint, delete the `#` in front of the hint you want to read, then
run the chunk (again). You can try it above.



------------ Let's get started! ------------------



# Loading, inspecting, and preparing the data

Follow a simplified workflow:

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

# 1. Define dyads, distinguishability, and research questions

Do not reproduce a predictor--outcome combination from the workshop examples,
such as self-efficacy predicting MVPA or received support predicting MVPA.
Any other substantively sensible combination is fine.

Place `dyadic-person-means.rds` in the same folder as this exercise file, then
load and inspect the raw aggregated dataset:

```{r}
raw_dyad_data <- readRDS("dyadic-person-means.rds")

raw_dyad_data
```


Choose your research question.

> [Insert your actor research question]

> [Insert your partner research question]



# 2. Prepare the data and validate the dyadic structure

Use `dyadMLM::prepare_dyad_data()`. For help with its arguments, run:

```{r eval=FALSE}
?dyadMLM::prepare_dyad_data
```

Adding `?` before a function name opens its help page in RStudio.

---

Prepare and validate the data for a **distinguishable** model now.

- Distinguish partners using `role = gender`.
- Set `include_arbitrary_member_contrast = TRUE` to use the same prepared data
  for distinguishable and exchangeable models later.
- Use `add_apim_gmc_predictors = TRUE` if you want the grand-mean-centered APIM
  predictors used in the tutorial.

Save the prepared data as `my_data`; later exercise code uses this name.

```{r}
# hint("prepare_data", "structure")
```

```{r}
# Prepare the data and save it as my_data.
```

Inspect the prepared data (like you did above with the raw data) and verify the
dyad IDs, member IDs, roles, row counts, missingness, and generated columns
make sense:

```{r}

```


# 3. Describe and visualize measures and observed relationships

Adapt the summary and plotting code from the tutorial to your focal variables.
At minimum, summarize the variables by role, inspect their distributions,
calculate observed dyadic correlations, and visualize partner similarity.

For role-specific summaries and dyadic correlations, create a wide view with
one row per couple using `tidyr::pivot_wider()`. Save it as `my_wide_data`;
continue using `my_data` for modeling.

Summaries:

```{r}
# hint("descriptive_summary", "functions")
# hint("descriptive_summary", "example_code")
```

```{r}
# Create my_wide_data and add your summaries here.
```

Visualize distributions by gender. The easiest way to visualize a distribution
is with `hist()`. You can select rows directly with `[rows, columns]`.

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

```{r}
# Pooled histogram

# Male histogram

# Female histogram
```


Describe partner similarity:

```{r}
# hint("partner_similarity", "correlations")
# hint("partner_similarity", "plot")
```

```{r}
# Correlations

# Plot
```


Describe what you see and how you interpret the summaries and plots:

[Insert your brief summary here]


# 4. Estimate the model

Inspect `names(my_data)` to find the generated columns you need in the formula.
Then estimate a model for distinguishable dyads:

```{r eval=FALSE}
# hint("distinguishable_model", "concept")
# hint("distinguishable_model", "structure")
# hint("distinguishable_model", "example_code")
```

```{r eval=FALSE}
my_model <- glmmTMB(

  YOUR_MODEL_FORMULA,

  data = my_data,
  family = gaussian()
)

summary(my_model)

```



# 5. Verify model adequacy, distinguishability, and robustness

### Adequacy

First, check estimation:

- no convergence warnings
- finite estimates, standard errors, and fit indices
- no variances near zero or correlations near $-1$ or $1$

Then, use DHARMa for residual diagnostics via simulation:

```{r eval=FALSE}
my_model_dharma_residuals <- DHARMa::simulateResiduals(
  fittedModel = my_model,
  n = 10000,

  # Re-simulate the random-effects block that represents the residuals.
  simulateREs = "unconditional",

  # Account for the correlation between partners' residuals.
  rotation = "estimated",
  seed = 123
)

plot(my_model_dharma_residuals)
```

Next, inspect each role separately:

```{r eval=FALSE}
female_dharma_residuals <- DHARMa::recalculateResiduals(
  my_model_dharma_residuals,
  sel = ~ .is_female == 1,
  rotation = NULL,
  seed = 123
)

male_dharma_residuals <- DHARMa::recalculateResiduals(
  my_model_dharma_residuals,
  sel = ~ .is_male == 1,
  rotation = NULL,
  seed = 123
)

# Show both role-specific plots side by side, then restore the previous layout.
old_par <- par(mfrow = c(1, 2))
DHARMa::plotResiduals(female_dharma_residuals, main = "Women")
DHARMa::plotResiduals(male_dharma_residuals, main = "Men")
par(old_par)
```

Interpreting these residual diagnostics takes practice. Consult the DHARMa
vignette section on
[interpreting residuals and recognizing misspecification](https://stat.ethz.ch/CRAN/web/packages/DHARMa/vignettes/DHARMa.html#interpreting-residuals-and-recognizing-misspecification-problems).
For now, briefly assess whether the diagnostics appear acceptable overall and
for each role.

If convergence problems or DHARMa patterns suggest a problematic outcome
scale, inspect a histogram. For a nonnegative, right-skewed outcome, try
`sqrt()` or `log1p()` (which handles zeros), check whether the histogram is less
skewed, refit the model, and rerun DHARMa.

For some outcomes, an appropriate non-Gaussian family may be preferable to a
transformation. Generalized families have different assumptions and require
case-specific choices, so they are not covered in this workshop. If the model
remains inadequate, reconsider its mean and covariance structure before
continuing.


### Distinguishability

Fit the equivalent model under exchangeability constraints. It should use the
same adequate specification and differ only in its exchangeability constraints.
Hint: you need a
`.member_contrast_arbitrary` column in the data. You may already have included
one earlier, or you may need to prepare the data again. Check your column names
with `names(my_data)`.


```{r}
# hint("exchangeable_model", "fixed_effects")
# hint("exchangeable_model", "random_effects")
# hint("exchangeable_model", "example_code")
# hint("outcome_scale", "matched_models")

```

```{r}
# Fit the equivalent exchangeable model here.
```

Compare the distinguishable and exchangeable models using
`dyadMLM::compare_nested_models()`.

```{r}

```

Insert your brief conclusion, and choose whether you want to continue with
the distinguishable or exchangeable model for the rest of the exercise:

[Insert conclusion]

[Insert decision and reasoning]

You may also test selected equality constraints rather than full exchangeability.
If you do, compare that model with the distinguishable model and include it in
your retention decision.

Save the model you retained as `retained_model`:

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


### If exchangeable: recheck adequacy

If you retained the distinguishable model unchanged, its adequacy was already
checked above. If you retained the exchangeable model, repeat the overall
DHARMa check:

```{r eval=FALSE}
retained_model_dharma_residuals <- DHARMa::simulateResiduals(
  fittedModel = retained_model,
  n = 10000,

  # Re-simulate the random-effects block that represents the residuals.
  simulateREs = "unconditional",

  # Account for the correlation between partners' residuals.
  rotation = "estimated",
  seed = 123
)

plot(retained_model_dharma_residuals)
```

### Robustness

Inspect unusual dyads in separate role panels. Use the same panels if you
retained the exchangeable model: the split simply keeps one observation per
dyad in each plot and does not imply role-specific model effects. Use the same
predictor and outcome versions as in your retained model. For example, use the
generated actor or partner predictor column and, if applicable, the transformed
outcome column.

Optionally, you can use the supplied helper that is loaded here from a separate
file `helpers.R`:

```{r}
# hint("robustness_mahalanobis", "columns")
```

```{r eval=FALSE}
# The function is defined in the separate helper file.
source("helpers.R")

women_mahalanobis_plot <- plot_mahalanobis(
  data = dplyr::filter(my_data, gender == "female"),
  x = "YOUR_ACTOR_PREDICTOR_COLUMN",
  y = "YOUR_OUTCOME_COLUMN",
  title = "Women"
)

men_mahalanobis_plot <- plot_mahalanobis(
  data = dplyr::filter(my_data, gender == "male"),
  x = "YOUR_ACTOR_PREDICTOR_COLUMN",
  y = "YOUR_OUTCOME_COLUMN",
  title = "Men"
)

see::plots(women_mahalanobis_plot, men_mahalanobis_plot)
```

Repeat the same for partner effects:

```{r}

```


Mahalanobis distance identifies unusual bivariate combinations, not formal
model influence, and the displayed regression line is not the adjusted APIM
effect. Use the labels to guide data inspection rather than to automatically
exclude dyads.

If a dyad was flagged for a substantive or procedural reason, optionally refit
the retained model without it:

```{r eval=FALSE}
# hint("robustness_refit", "ids")
```

```{r eval=FALSE}
data_without_selected_dyads <- dplyr::filter(
  my_data,
  !couple_id %in% DYAD_IDS_TO_EXCLUDE
)

model_without_selected_dyads <- update(
  retained_model,
  data = data_without_selected_dyads
)

summary(model_without_selected_dyads)
```

Do not exclude dyads merely because their removal changes a significance label.

Compare the focal estimates and their uncertainty with those from
`retained_model`. State whether any substantive conclusion changed:

> [Insert your robustness conclusion]


# 6. Visualize, interpret, and report results

Print a fixed-effects table for `retained_model` using
`parameters::model_parameters()`:

```{r}
# hint("report_results", "fixed_effects")
```

```{r}
# Add your fixed-effects table here.
```

Print a random-effects table:

```{r}
# hint("report_results", "random_effects")
```

```{r}
# Add your random-effects table here.
```

Plot the fixed effects:

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

```{r}
# Add your fixed-effects plot here.
```


Write a brief results paragraph that:

- describes the analysis sample and any exclusions or missing-data handling
- identifies the retained model and explains the distinguishability decision
- answers the actor and partner research questions using estimates and 95%
  confidence intervals, clearly stating the modeled outcome scale
- summarizes the fitted residual covariance, adequacy checks, and robustness
  analysis

Remember that the interpretation of the intercepts depends on whether you used
grand-mean-centered variables.

If you transformed the outcome, state the modeled scale and interpret estimates
on that scale—or use back-transformed predictions.

> [Insert your brief results paragraph]


You have now completed the APIM workflow!


## Optional: Compute DIM estimates from the exchangeable APIM

The transformation from exchangeable APIM coefficients to DIM estimates is
particularly simple.

Check the output from your **exchangeable** model, and enter your coefficients
here:

```{r eval=FALSE}

b_actor <- YOUR_ACTOR_COEFFICIENT
b_partner <- YOUR_PARTNER_COEFFICIENT

```

Recover the DIM estimates:

```{r eval=FALSE}
# hint("dim_transformation", "formulas")
```

```{r eval=FALSE}
b_mean <- YOUR_FORMULA
b_dev <- YOUR_FORMULA

cat(
  "The between-dyad effect is: ", round(b_mean, 2), "\n",
  "The within-dyad effect is: ", round(b_dev, 2)
)

```


Interpret the effects:

[Insert interpretation]


Optionally, fit a DIM or DSM directly using the package vignettes:

- [DIM vignette](https://pascal-kueng.github.io/dyadMLM/articles/dim.html)
- [DSM vignette](https://pascal-kueng.github.io/dyadMLM/articles/dsm.html)


# Intensive longitudinal dyads

You may continue with the
[intensive longitudinal exercises](04_exercises_ild.Rmd).
