03-analysis.Rmd

# Main analysis

This document contains the main analyses presented in the manuscript.

We first load the required R packages.

```{r packages}
#| results: 'hide'
#| cache: FALSE

library(knitr)
library(scales)
library(DT)
library(kableExtra)
library(ggpp)
library(bayestestR)
library(brms)
library(ragg)
library(here)
library(ggtext)
library(ggstance)
library(ggdist)
library(memoise)
library(cachem)
library(patchwork)
library(lavaan)
library(multidplyr)
library(tidyverse)
library(tidybayes)
library(janitor)
cd <- cache_disk("memoise")
```

Define plotting, table, and parallel processing options.

```{r settings}
#| cache: FALSE

# parallel computations
MAX_CORES <- as.numeric(Sys.getenv("MAX_CORES"))
if (is.na(MAX_CORES)) MAX_CORES <- parallel::detectCores(logical = FALSE)
cluster <- new_cluster(MAX_CORES)

# load packages on clusters
cluster_library(cluster, c("dplyr", "lavaan"))

# MCMC settings
options(mc.cores = 1)
if (require("cmdstanr")) options(brms.backend = "cmdstanr")

# Save outputs here
dir.create("models", FALSE)

# Functions in R/
source(here("R", "meta-analysis-helpers.R"))
```

We analyse the data cleaned previously.

```{r data-load}
data_path <- here("Data", "cleaned_data.rds")
if (file.exists(data_path)) {
  d <- read_rds(file = data_path)
} else {
  stop(str_glue("{data_path} doesn't exist, run `01-clean.Rmd` to create it."))
}

# Make wave a nicely labelled factor
d <- d %>%
  mutate(Wave = factor(wid, levels = 1:3, labels = paste0("Wave ", 1:3)))

# Rename game to fit titles in plots
d <- d %>% 
  mutate(Game = if_else(Game == "Gran Turismo Sport", "GT Sport", Game))
```

For all analyses, game time indicates average number of hours per day. The telemetry indicates total hours in the two week window, so we divide that by 14 to talk about average hours per day.

```{r hours-per-day-transform}
d <- d %>%
  mutate(
    Hours = Hours / 14,
    hours_est = hours_est / 14
  )
```

## Variables over time

Hours played per week is cut at a reasonable value for the figures:

```{r}
d %>%
  mutate(Hours_over_3 = Hours > 3) %>%
  tabyl(Hours_over_3) %>%
  adorn_pct_formatting() %>% 
  kbl(caption = "Table of hours excluded from figure") %>% 
  kable_classic(full_width = FALSE, html_font = Font)
```

```{r variables-over-time}
#| fig.height: 7
#| fig.cap: >
#|  Density plots of key variables over time.

tmp <- d %>%
  select(
    Game, Wave, pid,
    Hours, Affect, `Life satisfaction`,
    Intrinsic, Extrinsic
  ) %>%
  pivot_longer(Hours:Extrinsic) %>%
  drop_na(value) %>%
  filter(!(name == "Hours" & value > 3)) %>%
  mutate(name = if_else(name == "Hours", "Hours / day", name)) %>% 
  mutate(name = fct_inorder(name))
tmp %>%
  ggplot(
    aes(
      Wave,
      value
    )
  ) +
  scale_x_discrete(labels = 1:3, expand = expansion(c(0.1, .1))) +
  scale_y_continuous(
    "Value",
    breaks = pretty_breaks(),
    expand = expansion(.025)
  ) +
  geom_blank() +
  stat_halfeye(
    height = .02,
    normalize = "panels",
    slab_color = col1,
    slab_fill = alpha(col1, 0.3),
    slab_size = 0.5,
    adjust = 1.1,
    point_interval = NULL,
    show.legend = FALSE
  ) +
  stat_summary(
    fun.data = mean_cl_normal,
    fatten = 1.25
  ) +
  stat_summary(
    fun = mean,
    geom = "line",
    size = .33,
    group = 1
  ) +
  facet_grid(name ~ Game, scales = "free_y") +
  theme(
    legend.position = "none"
  )
```


```{r, include = FALSE}
agg_png(
  "Figures/Fig-variables-over-time.png", 
  width = W, height = W*0.7, units = "cm", res = 400
)
last_plot()
dev.off()
```

Then take a look at a simple model of change over time for each variable. Note that we can't use varying slopes with lmer because there's not enough data, so just random intercepts for players.

```{r model-variables-over-time}
library(lme4)
library(broom.mixed)
path <- "models/parameters_change_over_time.rds"
if (!file.exists(path)) {
  parameters_change_over_time <- d %>%
    select(
      Game, pid, wid,
      Hours, Intrinsic, Extrinsic,
      Affect, `Life satisfaction`
    ) %>%
    # Put intercept at first wave
    mutate(Wave = wid - 1) %>%
    pivot_longer(Hours:`Life satisfaction`, names_to = "Variable") %>%
    group_by(Variable) %>%
    summarise(
      lmer(value ~ Wave + (1 | pid) + (1 + Wave | Game), data = cur_data()) %>%
        tidy(., "fixed", conf.int = TRUE)
    )
  saveRDS(parameters_change_over_time, path)
} else {parameters_change_over_time <- readRDS(path)}

parameters_change_over_time %>%
  mutate(across(where(is.numeric), ~ format(round(.x, 2), nsmall = 2))) %>%
  mutate(Result = str_glue("{estimate}, 95%CI [{conf.low}, {conf.high}]")) %>%
  select(Variable, term, Result) %>%
  pivot_wider(names_from = term, values_from = Result) %>% 
  kbl(caption = "Change over time parameter estimates") %>% 
  kable_classic(full_width = FALSE, html_font = Font)
```

