SIM VS EMPERICAL VS PPC

library(reticulate)
# use_python("/opt/anaconda3/envs/hssm_venv/bin/python")
use_condaenv("hssm_venv", required = TRUE)

import pandas as pd
import numpy as np

source scripts with functions

source("Functions/packages.R")

load packages

loadpackages()
conflicts_prefer(rstatix::filter)
conflicts_prefer(dplyr::rename)
conflicts_prefer(dplyr::summarise)
conflicts_prefer(dplyr::summarize)
conflicts_prefer(dplyr::mutate)
conflicts_prefer(dplyr::arrange)
conflicts_prefer(dplyr::count)
conflicts_prefer(dplyr::filter)

GLOBAL SETTINGS (font, colors, paths)

poster_font <- "Book Antiqua"
if (requireNamespace("systemfonts", quietly = TRUE)) {
  if (!poster_font %in% systemfonts::system_fonts()$family) {
    poster_font <- "Palatino"
    message("Book Antiqua not found on this system -- falling back to Palatino")
  }
}
theme_poster_font <- theme(text = element_text(family = poster_font))
save_device <- if (requireNamespace("ragg", quietly = TRUE)) ragg::agg_png else "png"

# one color per SOURCE (orange = observed/empirical and blue = model-predicted,
# to stay consistent with the script-4 PPC figures)
colors_source <- c(
  "Simulated (power analysis)" = "#7570B3",
  "Empirical (pre-trip)"       = "#FF7F0E",
  "Predicted (GLMM)"           = "#1B9E77",
  "PPC (DDM)"                  = "#1F77B4"
)

# ---- EDIT PATHS HERE IF NEEDED (relative to project root) --------------------
sim_csv       <- "simed_data_02.18.26_mood_congruency.csv"
empirical_csv <- "pre-processed_full-ema_data/full_ema_pre-trip.csv"
ppc_csv       <- "post-model_outputs_cpc2026/ppc_sim_data_full_cpc2026.csv"
# ppc_csv is exported by pasting this at the bottom of script 4 and running once:
#   sim_data_full.to_csv(f"{OUTPUT_PATH}/ppc_sim_data_full_cpc2026.csv", index=False)
#-------------------------------------------------------------------------------

has_ppc <- file.exists(here::here(ppc_csv))
if (!has_ppc) message("PPC csv not found -- PPC rows/plots will be skipped. ",
                      "Export it from script 4 (see comment above).")

if (!dir.exists(here::here("figures"))) dir.create(here::here("figures"))

==============================================================================

LOAD + PREP EACH SOURCE

==============================================================================

1. SIMULATED (power analysis congruency experiment)

sim_raw <- read_csv(here::here(sim_csv))

# the sim df was re-indexed to a full 30 EMAs x 24 trials grid, so subjects
# with <30 EMAs have all-NA rows -> completed trials are rows with non-NA rt
sim_trials <- sim_raw %>% filter(!is.na(rt))

# one row per subject x EMA for session-level (affect) descriptives
sim_sessions <- sim_trials %>%
  distinct(subject_id, ema_num, affect_z_emalevel)

# PosNeg trials + chose_positive (same coding as poweranalyses script)
sim_posneg <- sim_trials %>%
  filter(delta_valence != 0) %>%
  mutate(
    chose_pos = as.integer(
      (delta_valence == -1 & response_0_1 == 0) |
        (delta_valence == 1 & response_0_1 == 1)
    ),
    rt_s = rt,                      # sim rt is already in SECONDS
    affect_z = affect_z_emalevel,
    source = "Simulated (power analysis)"
  )

tibble(
  step = "simulated congruency experiment",
  n_subjects = n_distinct(sim_trials$subject_id),
  n_sessions = nrow(sim_sessions),
  n_trials = nrow(sim_trials),
  n_posneg_trials = nrow(sim_posneg)
) %>% gt() %>% tab_options(table.width = pct(80))

step	n_subjects	n_sessions	n_trials	n_posneg_trials
simulated congruency experiment	70	2010	48240	32160

2. EMPIRICAL (pre-trip)

emas_pre <- read_csv(here::here(empirical_csv))

emas_pre$rt <- as.numeric(emas_pre$rt)
emas_pre$response <- as.numeric(emas_pre$response)
emas_pre$ema_number <- emas_pre$session_number

