---
title: "Chapter 10: Moderation, mediation, and missingness"
author: "Andy Field"
date: "2025-04-18"
date-modified: today
format: 
  html:
    theme: journal
    embed-resources: true
    toc: true
    code-overflow: wrap
execute: 
    message: false
    warning: false
---

```{r}
#| label: setup

library(easystats)
library(tidyverse)

# load data
# video game data
vids_tib  <- here::here("data/video_games.csv") |> 
  read_csv()
# infidelity data
infidelity_tib  <- here::here("data/lambert_2012.csv") |> 
  read_csv()
# fake news data
newz_tib  <- here::here("data/bronstein_2019.csv") |> 
  read_csv()
# more fake news data
newz_md_tib  <- here::here("data/bronstein_miss_2019.csv") |> 
  read_csv()


# or load from the discovr package
# vids_tib <- discovr::video_games
# infidelity_tib <- discovr::lambert_2012
# newz_tib <- discovr::bronstein_2019
# newz_md_tib <- discovr::bronstein_miss_2019
```


## Centering variables

```{r}
vids_cent_tib <- center(x = vids_tib, select = c(vid_game, caunts))
```

## Fitting the model

```{r}
aggress_lm <- lm(aggress ~ caunts + vid_game + caunts:vid_game, data = vids_cent_tib)
```

Or

```{r eval = FALSE}
aggress_lm <- lm(aggress ~ caunts*vid_game, data = vids_cent_tib)
```

Viewing the model:

```{r}
# fit
model_performance(aggress_lm) |> 
  display()
# test of fit
test_wald(aggress_lm) |> 
  display()
# parameter estimates
model_parameters(aggress_lm) |> 
  display(digits = 3)
```


Fit a robust model:

```{r}
model_parameters(aggress_lm, vcov = "HC4") |> 
  display(digits = 3)
```


## Simple slopes and Johnson-Neyman interval

```{r}
# Johnson-Neyman interval
vid_jn <- estimate_slopes(model = aggress_lm,
                          trend = "vid_game",
                          by = "caunts",
                          length = 100)
summary(vid_jn) |> 
  display()
plot(vid_jn) + theme_bw()

# simple slopes
vid_ss <- estimate_slopes(model = aggress_lm, trend = "vid_game", by = "caunts=[sd]")
display(vid_ss)

# plot simple slopes:
estimate_means(model = aggress_lm, by = c("vid_game", "caunts=[sd]")) |> 
  plot() +
  labs(x = " Video game use per week (centred) ", y = "Aggression", colour = "Callous traits (centred)", fill = "Callous traits (centred)") +
  theme_minimal()

```


# Mediation
## Self-test: the bad way

```{r}

m1 <- lm(phys_inf ~ ln_porn, data = infidelity_tib)
m2 <- lm(commit ~ ln_porn, data = infidelity_tib)
m3 <- lm(phys_inf ~ ln_porn + commit, data = infidelity_tib)

model_parameters(m1) |>  display()
model_parameters(m2) |>  display()
model_parameters(m3) |>  display()
```

## The better method

```{r}
infidelity_mod <- 'phys_inf ~ c*ln_porn + b*commit
                   commit ~ a*ln_porn

                   indirect_effect := a*b
                   total_effect := c + (a*b)
                   '

infidelity_fit <- lavaan::sem(infidelity_mod,
                              data = infidelity_tib,
                              missing = "FIML",
                              estimator = "MLR")
model_performance(infidelity_fit) |>
  display() 
model_parameters(infidelity_fit) |>  
  display()
```

Non-robust version for kicks and giggles:

```{r}
infidelity_nonrob <- lavaan::sem(infidelity_mod, data = infidelity_tib, missing = "FIML")

model_parameters(infidelity_nonrob) |>  
  display()
```

Standardized solution:

```{r}
model_parameters(infidelity_nonrob, standardize = TRUE) |>  
  display()
```


## Two mediators

```{r}
newz_mod <- 'fake_newz ~ c*delusionz + b1*thinkz_open + b2*thinkz_anal
                thinkz_open ~ a1*delusionz
                thinkz_anal ~ a2*delusionz

                indirect_open := a1*b1
                indirect_anal := a2*b2
                total_effect := c + (a1*b1) + (a2*b2)
                thinkz_open ~~ thinkz_anal
                '

newz_fit <- lavaan::sem(newz_mod, data = newz_tib, se = "bootstrap")
model_parameters(newz_fit) |> 
  display()
```

# Missing data

## Look for patterns of missingness

```{r}
newz_mdpattern <- newz_md_tib |>
  select(-id) |> 
  mice::md.pattern(rotate.names = TRUE)

display(newz_mdpattern)

```

```{r}
mice::ic(newz_md_tib) |> 
  display()
```


## Listwise deletion (don't do it!)

```{r}

fake_lm <- lm(fake_newz ~ 0 + delusionz + thinkz_open + thinkz_anal, data = newz_md_tib)
n_obs(fake_lm)
model_parameters(fake_lm) |> display()
```

## FIML

```{r}

newz_lm <- 'fake_newz ~ delusionz + thinkz_open + thinkz_anal'

# If you want the intercept include 1
# newz_lm <- 'fake_newz ~ 1 + delusionz + thinkz_open + thinkz_anal'

newz_lm_fit <- lavaan::sem(newz_lm, data = newz_md_tib, missing = "fiml.x")

model_parameters(newz_lm_fit) |> 
  display()
```

## Multiple imputation

### Find auxillary variables

```{r}
aux_vars <- mice::quickpred(newz_md_tib, exclude = "id", mincor = 0.1)
display(aux_vars)
```


### Imputation

```{r}

newz_imps <- mice::mice(newz_md_tib, pred = aux_vars, method = "pmm", m = 20)

newz_mi_fit <- with(newz_imps, lm(fake_newz ~ 0 + delusionz + thinkz_open + thinkz_anal))
newz_pool <- mice::pool(newz_mi_fit)

model_parameters(newz_pool) |> 
  display()
```