## Simple correlation

Before more informative modelling, we estimate simple bivariate regressions between the key variables. 

```{r simple-regression-coefficients}
#| fig.height: 7
#| fig.cap: >
#|  Unstandardised bivariate regression coefficients of models predicting well-being. Columns indicate predictors, and the two rows are the different well-being outcome variables.

tmp <- d %>% 
  arrange(pid, wid) %>% 
  group_by(pid) %>% 
  mutate(h1 = lag(Hours))

fit <- lmer(Affect ~ h1 + (1 | pid) + (h1 | Game), data = tmp)
summary(fit)

d %>%
  select(
    Game, pid, wid, Wave,
    Hours, hours_est, Intrinsic, Extrinsic,
    Affect, `Life satisfaction`
  ) %>%
  pivot_longer(
    c(Affect, `Life satisfaction`),
    names_to = "Outcome", values_to = "Outcome_value"
  ) %>%
  pivot_longer(
    c(Hours, hours_est, Intrinsic, Extrinsic),
    names_to = "Predictor", values_to = "Predictor_value"
  ) %>%
  group_by(Game, Outcome, Predictor, Wave) %>%
  summarise(
    tidy(
      lm(
        Outcome_value ~ Predictor_value, 
        data = cur_data()
      ), 
      conf.int = TRUE
    ),
  ) %>%
  filter(term != "(Intercept)") %>%
  ggplot(aes(estimate, Game, col = Wave)) +
  geom_vline(xintercept = 0, lty = 2, size = .25) +
  scale_x_continuous(
    "Bivariate regression coefficient (95%CI)",
    breaks = pretty_breaks()
  ) +
  geom_pointrangeh(
    aes(xmin = conf.low, xmax = conf.high),
    size = .4,
    position = position_dodge2v(.35)
  ) +
  facet_grid(Outcome ~ Predictor) +
  theme(
    legend.position = "bottom",
    axis.title.y = element_blank()
  )
```

## RICLPM

Wrangle data to a format where lavaan model is easier to map to variable pairs: Wide format with different rows for outcomes (well-being) and predictors (hours, needs, motivations).

```{r riclpm-data-transform}
d_riclpm_long <- d %>%
  select(
    Game, pid, wid,
    Hours, hours_est,
    Intrinsic, Extrinsic,
    Affect, `Life satisfaction`, AN, AP
  ) %>%
  # Long format on well-being scores (outcomes)
  pivot_longer(
    c(Affect, `Life satisfaction`, AN, AP),
    names_to = "y_var", values_to = "y"
  ) %>%
  # Long format on other variables (predictors)
  pivot_longer(
    c(Hours, hours_est, Intrinsic, Extrinsic),
    names_to = "x_var", values_to = "x"
  )
d_riclpm_wide <- d_riclpm_long %>%
  pivot_wider(names_from = wid, values_from = c(x, y), names_sep = "")
write_rds(d_riclpm_wide, here("Temp", "d_riclpm_wide.rds"))
```

### Cross-lagged scatterplots

```{r riclpm-crosslag-scatter-hours}
#| fig.height: 5
#| fig.cap: >
#|  Scatterplots of well-being (rows indicate variables and waves) on average hours played per day during the previous wave. Regression lines are GAM fitted lines and 95%CIs.

Order <- tibble(
  y_var = rep(c("Affect", "Life satisfaction"), each = 2),
  wid = c(2, 3, 2, 3),
  Panel = c(
    "Affect[plain('[Wave 2]')]",
    "Affect[plain('[Wave 3]')]",
    "LS[plain('[Wave 2]')]",
    "LS[plain('[Wave 3]')]"
  ),
  Panel_label = fct_inorder(Panel)
)

d_riclpm_plots <- d_riclpm_long %>%
  # Don't show AN and AP
  filter(!(y_var %in% c("AN", "AP"))) %>% 
  # Take out hours per day over 3 for these plots
  mutate(x = if_else(str_detect(x_var, "ours") & x > 3, NaN, x)) %>%
  arrange(Game, pid, x_var, wid) %>%
  group_by(Game, pid, x_var, y_var) %>%
  mutate(lag_x = lag(x)) %>%
  left_join(Order) %>%
  ungroup() %>%
  filter(wid > 1)

d_riclpm_plots %>%
  filter(x_var == "Hours") %>%
  ggplot(aes(lag_x, y)) +
  scale_x_continuous(
    "Hours played per day at previous wave",
    breaks = pretty_breaks(3)
  ) +
  scale_y_continuous(
    "Well-being at current wave",
    breaks = pretty_breaks()
  ) +
  geom_point(size = .2, alpha = .2, shape = 1, col = col1) +
  geom_smooth(
    method = "gam", size = .4,
    color = "black",
    alpha = .33, show.legend = FALSE
  ) +
  facet_grid(
    Panel_label ~ Game,
    scales = "free_y",
    labeller = labeller(.rows = label_parsed)
  ) +
  theme(
    aspect.ratio = 1,
    panel.grid = element_blank()
  )
```