emas_pre$delta_valence <- ifelse(emas_pre$valence %in% c("PosPos", "NegNeg"), 0,
                          ifelse(emas_pre$right_valence == "positive", 1, -1))

emas_pre$chose_pos  <- ifelse(emas_pre$delta_valence == 0, NA,
                        ifelse(emas_pre$delta_valence == 1,
                               ifelse(emas_pre$response == 1, 1, 0),
                               ifelse(emas_pre$response == 0, 1, 0)))

emas_pre$delta_posneg_sb  <- ifelse(emas_pre$delta_valence == 0, NA,
                              ifelse(emas_pre$delta_valence == 1,
                                     emas_pre$right_base_rating - emas_pre$left_base_rating,
                                     emas_pre$left_base_rating - emas_pre$right_base_rating))

emas_pre$is_noresponse <- emas_pre$valid_response == FALSE | is.na(emas_pre$valid_response)
emas_pre$under_200ms <- emas_pre$rt < 200
emas_pre$somehow_over_5000ms <- emas_pre$rt > 5000

emas_pre$rt_quality <- case_when(
  emas_pre$is_noresponse ~ "no_response",
  emas_pre$under_200ms ~ "under_200ms",
  emas_pre$somehow_over_5000ms ~ "somehow_over_5000ms",
  TRUE ~ "valid"
)

emas_pre  <- emas_pre  %>% filter(rt_quality == "valid", !is.na(response), !is.na(rt))

emas_pre <- emas_pre %>%
         filter(!is.na(affect_z),
         !is.na(rt),
         rt >= 200,
         rt <= 5000)

emas_pre <- emas_pre %>% mutate(rt_s = rt / 1000)   # empirical rt is in MS -> convert to s

# one row per subject x EMA for session-level (affect) descriptives
emp_sessions <- emas_pre %>%
  distinct(PUNS_ID, ema_number, affect, affect_z)

posneg_pre <- emas_pre %>%
  filter(!is.na(chose_pos)) %>%
  mutate(source = "Empirical (pre-trip)")

tibble(
  step = "empirical pre-trip (post-QC)",
  n_subjects = n_distinct(emas_pre$PUNS_ID),
  n_sessions = nrow(emp_sessions),
  n_trials = nrow(emas_pre),
  n_posneg_trials = nrow(posneg_pre)
) %>% gt() %>% tab_options(table.width = pct(80))

step	n_subjects	n_sessions	n_trials	n_posneg_trials
empirical pre-trip (post-QC)	26	519	12172	6149

3. PREDICTED (GLMM) – refit script 3 models, keep fitted values

# models copied unchanged from script 3
choice_model_pretrip <- glmer(chose_pos ~ affect_z + delta_posneg_sb + ema_number +
                                (1 + affect_z + delta_posneg_sb + ema_number || PUNS_ID),
                              data = posneg_pre, family = binomial)

posneg_pre <- posneg_pre %>%
  group_by(PUNS_ID) %>%
  mutate(rt_z = as.numeric(scale(rt_s))) %>%
  ungroup()

poly_matrix_rt <- poly(posneg_pre$delta_posneg_sb, degree = 2)
posneg_pre$delta_posneg_sb_lin  <- poly_matrix_rt[, 1]
posneg_pre$delta_posneg_sb_quad <- poly_matrix_rt[, 2]

rt_model_pretrip <- lmerTest::lmer(
  rt_z ~ affect_z * delta_posneg_sb_lin + delta_posneg_sb_quad + ema_number +
    (1 + affect_z + delta_posneg_sb_lin + delta_posneg_sb_quad || PUNS_ID),
  data = posneg_pre)

# fitted values back on the response scales
# (rt_z -> seconds using the same grand mean/sd back-transform as script 3)
pt_rt_stats <- posneg_pre %>%
  group_by(PUNS_ID) %>%
  summarise(pt_mean = mean(rt_s, na.rm = TRUE),
            pt_sd = sd(rt_s, na.rm = TRUE), .groups = "drop")
grand_mean <- mean(pt_rt_stats$pt_mean)
grand_sd <- mean(pt_rt_stats$pt_sd)

posneg_pre$pred_p_chose_pos_glmm <- fitted(choice_model_pretrip)
posneg_pre$pred_rt_z_glmm <- fitted(rt_model_pretrip)
posneg_pre$pred_rt_s_glmm <- posneg_pre$pred_rt_z_glmm * grand_sd + grand_mean