```{r, include = FALSE}
agg_png(
  "Figures/Fig-cl-scatter-hours.png", 
  width = W, height = W*0.6, units = "cm", res = 400
)
last_plot()
dev.off()
```

```{r riclpm-crosslag-scatter-hours-estimated}
#| fig.height: 5
#| fig.cap: >
#|  Scatterplots of well-being (rows indicate variables and waves) on estimated average hours played per day during the previous wave. Regression lines are GAM fitted lines and 95%CIs.

# Plots with other predictors
p_tmp_0 <- last_plot() %+%
  filter(d_riclpm_plots, x_var == "hours_est") +
  scale_x_continuous(
    "Estimated hours played per day at previous wave",
    breaks = 0:3
  )
p_tmp_0
```

```{r riclpm-crosslag-scatter-intrinsic}
#| fig.height: 5
#| fig.cap: >
#|  Scatterplots of well-being (rows indicate variables and waves) on intrinsic motivation during the previous wave. Regression lines are GAM fitted lines and 95%CIs.

p_tmp_1 <- last_plot() %+%
  filter(d_riclpm_plots, x_var == "Intrinsic") +
  scale_x_continuous(
    "Intrinsic motivation at previous wave",
    breaks = pretty_breaks()
  )
p_tmp_1
```

```{r riclpm-crosslag-scatter-extrinsic}
#| fig.height: 5
#| fig.cap: >
#|  Scatterplots of well-being (rows indicate variables and waves) on extrinsic motivation during the previous wave. Regression lines are GAM fitted lines and 95%CIs.

p_tmp_2 <- last_plot() %+%
  filter(d_riclpm_plots, x_var == "Extrinsic") +
  scale_x_continuous(
    "Extrinsic motivation at previous wave",
    breaks = pretty_breaks()
  )
p_tmp_2
```

```{r riclpm-crosslag-scatter-motivations}
#| fig.height: 9
#| fig.cap: >
#|  Scatterplots of well-being (rows indicate variables and waves) on intrinsic and extrinsic motivation during the previous wave. Regression lines are GAM fitted lines and 95%CIs.

p_tmp_1 / p_tmp_2
```


```{r, include = FALSE}
agg_png(
  "Figures/Fig-cl-scatter-motivations.png", 
  width = W, height = W*1.1, units = "cm", res = 400
)
last_plot()
dev.off()
```

### Fit models

First fit RICLPM to all data. Separately to the two well-being outcomes, and each predictor (hours, subjective scales). Constrain (cross)lagged parameters

Syntax based on <https://jeroendmulder.github.io/RI-CLPM/lavaan.html>

```{r riclpm-define-model}
riclpm_1 <- "
  # Create between components (random intercepts)
  RIx =~ 1*x1 + 1*x2 + 1*x3
  RIy =~ 1*y1 + 1*y2 + 1*y3

  # Create within-person centered variables
  wx1 =~ 1*x1
  wx2 =~ 1*x2
  wx3 =~ 1*x3
  wy1 =~ 1*y1
  wy2 =~ 1*y2
  wy3 =~ 1*y3

  # Estimate the lagged effects between the within-person centered variables (constrained).
  wx2 ~ bx*wx1 + gx*wy1
  wy2 ~ gy*wx1 + by*wy1
  wx3 ~ bx*wx2 + gx*wy2
  wy3 ~ gy*wx2 + by*wy2

  # Estimate the covariance between the within-person centered
  # variables at the first wave.
  wx1 ~~ wy1 # Covariance

  # Estimate the covariances between the residuals of the
  # within-person centered variables (the innovations).
  wx2 ~~ wy2
  wx3 ~~ wy3

  # Estimate the variance and covariance of the random intercepts.
  RIx ~~ RIx
  RIy ~~ RIy
  RIx ~~ RIy

  # Estimate the (residual) variance of the within-person centered variables.
  wx1 ~~ wx1 # Variances
  wy1 ~~ wy1
  wx2 ~~ wx2 # Residual variances
  wy2 ~~ wy2
  wx3 ~~ wx3
  wy3 ~~ wy3
"
```

```{r riclpm-define-model-alt}
# This is the same but was suggested to not include autoregressions
riclpm_2 <- "
  # Create between components (random intercepts)
  RIx =~ 1*x1 + 1*x2 + 1*x3
  RIy =~ 1*y1 + 1*y2 + 1*y3

  # Create within-person centered variables
  wx1 =~ 1*x1
  wx2 =~ 1*x2
  wx3 =~ 1*x3
  wy1 =~ 1*y1
  wy2 =~ 1*y2
  wy3 =~ 1*y3

  # Estimate the lagged effects between the within-person centered variables (constrained). Here, we removed the autoregressions.
  wx2 ~ gx*wy1
  wy2 ~ gy*wx1
  wx3 ~ gx*wy2
  wy3 ~ gy*wx2

  # Estimate the covariance between the within-person centered
  # variables at the first wave.
  wx1 ~~ wy1 # Covariance

  # Estimate the covariances between the residuals of the
  # within-person centered variables (the innovations).
  wx2 ~~ wy2
  wx3 ~~ wy3
  
  # Here, we also included the correlated residuals as suggested:
  wx2 ~~ wx1
  wx3 ~~ wx2
  wy2 ~~ wy1
  wy3 ~~ wy2

  # Estimate the variance and covariance of the random intercepts.
  RIx ~~ RIx
  RIy ~~ RIy
  RIx ~~ RIy

  # Estimate the (residual) variance of the within-person centered variables.
  wx1 ~~ wx1 # Variances
  wy1 ~~ wy1
  wx2 ~~ wx2 # Residual variances
  wy2 ~~ wy2
  wx3 ~~ wx3
  wy3 ~~ wy3
"
```

```{r riclpm-estimate-pooled}
cluster_library(cluster, "purrr")
cluster_copy(cluster, c("riclpm_1", "riclpm_2"))
path <- "models/riclpm-pooled.rds"
if (!file.exists(path)) {
  fit_riclpm_all <- d_riclpm_wide %>%
    group_by(y_var, x_var) %>%
    nest() %>% 
    crossing(
      nesting(
        lvn = c(riclpm_1, riclpm_2), 
        model = c("1", "2")
      )
    ) %>% 
    partition(cluster) %>%
    mutate(
      fit = map2(
        data, lvn,
        ~lavaan(
          .y,
          data = .x,
          missing = "ml",
          meanstructure = TRUE,
          int.ov.free = TRUE
        )
      )
    ) %>%
    collect() %>% 
    ungroup() %>% 
    select(-data, -lvn)
  saveRDS(fit_riclpm_all, path)
} else {fit_riclpm_all <- readRDS(path)}
```

Fit the above model separately to each game.

```{r riclpm-estimate-separate}
path <- "models/riclpm-games.rds"
if (!file.exists(path)) {
  fit_riclpm_games <- d_riclpm_wide %>%
    group_by(y_var, x_var, Game) %>%
    nest() %>% 
    crossing(
      nesting(
        lvn = c(riclpm_1, riclpm_2), 
        model = c("1", "2")
      )
    ) %>% 
    partition(cluster) %>%
    mutate(
      fit = map2(
        data, lvn,
        ~lavaan(
          .y,
          data = .x,
          missing = "ml",
          meanstructure = TRUE,
          int.ov.free = TRUE
        )
      )
    ) %>%
    collect() %>% 
    ungroup() %>% 
    select(-data, -lvn)
  saveRDS(fit_riclpm_games, path)
} else {fit_riclpm_games <- readRDS(path)}
```

```{r riclpm-parameters}
get_lavaan_pars <- function(x) {
  bind_rows(
    parameterestimates(x) %>%
      mutate(Type = "Unstandardized"),
    standardizedsolution(x) %>%
      rename(est = est.std) %>%
      mutate(Type = "Standardized")
  ) %>%
    as_tibble() %>%
    unite("Parameter", c(lhs, op, rhs), sep = " ", remove = FALSE)
}

pars_riclpm <- bind_rows(
  "Pooled" = fit_riclpm_all,
  "Independent" = fit_riclpm_games,
  .id = "Fit"
) %>% 
  mutate(pars = map(fit, get_lavaan_pars)) %>%
  select(-fit)

# Take only parameters of interest
pars_riclpm <- pars_riclpm %>% 
  unnest(pars) %>% 
    filter(
    (str_detect(Parameter, " ~ ") & !str_detect(Parameter, "3")) |
      Parameter %in% c("RIx ~~ RIy", "wx1 ~~ wy1", "wx2 ~~ wy2", "wx3 ~~ wy3")
  ) %>% 
  select(-c(lhs:label, z))

pars_riclpm %>%
  filter(Fit == "Pooled", model == "1") %>% 
  filter(
    str_detect(Parameter, " ~ "),
    Type == "Unstandardized"
  ) %>%
  mutate(across(where(is.numeric), ~ format(round(.x, 2), nsmall = 2))) %>%
  mutate(Result = str_glue("{est}, [{ci.lower}, {ci.upper}]")) %>%
  select(y_var, x_var, Parameter, Result) %>%
  pivot_wider(names_from = Parameter, values_from = Result) %>%
  arrange(y_var, x_var) %>% 
  kbl(
    caption = "RICLPM regression parameters from pooled model."
  ) %>% 
  kable_classic(full_width = FALSE, html_font = Font)
```

### Table of parameters

This table displays the lavaan estimates (game specific and pooled) for model 1 (full RICLPM) and 2 (RICLPM without autocorrelations).