glmm_pred <- posneg_pre %>%
  transmute(
    subject = as.character(PUNS_ID),
    affect_z,
    rt_s = pred_rt_s_glmm,
    p_chose_pos = pred_p_chose_pos_glmm,
    source = "Predicted (GLMM)"
  )

4. PPC (DDM posterior predictive sims from script 4)

if (has_ppc) {
  ppc_raw <- read_csv(here::here(ppc_csv))

  # pred_rt is in SECONDS; pred_choice = 1 means chose right
  # recover chose_pos on PosNeg trials using delta_valence
  ppc_posneg <- ppc_raw %>%
    filter(delta_valence != 0) %>%
    mutate(
      chose_pos = as.integer(
        (delta_valence == -1 & pred_choice == 0) |
          (delta_valence == 1 & pred_choice == 1)
      ),
      rt_s = abs(pred_rt),   # hssm sims can code rt sign by response boundary
      source = "PPC (DDM)"
    )

  tibble(
    step = "ppc sims (cpc_model, 100 sims per trial)",
    n_subjects = n_distinct(ppc_posneg$subject),
    n_sims = n_distinct(ppc_posneg$sim_num),
    n_posneg_rows = nrow(ppc_posneg)
  ) %>% gt() %>% tab_options(table.width = pct(80))
}

step	n_subjects	n_sims	n_posneg_rows
ppc sims (cpc_model, 100 sims per trial)	25	100	601000

==============================================================================

COMPLIANCE

==============================================================================

Empirical pre-trip compliance

pt_compliance <- emas_pre %>%
  group_by(PUNS_ID) %>%
  summarise(
    n_emas = n_distinct(ema_number),
    compliance = n_emas / 30
  ) %>%
  ungroup() %>%
  arrange(desc(compliance))

pt_compliance %>%
  gt() %>%
  fmt_percent(columns = compliance, decimals = 1) %>%
  tab_header(title = "Pre-trip EMA compliance by participant") %>%
  tab_options(table.width = pct(80))

Pre-trip EMA compliance by participant
PUNS_ID	n_emas	compliance
1033	29	96.7%
1044	29	96.7%
1015	28	93.3%
1086	28	93.3%
1032	25	83.3%
1107	25	83.3%
1021	24	80.0%
1050	24	80.0%
1083	24	80.0%
1068	21	70.0%
1088	21	70.0%
1017	20	66.7%
1060	19	63.3%
1071	19	63.3%
1027	18	60.0%
1067	18	60.0%
1070	18	60.0%
1065	17	56.7%
1024	16	53.3%
1096	16	53.3%
1016	14	46.7%
1035	14	46.7%
1094	14	46.7%
1026	13	43.3%
1080	13	43.3%
1045	12	40.0%

Simulated compliance

sim_compliance <- sim_trials %>%
  group_by(subject_id) %>%
  summarise(n_emas = n_distinct(ema_num),
            compliance = n_emas / 30,
            .groups = "drop")

# by design: 40 pts x 30 EMAs, 20 pts x 28, 10 pts x 25
sim_compliance %>%
  count(n_emas, name = "n_pts") %>%
  mutate(compliance = n_emas / 30) %>%
  gt() %>%
  fmt_percent(columns = compliance, decimals = 1) %>%
  tab_header(title = "Simulated compliance structure (power analysis design)") %>%
  tab_options(table.width = pct(80))

Simulated compliance structure (power analysis design)
n_emas	n_pts	compliance
25	10	83.3%
28	20	93.3%
30	40	100.0%

Compliance: simulated vs. empirical summary

compliance_compare <- bind_rows(
  sim_compliance %>%
    summarise(source = "Simulated (power analysis)",
              n_pts = n(),
              mean_compliance = mean(compliance),
              median_compliance = median(compliance),
              min_compliance = min(compliance),
              max_compliance = max(compliance)),
  pt_compliance %>%
    summarise(source = "Empirical (pre-trip)",
              n_pts = n(),
              mean_compliance = mean(compliance),
              median_compliance = median(compliance),
              min_compliance = min(compliance),
              max_compliance = max(compliance))
)

compliance_compare %>%
  gt() %>%
  fmt_percent(columns = mean_compliance:max_compliance, decimals = 1) %>%
  tab_header(title = "Compliance: simulated (power analysis) vs. empirical") %>%
  tab_options(table.width = pct(80))