```{r}
# Interactive table so reader can find what they're looking for
pars_riclpm %>% 
  mutate(
    pvalue = pvalue(pvalue),
    across(c(model, x_var, y_var, Game, Parameter, Type, Fit), factor),
    across(where(is.numeric), ~format(round(.x, 2), nsmall = 2)),
    Estimate = str_glue("{est} [{ci.lower}, {ci.upper}]")
  ) %>% 
  select(
    model, Fit, x_var, y_var, Game, Parameter, Type,
    Estimate, pvalue
  ) %>% 
  datatable(
    filter = "top",
    class = "display",
    rownames = FALSE
  ) %>% 
  formatStyle(TRUE, `font-size` = '12px')
```

## Meta-analyses

Number of observations/participants per model

```{r table-n-per-meta-analysis}
fit_riclpm_games %>%
  filter(model == 1) %>% 
  mutate(N = map_dbl(fit, nobs) %>% comma()) %>%
  select(-fit, -model) %>%
  pivot_wider(
    names_from = x_var,
    values_from = N,
    names_glue = "{x_var} {.value}"
  ) %>%
  rename(Outcome = y_var) %>%
  arrange(Outcome, Game) %>%
  kbl(caption = "Sample sizes per RICLPM") %>% 
  kable_classic(full_width = FALSE, html_font = Font)
```

Then conduct all meta-analyses (for each model, y_var-x_var pair, parameter of interest, and type of parameter (standardized, unstandardized)) of the game-specific RICLPM parameters.

```{r fit-meta-analyses}
#| results: 'hide'
#| cache: FALSE

# Select variables from RICLPM parameters to data for MA
d_ma <- pars_riclpm %>%
  filter(
    Fit == "Independent",
    str_detect(Parameter, " ~ ") | 
      Parameter %in% c("wx1 ~~ wy1", "RIx ~~ RIy")
  ) %>%
  select(
    model, Type, y_var, 
    x_var, Game, Parameter, 
    est, se, ci.lower, ci.upper
  )
path <- "models/all-meta-analyses.rds"
if (!file.exists(path)) {
  cluster_copy(cluster, "get_ma_post")
  cluster_library(
    cluster, 
    c("brms", "purrr", "cmdstanr", "rstan", "bayestestR", "forcats")
    )
  fit_ma <- d_ma %>%
    mutate(
      name = str_glue("models/brm-{model}-{Type}-{y_var}-{x_var}-{Parameter}")
    ) %>% 
    group_by(model, x_var, y_var, Parameter, Type, name) %>%
    nest() %>% 
    partition(cluster) %>%
    mutate(
      # Fit model
      fit = map2(
        data, name,
        ~brm(
          bf(est | se(se) ~ 1 + (1 | Game)),
          data = .x,
          backend = "cmdstanr",
          prior = prior(student_t(7, 0, 0.25), class = "sd", group = "Game"),
          control = list(adapt_delta = .9999, max_treedepth = 15),
          iter = 15000, warmup = 10000,
          file = .y
        )
      ),
      # Table of posterior samples
      posterior = map(fit, get_ma_post),
      # Table of posterior summary
      summary = map(
        posterior,
        ~ describe_posterior(
          .x, 
          centrality = "mean", 
          ci = 0.95, 
          ci_method = "eti",
          test = "pd"
        ) %>%
          rename(Game = Parameter) %>% 
          mutate(Game = fct_relevel(Game, "Average")) %>% 
          as_tibble()
      ),
      # 
      sd = map(
        fit, 
        ~posterior_summary(.x, variable = "sd_Game__Intercept") %>% 
          as_tibble()
      )
    ) %>% 
    collect() %>%
    ungroup() %>% 
    select(-name, -fit)
  saveRDS(fit_ma, path, compress = FALSE)
} else {fit_ma <- readRDS(path)}

# Arrange for tables etc
fit_ma <- fit_ma %>% 
  mutate(
    model = factor(model, levels = 1:2, labels = c("AR", "No AR")),
    Type = factor(Type, levels = c("Unstandardized", "Standardized")),
    y_var = factor(
      y_var, levels = c("Affect", "Life satisfaction", "AN", "AP")
    ),
    x_var = factor(
      x_var, levels = c("Hours", "Intrinsic", "Extrinsic", "hours_est")
    ),
    Parameter = fct_relevel(Parameter, "wy2 ~ wx1")
  ) %>% 
  arrange(model, Type, y_var, x_var, Parameter)
```

### Table of parameters

All the resulting parameters are listed in the table below.

```{r meta-analysis-parameter-table}
fit_ma %>% 
  select(model:Parameter, summary) %>% 
  unnest(summary) %>% 
  mutate(
    across(where(is.numeric), ~format(round(.x, 2), nsmall = 2)),
    Estimate = str_glue("{Mean} [{CI_low}, {CI_high}]")
  ) %>% 
  select(model, Type, y_var, x_var, Parameter, Game, Estimate, pd) %>% 
  datatable(
    filter = "top",
    class = "display",
    rownames = FALSE
  ) %>% 
  formatStyle(TRUE, `font-size` = '12px')
```

### Hours <-> WB

Draw a forest plot of meta-analyses with hours. First, define informative labels and correct order for them, and load a function for drawing forest plots.

```{r forest-plot-labels}
# Table of nice labels in proper order
Order <- distinct(fit_ma, x_var, y_var, Parameter) %>%
  ungroup() %>%
  # Meta-analyses only for regression parameters
  filter(str_detect(Parameter, " ~ ")) %>%
  mutate(
    x_lab = case_when(
      x_var == "Hours" ~ "Play",
      x_var == "Intrinsic" ~ "Intrinsic",
      x_var == "Extrinsic" ~ "Extrinsic"
    ),
    y_lab = case_when(
      y_var == "Affect" ~ "Affect",
      y_var == "Life satisfaction" ~ "Life~satisfaction"
    ),
    Panel_label = case_when(
      Parameter == "wx2 ~ wx1" ~
        str_glue('{x_lab}[plain("[t-1]")]%->%{x_lab}[plain("[t]")]'),
      Parameter == "wy2 ~ wy1" ~
        str_glue('{y_lab}[plain("[t-1]")]%->%{y_lab}[plain("[t]")]'),
      Parameter == "wy2 ~ wx1" ~
        str_glue('{x_lab}[plain("[t-1]")]%->%{y_lab}[plain("[t]")]'),
      Parameter == "wx2 ~ wy1" ~
        str_glue('{y_lab}[plain("[t-1]")]%->%{x_lab}[plain("[t]")]')
    )
  ) %>%
  mutate(
    Parameter_order = factor(
      Parameter,
      levels = c("wy2 ~ wx1", "wx2 ~ wy1")
    )
  ) %>%
  arrange(Parameter_order, y_var) %>%
  mutate(Panel_label = fct_inorder(Panel_label))

# Only these are needed
forest_plot_data <- fit_ma %>% 
  filter(
    model == "AR", 
    Type == "Unstandardized",
    y_var %in% c("Affect", "Life satisfaction"),
    x_var %in% c("Hours", "Intrinsic", "Extrinsic"),
    Parameter %in% c("wy2 ~ wx1", "wx2 ~ wy1")
  ) %>% 
  left_join(Order) %>% 
  select(-model, -Type)
```

```{r forest-hours-wb}
#| fig.height: 5
#| fig.cap: "Plot"

forest_plot_data %>% 
  filter(x_var == "Hours") %>% 
  forest_plot()
```


```{r, include = FALSE}
agg_png(
  "Figures/Fig-forest-hours-wb.png", 
  width = W, height = W*0.6, units = "cm", res = 400
)
last_plot()
dev.off()
```

```{r forest-hours-wb-lavaan}
#| fig.height: 5
#| fig.cap: >
#|  Figure as above but also displaying the underlying lavaan fits' point estimates and 95%CIs.

# Also see in comparison to lavaan model estimates
forest_plot_data %>% 
  filter(x_var == "Hours") %>% 
  forest_plot(lavaan = TRUE)
```

### Experiences <-> WB

```{r forest-motivations-wb-1}
#| fig.height: 5
#| fig.cap: >
#|  Forest plots of meta-analytic raw cross-lagged regression coefficients from models examining motivations and well-being (motivations to well-being). Shaded curves indicate approximate posterior densities, which are numerically summarised by means and 95%CIs in the right margin. Numbers below average effects indicate posterior probabilities of direction.

a <- forest_plot_data %>% 
  filter(x_var %in% c("Intrinsic"), Parameter == "wy2 ~ wx1") %>% 
  forest_plot() +
  theme(axis.title.x = element_blank())

b <- forest_plot_data %>% 
  filter(x_var %in% c("Extrinsic"), Parameter == "wy2 ~ wx1") %>% 
  forest_plot()

a/b
```


```{r, include = FALSE}
agg_png(
  "Figures/Fig-forest-exp-wb.png", 
  width = W, height = W*0.6, units = "cm", res = 400
)
last_plot()
dev.off()
```

### Some post-hoc calculations

The model parameters report effects of one-unit (e.g. one hour per day) changes in the predictor on the outcome. We also expand on this by first asking what is the effect of an average deviation in a predictor on the outcome. That is, for example, we calculate the average range of players' average daily hours played, and then answer what the effect of an average player moving from their smallest to greatest amount of play would be on well-being.

```{r average-range-effects}
# Average ranges of predictors
avg_ranges <- d %>%
  # Calculate range for everyone (returns -Inf if all values missing)
  group_by(Game, pid) %>%
  summarise(
    across(
      c(Intrinsic, Extrinsic, Hours, hours_est),
      .fns = ~max(.x, na.rm = TRUE) - min(.x, na.rm = TRUE)
    )
  ) %>% 
  ungroup() %>% 
  # Calculate grand mean range excluding -Inf values
  summarise(
    across(
      Intrinsic:hours_est, 
      ~mean(if_else(.x==-Inf, NaN, .x), na.rm = TRUE)
    )
  ) %>% 
  pivot_longer(
    everything(),
    names_to = "x_var",
    values_to = "Mean_range"
  )
fit_ma %>%
  filter(
    model == "AR",
    Type == "Unstandardized",
    Parameter == "wy2 ~ wx1",
    x_var != "hours_est",
    y_var %in% c("Affect", "Life satisfaction")
  ) %>%
  left_join(avg_ranges) %>% 
  mutate(
    hypothesis = map2(
      posterior,
      Mean_range,
      ~ hypothesis(.x, str_glue("Average * {.y} = 0")) %>% .[[1]]
    )
  ) %>%
  unnest(hypothesis) %>%
  select(x_var, y_var, Mean_range, Parameter, Estimate:CI.Upper) %>%
  kbl(digits = 3, caption = "ES of average hours per week on WB") %>% 
  kable_classic(full_width = FALSE, html_font = Font)
```