Compliance: simulated (power analysis) vs. empirical
source	n_pts	mean_compliance	median_compliance	min_compliance	max_compliance
Simulated (power analysis)	70	95.7%	100.0%	83.3%	100.0%
Empirical (pre-trip)	26	66.5%	63.3%	40.0%	96.7%

==============================================================================

AFFECT + AFFECT_Z (session level)

==============================================================================

Raw affect (empirical) vs. generating values used in the sims

# empirical: person means / between-person SD of means / mean within-person SD
emp_affect_structure <- emp_sessions %>%
  group_by(PUNS_ID) %>%
  summarise(pt_mean_affect = mean(affect, na.rm = TRUE),
            pt_sd_affect = sd(affect, na.rm = TRUE),
            .groups = "drop") %>%
  summarise(
    mean_of_pt_means = mean(pt_mean_affect, na.rm = TRUE),
    between_pt_sd = sd(pt_mean_affect, na.rm = TRUE),
    mean_within_pt_sd = mean(pt_sd_affect, na.rm = TRUE)
  )

affect_structure_compare <- bind_rows(
  tibble(source = "Simulated (generating values, from TREAD)",
         mean_of_pt_means = 10.675,
         between_pt_sd = 1.442,
         mean_within_pt_sd = 1.473),
  emp_affect_structure %>% mutate(source = "Empirical (pre-trip)")
) %>%
  select(source, everything())

affect_structure_compare %>%
  gt() %>%
  fmt_number(columns = -source, decimals = 3) %>%
  tab_header(title = "Raw affect structure: sim generating values vs. empirical") %>%
  tab_options(table.width = pct(80))

Raw affect structure: sim generating values vs. empirical
source	mean_of_pt_means	between_pt_sd	mean_within_pt_sd
Simulated (generating values, from TREAD)	10.675	1.442	1.473
Empirical (pre-trip)	6.871	1.047	1.544

affect_z descriptives (simulated vs. empirical)

affect_z_compare <- bind_rows(
  sim_sessions %>%
    transmute(source = "Simulated (power analysis)", affect_z = affect_z_emalevel),
  emp_sessions %>%
    transmute(source = "Empirical (pre-trip)", affect_z)
)