We then turn to a second question, how much would the average player need to increase their play to result in a subjectively noticeable effect. Anvari and Lakens estimated that a 2% movement might be a lower limit to a subjectively noticeable change in well-being using a similar measurement scale. 

We calculate, for each predictor and outcome, how much the predictor would need to change in order to elicit a change in the outcome greater than the subjective threshold. e.g. How many more hours would one need to play to "feel it".

```{r subjective-threshold-estimates}
# Subjective threshold for each outcome
subjective_thresholds <- tibble(
  y_var = c("Affect", "Life satisfaction"),
  Threshold = 0.02 * c(13, 11)
)

# Table of effects and thresholds
fit_ma %>%
  filter(
    model == "AR",
    Type == "Unstandardized",
    Parameter == "wy2 ~ wx1",
    x_var != "hours_est",
    y_var %in% c("Affect", "Life satisfaction")
  ) %>% 
  select(x_var, y_var, summary) %>% 
  arrange(x_var, y_var) %>% 
  unnest(summary) %>% 
  select(x_var, y_var, Game, Effect = Mean) %>% 
  filter(Game == "Average") %>% 
  left_join(subjective_thresholds) %>% 
  mutate(`Change needed` = Threshold / abs(Effect)) %>% 
  kbl(
    digits = 2, 
    caption = "Estimated required change in X to elicit a subjectively noticeable change in well-being outcomes"
  ) %>% 
  kable_classic(full_width = FALSE, html_font = Font)
```

How much of posterior density is within ROPE defined by A&L subjectively noticeable limits?

```{r subjective-threshold-rope}
# Table of percentages of effects' posterior distributions with ROPE (as defined by thresholds for being subjectively noticed)
fit_ma %>%
  filter(
    model == "AR",
    Type == "Unstandardized",
    Parameter == "wy2 ~ wx1",
    x_var != "hours_est",
    y_var %in% c("Affect", "Life satisfaction")
  ) %>%
  left_join(subjective_thresholds) %>% 
  mutate(
    out = map2(
      posterior,
      Threshold,
      ~ describe_posterior(.x, rope_range = c(-.y, .y), rope_ci = 1) %>% 
        rename(Game = Parameter) %>% 
        filter(Game == "Average")
    )
  ) %>% 
  select(x_var, y_var, out, Threshold) %>% 
  unnest(out) %>% 
  select(x_var, y_var, Threshold, Game, CI_low, CI_high, ROPE_Percentage) %>%
  mutate(ROPE_Percentage = percent(ROPE_Percentage, .1)) %>%
  kbl(
    digits = 2,
    caption = "ROPE percentages below threshold"
  ) %>% 
  kable_classic(full_width = FALSE, html_font = Font)
```

## Subjective time

Comparing subjective hours estimates to objective play behaviour is complicated by the fact that a missing value in objective behaviour indicates zero hours of play, whereas a missing value in estimated time can indicate either an estimated lack of play, or a missing response (e.g. player declined to respond to this item). Therefore any comparison here is done only on individuals who provided subjective estimate data.

Histograms of hours played (behaviour: top, strong; estimated: bottom, light). Only person-waves with both variables present are included.

```{r time-subjective-objective-histogram}
#| fig.height: 7
#| fig.cap: >
#|  Bivariate histograms of average hours played per day by wave (columns) and game (rows). Histograms above zero, with strong colour, reflect behavioural telemetry log data. Histograms below zero, in light colour, indicate subjective estimates of play. Histogram heights are normalized counts. Small triangles indicate means.

tmp <- d %>%
  drop_na(Hours, hours_est) %>%
  group_by(Game, pid, Wave) %>%
  summarise(
    Hours = sum(Hours, na.rm = T),
    hours_est = sum(hours_est, na.rm = T)
  ) %>%
  group_by(Game, Wave) %>%
  summarise(
    across(c(Hours, hours_est), .fns = list(m = mean, se = ~ sd(.x) / sqrt(length(.x))))
  )

d %>%
  select(Game, pid, Wave, Hours, hours_est) %>%
  drop_na(Hours, hours_est) %>%
  filter(Hours < 3, hours_est < 3) %>%
  ggplot(aes()) +
  scale_y_continuous(
    "Normalised count",
    breaks = pretty_breaks()
  ) +
  scale_x_continuous(breaks = pretty_breaks()) +
  geom_histogram(
    aes(x = Hours, y = stat(ncount)),
    bins = 30, fill = col1, col = "white", alpha = .75
  ) +
  geom_histogram(
    aes(x = hours_est, y = stat(ncount) * -1),
    bins = 30, fill = col1, col = "white", alpha = .5
  ) +
  geom_pointrangeh(
    data = tmp, shape = 25, fill = col1, alpha = 1,
    aes(
      x = Hours_m,
      xmin = Hours_m - Hours_se,
      xmax = Hours_m + Hours_se,
      y = .075,
    )
  ) +
  geom_pointrangeh(
    data = tmp, shape = 24, fill = col1, alpha = .75,
    aes(
      x = hours_est_m,
      xmin = hours_est_m - hours_est_se,
      xmax = hours_est_m + hours_est_se,
      y = -.075
    )
  ) +
  facet_grid(Game ~ Wave, scales = "free")
```

A scatterplot of subjective on objective time. Only person-waves with both variables present are included.

```{r time-subjective-objective-scatterplot}
#| fig.height: 4
#| fig.cap: >
#|  Scatterplots of subjectively estimated average hours of daily play on objective behavioural hours played. Dashed line indicates identity, solid lines are simple regression lines. Small blue points are individual participants, and solid dark points are bivariate sample means.

d %>%
  drop_na(Hours, hours_est) %>%
  filter(Hours < 3, hours_est < 3) %>%
  select(Game, pid, Wave, Objective = Hours, Subjective = hours_est) %>%
  ggplot(aes(Objective, Subjective)) +
  scale_x_continuous(
    "Hours played per day (behaviour)",
    breaks = c(0, 1, 2, 3), labels = c(0, 1, 2, 3)
  ) +
  scale_y_continuous(
    "Hours played per day (estimate)",
    breaks = pretty_breaks()
  ) +
  geom_point(size = .2, alpha = .2, shape = 1, col = col1) +
  geom_abline(intercept = 0, slope = 1, lty = 2, size = .25) +
  stat_centroid() +
  geom_smooth(
    method = "lm", size = .4,
    color = "black",
    alpha = .33, show.legend = FALSE
  ) +
  facet_grid(Wave ~ Game) +
  theme(aspect.ratio = 1)
```

Wave-person differences. Difference indicates estimate - behaviour, so positive numbers indicate overestimates.

```{r time-subjective-objective-differences}
#| fig.height: 7
#| fig.cap: >
#|  Density curves of differences between subjectively estimated and objectively logged average daily play time.

d %>%
  select(Game, pid, Wave, Hours, hours_est) %>%
  mutate(Difference = hours_est - Hours) %>%
  filter(Hours < 3, hours_est < 3) %>%
  drop_na(Difference) %>%
  ggplot(aes(Difference, Game)) +
  scale_x_continuous(
    "Estimation error (hours/day)",
    breaks = pretty_breaks()
  ) +
  coord_cartesian(xlim = c(-0.5, 0.5)) +
  geom_vline(xintercept = 0, lty = 2, size = .25) +
  stat_halfeye(point_interval = NULL) +
  stat_summaryh(fun.data = mean_cl_boot_h) +
  facet_wrap("Wave", ncol = 1)

# Overall confidence interval for estimation error
d %>%
  select(Game, pid, Wave, Hours, hours_est) %>%
  mutate(Difference = hours_est - Hours) %>%
  filter(Hours < 3, hours_est < 3) %>%
  drop_na(Difference) %>%
  group_by(Wave) %>%
  summarise(mean_cl_boot(Difference)) %>%
  kbl(digits = 2, caption = "Estimation errors (Estimated - Behaviour)") %>% 
  kable_classic(full_width = FALSE, html_font = Font)
```

### Meta-analysis

```{r forest-subjective-wb}
fit_ma %>% 
  filter(
    model == "AR",
    Type == "Unstandardized",
    x_var == "Hours",
    y_var %in% c("Affect", "Life satisfaction"),
    Parameter == "wy2 ~ wx1"
  ) %>% 
  select(-posterior) %>% 
  unnest(summary) %>% 
  mutate(
    pd = percent(pd, .1),
    across(where(is.numeric), ~format(round(.x, 2), nsmall = 2)),
    Estimate = str_glue("{Mean} [{CI_low}, {CI_high}]")
  ) %>% 
  select(y_var, Game, Estimate, pd) %>% 
  kbl(caption = "Meta-analytic parameters of subjective time on WB") %>% 
  row_spec(c(8, 16), TRUE) %>% 
  kable_classic(full_width = FALSE, html_font = Font)
```

## Variability

```{r}
fit_ma %>% 
  filter(
    model == "AR",
    Type == "Unstandardized",
    y_var %in% c("Affect", "Life satisfaction"),
    x_var == "Hours",
    Parameter == "wy2 ~ wx1"
  ) %>% 
  select(y_var, sd) %>% 
  unnest(sd) %>% 
  kbl(
    caption = "Between game standard deviations of hours -> WB", digits = 2
  ) %>% 
  kable_classic(full_width = FALSE, html_font = Font)
```

## System information

```{r}
sessionInfo()
```