affect_z_compare %>%
  group_by(source) %>%
  summarise(
    n_sessions = n(),
    mean = mean(affect_z, na.rm = TRUE),
    sd = sd(affect_z, na.rm = TRUE),
    median = median(affect_z, na.rm = TRUE),
    min = min(affect_z, na.rm = TRUE),
    max = max(affect_z, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  gt() %>%
  fmt_number(columns = mean:max, decimals = 3) %>%
  tab_header(title = "affect_z (session level): simulated vs. empirical") %>%
  tab_options(table.width = pct(80))

affect_z (session level): simulated vs. empirical
source	n_sessions	mean	sd	median	min	max
Empirical (pre-trip)	519	0.000	1.000	0.107	−3.806	2.542
Simulated (power analysis)	2010	0.000	1.000	−0.003	−3.150	3.417

affect_z distributions

p_affect_z <- ggplot(affect_z_compare, aes(x = affect_z, fill = source, color = source)) +
  geom_density(alpha = 0.3, linewidth = 1) +
  scale_fill_manual(values = colors_source, name = "Source") +
  scale_color_manual(values = colors_source, name = "Source") +
  labs(x = "affect_z (within-person z-scored affect)", y = "Density",
       title = "Session-level affect_z: simulated vs. empirical") +
  theme_minimal(base_size = 16) +
  theme(panel.grid = element_blank(),
        panel.border = element_rect(colour = "black", fill = NA, linewidth = 1),
        plot.title = element_text(hjust = 0.5)) +
  theme_poster_font
p_affect_z

ggsave(here::here("figures", "compare_affect_z_sim_vs_empirical.png"),
       plot = p_affect_z, width = 9, height = 6, dpi = 300, bg = "white",
       device = save_device)

Raw affect distribution + trajectory over sessions (empirical)

p_affect_raw <- ggplot(emp_sessions, aes(x = affect)) +
  geom_histogram(bins = 25, fill = colors_source[["Empirical (pre-trip)"]],
                 color = "white", alpha = 0.85) +
  labs(x = "Affect (raw, session level)", y = "Count",
       title = "Empirical raw affect (pre-trip sessions)") +
  theme_minimal(base_size = 16) +
  theme(panel.grid = element_blank(),
        panel.border = element_rect(colour = "black", fill = NA, linewidth = 1),
        plot.title = element_text(hjust = 0.5)) +
  theme_poster_font

affect_by_session <- emp_sessions %>%
  group_by(ema_number) %>%
  summarise(mean_affect = mean(affect, na.rm = TRUE),
            se = sd(affect, na.rm = TRUE) / sqrt(n()),
            .groups = "drop")

p_affect_time <- ggplot(affect_by_session, aes(x = ema_number, y = mean_affect)) +
  geom_ribbon(aes(ymin = mean_affect - se, ymax = mean_affect + se),
              alpha = 0.25, fill = colors_source[["Empirical (pre-trip)"]]) +
  geom_line(color = colors_source[["Empirical (pre-trip)"]], linewidth = 1) +
  geom_point(color = colors_source[["Empirical (pre-trip)"]], size = 2) +
  labs(x = "EMA number (pre-trip)", y = "Mean affect (raw)",
       title = "Empirical raw affect across pre-trip sessions") +
  theme_minimal(base_size = 16) +
  theme(panel.grid = element_blank(),
        panel.border = element_rect(colour = "black", fill = NA, linewidth = 1),
        plot.title = element_text(hjust = 0.5)) +
  theme_poster_font

p_affect_raw + p_affect_time

==============================================================================

RT (PosNeg trials, seconds, all four sources)

==============================================================================

rt_all_sources <- bind_rows(
  sim_posneg %>% select(source, rt_s, affect_z, chose_pos),
  posneg_pre %>% select(source, rt_s, affect_z, chose_pos),
  glmm_pred %>% transmute(source, rt_s, affect_z, chose_pos = NA_integer_),
  if (has_ppc) ppc_posneg %>% select(source, rt_s, affect_z, chose_pos) else NULL
) %>%
  mutate(source = factor(source, levels = names(colors_source)))

rt_all_sources %>%
  group_by(source) %>%
  summarise(
    n_rows = n(),
    mean_rt_s = mean(rt_s, na.rm = TRUE),
    sd_rt_s = sd(rt_s, na.rm = TRUE),
    median_rt_s = median(rt_s, na.rm = TRUE),
    min_rt_s = min(rt_s, na.rm = TRUE),
    max_rt_s = max(rt_s, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  gt() %>%
  fmt_number(columns = mean_rt_s:max_rt_s, decimals = 3) %>%
  tab_header(title = "RT (s), PosNeg trials: simulated vs. empirical vs. GLMM-predicted vs. PPC",
             subtitle = "GLMM row = fitted values (conditional means), so its SD is expectedly smaller than raw/simulated RTs; PPC rows pool 100 simulations per trial") %>%
  tab_options(table.width = pct(90))

RT (s), PosNeg trials: simulated vs. empirical vs. GLMM-predicted vs. PPC
GLMM row = fitted values (conditional means), so its SD is expectedly smaller than raw/simulated RTs; PPC rows pool 100 simulations per trial
source	n_rows	mean_rt_s	sd_rt_s	median_rt_s	min_rt_s	max_rt_s
Simulated (power analysis)	32160	1.500	0.819	1.247	0.591	11.440
Empirical (pre-trip)	6149	1.724	0.819	1.510	0.221	4.982
Predicted (GLMM)	6149	1.725	0.167	1.728	1.117	2.296
PPC (DDM)	601000	1.744	1.081	1.398	0.585	16.892

RT distributions

p_rt_dens <- ggplot(rt_all_sources, aes(x = rt_s, fill = source, color = source)) +
  geom_density(alpha = 0.25, linewidth = 1) +
  coord_cartesian(xlim = c(0, 5)) +
  scale_fill_manual(values = colors_source, name = "Source") +
  scale_color_manual(values = colors_source, name = "Source") +
  labs(x = "Reaction Time (s)", y = "Density",
       title = "RT distributions by source (PosNeg trials)") +
  theme_minimal(base_size = 16) +
  theme(panel.grid = element_blank(),
        panel.border = element_rect(colour = "black", fill = NA, linewidth = 1),
        plot.title = element_text(hjust = 0.5)) +
  theme_poster_font
p_rt_dens

ggsave(here::here("figures", "compare_rt_distributions_by_source.png"),
       plot = p_rt_dens, width = 9, height = 6, dpi = 300, bg = "white",
       device = save_device)

p_rt_box <- ggplot(rt_all_sources, aes(x = source, y = rt_s, fill = source)) +
  geom_boxplot(outlier.alpha = 0.1, width = 0.6) +
  coord_cartesian(ylim = c(0, 5)) +
  scale_fill_manual(values = colors_source, name = "Source") +
  labs(x = NULL, y = "Reaction Time (s)",
       title = "RT by source (PosNeg trials)") +
  theme_minimal(base_size = 16) +
  theme(panel.grid = element_blank(),
        panel.border = element_rect(colour = "black", fill = NA, linewidth = 1),
        plot.title = element_text(hjust = 0.5),
        legend.position = "none",
        axis.text.x = element_text(angle = 20, hjust = 1)) +
  theme_poster_font
p_rt_box

==============================================================================

CHOICE: p(chose positive)

==============================================================================

Overall p(chose pos) by source

p_chose_pos_by_source <- bind_rows(
  sim_posneg %>%
    summarise(source = "Simulated (power analysis)",
              p_chose_pos = mean(chose_pos, na.rm = TRUE),
              n = n()),
  posneg_pre %>%
    summarise(source = "Empirical (pre-trip)",
              p_chose_pos = mean(chose_pos, na.rm = TRUE),
              n = n()),
  glmm_pred %>%
    summarise(source = "Predicted (GLMM)",
              p_chose_pos = mean(p_chose_pos, na.rm = TRUE),
              n = n()),
  if (has_ppc) ppc_posneg %>%
    summarise(source = "PPC (DDM)",
              p_chose_pos = mean(chose_pos, na.rm = TRUE),
              n = n()) else NULL
)

p_chose_pos_by_source %>%
  gt() %>%
  fmt_number(columns = p_chose_pos, decimals = 3) %>%
  tab_header(title = "Overall p(chose positive) on PosNeg trials, by source") %>%
  tab_options(table.width = pct(80))

Overall p(chose positive) on PosNeg trials, by source
source	p_chose_pos	n
Simulated (power analysis)	0.574	32160
Empirical (pre-trip)	0.804	6149
Predicted (GLMM)	0.804	6149
PPC (DDM)	0.801	601000

p(chose pos) by affect bin, by source

# same +/- 0.6 SD binning as script 3
bin_affect <- function(df) {
  df %>%
    mutate(affect_z_bin = case_when(
      affect_z < -0.6 ~ "unpleasant",
      affect_z >  0.6 ~ "pleasant",
      TRUE ~ "neutral"
    ),
    affect_z_bin = factor(affect_z_bin, levels = c("unpleasant", "neutral", "pleasant")))
}

choice_by_affect_source <- bind_rows(
  sim_posneg %>% bin_affect() %>%
    group_by(affect_z_bin) %>%
    summarise(p = mean(chose_pos, na.rm = TRUE),
              se = sd(chose_pos, na.rm = TRUE) / sqrt(n()), .groups = "drop") %>%
    mutate(source = "Simulated (power analysis)"),
  posneg_pre %>% bin_affect() %>%
    group_by(affect_z_bin) %>%
    summarise(p = mean(chose_pos, na.rm = TRUE),
              se = sd(chose_pos, na.rm = TRUE) / sqrt(n()), .groups = "drop") %>%
    mutate(source = "Empirical (pre-trip)"),
  glmm_pred %>% bin_affect() %>%
    group_by(affect_z_bin) %>%
    summarise(p = mean(p_chose_pos, na.rm = TRUE),
              se = sd(p_chose_pos, na.rm = TRUE) / sqrt(n()), .groups = "drop") %>%
    mutate(source = "Predicted (GLMM)"),
  if (has_ppc) ppc_posneg %>% bin_affect() %>%
    group_by(affect_z_bin) %>%
    summarise(p = mean(chose_pos, na.rm = TRUE),
              se = sd(chose_pos, na.rm = TRUE) / sqrt(n()), .groups = "drop") %>%
    mutate(source = "PPC (DDM)") else NULL
) %>%
  mutate(source = factor(source, levels = names(colors_source)))

p_choice_affect <- ggplot(choice_by_affect_source,
       aes(x = affect_z_bin, y = p, color = source, group = source)) +
  geom_point(size = 3, position = position_dodge(width = 0.3)) +
  geom_line(linewidth = 1, position = position_dodge(width = 0.3)) +
  geom_errorbar(aes(ymin = p - se, ymax = p + se), width = 0.15,
                position = position_dodge(width = 0.3)) +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "gray40") +
  scale_color_manual(values = colors_source, name = "Source") +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, by = 0.2)) +
  labs(x = "Affect bin (±0.6 SD)", y = "P(Choose Positive)",
       title = "p(chose positive) by affect, by source (PosNeg trials)") +
  theme_minimal(base_size = 16) +
  theme(panel.grid = element_blank(),
        panel.border = element_rect(colour = "black", fill = NA, linewidth = 1),
        plot.title = element_text(hjust = 0.5)) +
  theme_poster_font
p_choice_affect

ggsave(here::here("figures", "compare_choice_by_affect_by_source.png"),
       plot = p_choice_affect, width = 9, height = 6, dpi = 300, bg = "white",
       device = save_device)

==============================================================================

GRAND SUMMARY TABLE (one row per source)

==============================================================================

grand_summary <- bind_rows(
  sim_posneg %>%
    summarise(source = "Simulated (power analysis)",
              n_pts = n_distinct(subject_id),
              n_sessions = nrow(sim_sessions),
              mean_compliance = mean(sim_compliance$compliance),
              mean_affect_z = mean(affect_z, na.rm = TRUE),
              sd_affect_z = sd(affect_z, na.rm = TRUE),
              mean_rt_s = mean(rt_s, na.rm = TRUE),
              sd_rt_s = sd(rt_s, na.rm = TRUE),
              p_chose_pos = mean(chose_pos, na.rm = TRUE)),
  posneg_pre %>%
    summarise(source = "Empirical (pre-trip)",
              n_pts = n_distinct(PUNS_ID),
              n_sessions = nrow(emp_sessions),
              mean_compliance = mean(pt_compliance$compliance),
              mean_affect_z = mean(affect_z, na.rm = TRUE),
              sd_affect_z = sd(affect_z, na.rm = TRUE),
              mean_rt_s = mean(rt_s, na.rm = TRUE),
              sd_rt_s = sd(rt_s, na.rm = TRUE),
              p_chose_pos = mean(chose_pos, na.rm = TRUE)),
  glmm_pred %>%
    summarise(source = "Predicted (GLMM)",
              n_pts = n_distinct(subject),
              n_sessions = NA_integer_,
              mean_compliance = NA_real_,
              mean_affect_z = NA_real_,
              sd_affect_z = NA_real_,
              mean_rt_s = mean(rt_s, na.rm = TRUE),
              sd_rt_s = sd(rt_s, na.rm = TRUE),
              p_chose_pos = mean(p_chose_pos, na.rm = TRUE)),
  if (has_ppc) ppc_posneg %>%
    summarise(source = "PPC (DDM)",
              n_pts = n_distinct(subject),
              n_sessions = NA_integer_,
              mean_compliance = NA_real_,
              mean_affect_z = NA_real_,
              sd_affect_z = NA_real_,
              mean_rt_s = mean(rt_s, na.rm = TRUE),
              sd_rt_s = sd(rt_s, na.rm = TRUE),
              p_chose_pos = mean(chose_pos, na.rm = TRUE)) else NULL
)

grand_summary %>%
  gt() %>%
  fmt_percent(columns = mean_compliance, decimals = 1) %>%
  fmt_number(columns = c(mean_affect_z, sd_affect_z, mean_rt_s, sd_rt_s, p_chose_pos),
             decimals = 3) %>%
  sub_missing(missing_text = "--") %>%
  tab_header(title = "Grand summary by source (PosNeg trials)",
             subtitle = "affect_z / compliance blank for GLMM + PPC rows since affect enters those as a predictor") %>%
  tab_options(table.width = pct(95))

Grand summary by source (PosNeg trials)
affect_z / compliance blank for GLMM + PPC rows since affect enters those as a predictor
source	n_pts	n_sessions	mean_compliance	mean_affect_z	sd_affect_z	mean_rt_s	sd_rt_s	p_chose_pos
Simulated (power analysis)	70	2010	95.7%	0.000	1.000	1.500	0.819	0.574
Empirical (pre-trip)	26	519	66.5%	0.000	0.998	1.724	0.819	0.804
Predicted (GLMM)	26	–	–	–	–	1.725	0.167	0.804
PPC (DDM)	25	–	–	–	–	1.744	1.081	0.801