---
title: "R Notebook"
output: html_notebook
---
```{r}
library(tidyverse)
library(tidybayes)
library(rstan)
library(ggplot2)
library(ggpattern)
library(shinystan)
library(loo)
library(patchwork)
library(HDInterval)

mytheme <-
  theme(
    plot.title = element_text(
      hjust = 0.5,
      face = "bold",
      size = (25),
      colour = "steelblue4"
    ),
    axis.title = element_text(
      size = (25),
      colour = "steelblue4",
      face = "bold"
    ),
    axis.text = element_text(colour = "cornflowerblue", size = (25)),
    axis.title.y = element_text(vjust = 0.5),
    panel.grid = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(colour = "black"),
    strip.text = element_text(size = 25),
    legend.title = element_text(size = 25),
    legend.text = element_text(size=25),
    legend.key.size = unit(2, 'cm')
  )

cb1 <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

log_to_prob <- function(logit){
  return(1 / (1 + exp(-logit)))
}

prob_to_log <- function(p){
  return(log(p / (1 - p)))
}

softmax_function <- function(values) {
  exp_values <- exp(values)
  p_values <- exp_values / sum(exp_values)
  return(p_values)
}
```

## Load in the datasets 

```{r}
Player_two <- read.csv("Player_two.csv")
```

## A little extra tidying 

Perform a little bit of extra tidying for Player two. 

```{r}
Player_two <- Player_two %>%
  mutate(Answer_difference = (Second_answer - First_answer),
         Information_difference = Social_info - First_answer) %>%
  filter(Information_difference != 0) %>%
  mutate(Condition_recoded = case_when(
           Condition == "Cooperative" ~ 1,
           Condition == "Competitive" ~ 2
         ),
         Condition_recoded = factor(Condition_recoded),
         Spied = ifelse(Type == "Spied", 1, 0),
         Advice = ifelse(Type == "Advice", 1, 0),
         participant_id = match(participant_id, unique(participant_id)),
         participant_id = factor(participant_id)) 

table(Player_two$Spied, Player_two$Type)
table(Player_two$Advice, Player_two$Type)
```

Divide by 20 to get the condition number, as there are 20 rows of data per participant. 

```{r}
table(Player_two$Condition) / 20
```

## Player two descriptives

Participant's use of social information - their "shift" strategy

```{r}
Shift <- Player_two %>%
  mutate(Shift = (Second_answer - First_answer) / (Social_info - First_answer)) %>%
  filter(Shift < 1.5) %>%
  filter(Shift > -0.5) %>%
  select(Shift, Condition)

Shift %>%
  mutate(Condition = factor(Condition, levels = c("Cooperative", "Competitive"))) %>%
  ggplot(aes(x = Shift)) + 
  geom_histogram(bins = 25, fill = "#999999", col = "black") + 
  labs(x = "Proportion shift") + 
  mytheme + 
  theme(axis.title.y = element_blank()) + 
  facet_grid(~Condition)

ggsave("Information_shift.png", width = 12, height = 6, dpi = 300)
```

### Model two

Reminder of the conditions:

1 = Cooperative
2 = Competitive
```{r}
Shift_by_info <- Player_two %>%
  mutate(Shift = (Second_answer - First_answer) / (Social_info - First_answer)) %>%
  filter(Shift < 1.5, Shift > -0.5) %>%
  filter(Type %in% c("Advice", "Spied")) %>%
  mutate(`Information type` = factor(Type, levels = c("Advice", "Spied"))) %>%
  select(Shift, `Information type`)

ggplot(Shift_by_info, aes(x = Shift)) +
  geom_histogram(bins = 25, fill = "#999999", col = "black") +
  labs(x = "Proportion shift") +
  mytheme +
  theme(axis.title.y = element_blank()) +
  facet_grid(~`Information type`)

ggsave("Information_shift_by_type.png", width = 12, height = 6, dpi = 300)
```
```{r}
Shift_by_info <- Player_two %>%
  mutate(Shift = (Second_answer - First_answer) / (Social_info - First_answer)) %>%
  filter(Shift < 1.5, Shift > -0.5) %>%
  filter(Type %in% c("Advice", "Spied")) %>%
  mutate(`Information type` = factor(Type, levels = c("Advice", "Spied"))) %>%
  select(Shift, `Information type`)

# Extract the two samples
advice_shifts <- Shift_by_info %>% filter(`Information type` == "Advice") %>% pull(Shift)
spied_shifts  <- Shift_by_info %>% filter(`Information type` == "Spied") %>% pull(Shift)

# Generate posterior over difference in means
n_draws <- 10000

samples_advice <- rnorm(n_draws, mean(advice_shifts), sd(advice_shifts) / sqrt(length(advice_shifts)))
samples_spied  <- rnorm(n_draws, mean(spied_shifts),  sd(spied_shifts) / sqrt(length(spied_shifts)))

diff_samples <- samples_spied - samples_advice

# Summarize posterior
tibble(Diff = diff_samples) %>%
  summarise(
    `Mean Diff` = mean(Diff),
    `95% CI Lower` = quantile(Diff, 0.025),
    `95% CI Upper` = quantile(Diff, 0.975),
    `P(Spied > Advice)` = mean(Diff > 0)
  )
```

```{r}
B_stan <- Player_two %>%
  compose_data(y = Answer_difference,
               Info_type = Type,
               Condition = Condition_recoded,
               ppt = participant_id,
               Advice = Advice,
               Spied = Spied)

table("Type" = B_stan$Type, "Advice" = B_stan$Advice)
```

Non centered parameterisation was necessary for this model to run.

```{stan output.var="B_Model_1"}
data {
  int<lower=1> n;
  vector[n] Answer_difference; 
  vector[n] Information_difference;
  int n_participant_id; 
  int Advice[n];
  int Condition[n];
  int ppt[n];
}
parameters {
  real<lower=0,upper=5> sigma; 
  real<lower=0,upper=5> sigma_participant;
  vector[2] b_condition;
  real b_advice;
  vector[n_participant_id] participant_effect;
}
transformed parameters {
  vector[n] influence;
  vector[n] mu; 
  for(i in 1:n){
   influence[i] = b_condition[Condition[i]] * (1 + b_advice * Advice[i]) + participant_effect[ppt[i]];
   mu[i] =  influence[i] * Information_difference[i];
  }  
}
model {
  Answer_difference ~ normal(mu, sigma);
  sigma ~ exponential(1); 
  b_advice ~ normal(0,0.5);
  b_condition ~ normal(0,0.5);
  participant_effect ~ normal(0, sigma_participant);
  sigma_participant ~ exponential(1);
}
```

```{r}
ModelB_1 <- sampling(B_Model_1, iter = 5000, warmup = 1500, chains = 4, cores = 4, data = B_stan)
ModelB_1 

# summary <- summary(ModelB_1)
# rhat <- summary$summary[, "Rhat"]
# max(rhat)

#launch_shinystan(ModelB_1)
```

### Generating predictions

```{r}
Model_1B <- extract(ModelB_1)

# Conditions
Cooperative_condition <- Model_1B$b_condition[,1]
Competitive_condition <- Model_1B$b_condition[,2]

# Predictors
Advice <- Model_1B$b_advice
```

Showing parameter intervals 

```{r}
mean_qi(Cooperative_condition)
mean_qi(Competitive_condition)

mean_qi(Advice)
```

Contrasts between parameters 

```{r}
mean_qi(Competitive_condition - Cooperative_condition)
```

Generating predictions

```{r}
Influence_Predictions <- tibble(
  Cooperative_Spied = Cooperative_condition,
  Cooperative_Advice = Cooperative_condition * (1 + Advice),
  
  Competitive_Spied = Competitive_condition,
  Competitive_Advice = Competitive_condition * (1 + Advice)
)

Plotting_influence_predictions <- Influence_Predictions %>%
  pivot_longer(everything(), names_to = "Combination", values_to = "Estimate") %>%
  separate(Combination, into = c("Condition", "Information type"))
```

Putting them into a plot

```{r}
Plotting_influence_predictions %>%
  mutate(Condition = factor(Condition, levels = c("Cooperative", "Competitive"))) %>%
  ggplot(aes(x = Condition, y = Estimate, fill = `Information type`)) + 
  stat_eye(position = "dodge") +
  scale_fill_manual(values = cb1) +
  coord_cartesian(ylim = c(-0.05, 1)) +  # ← More breathing room on both ends!
  mytheme + 
  labs(y = "Social influence") + 
  theme(axis.title.x = element_blank(),
        panel.spacing = unit(1, "lines"),
        axis.text.x = element_text(size = 17))

ggsave("Influence_plot.png", width = 12, height = 8, dpi = 300)
```

Predicted influence values

```{r}
Plotting_influence_predictions %>%
  group_by(Condition, `Information type`) %>%
  summarise(
    Mean = mean(Estimate),
    `Upper 95%` = HDInterval::hdi(Estimate, credMass = 0.95)[2],
    `Lower 95%` = HDInterval::hdi(Estimate, credMass = 0.95)[1],
    `P(Estimate > 0)` = mean(Estimate > 0),
    .groups = 'drop'
  )

# Extract the samples from your Influence_Predictions tibble
Coop_Advice <- Influence_Predictions$Cooperative_Advice
Comp_Advice <- Influence_Predictions$Competitive_Advice

Coop_Spied <- Influence_Predictions$Cooperative_Spied
Comp_Spied <- Influence_Predictions$Competitive_Spied

# Calculate the probability that Cooperative Advice > Competitive Advice
p_coop_advice_beats_comp_advice <- mean(Coop_Advice > Comp_Advice)

# Show the result
p_coop_advice_beats_comp_advice

p_coop_spied_beats_comp_spied <- mean(Coop_Spied > Comp_Spied)

p_coop_spied_beats_comp_spied
```

```{r}
library(tidyverse)
library(ggpattern)


# Combine all CSVs
combined_df <- bind_rows(
  read_csv("Player_two_gpt-4o_direct.csv") %>% mutate(Model = "GPT-4o (Direct)", Color = "forestgreen", Pattern = "none"),
  read_csv("Player_two_gpt-4o_cot.csv") %>% mutate(Model = "GPT-4o (CoT)", Color = "forestgreen", Pattern = "stripe"),
  read_csv("Player_two_claude-3-5-sonnet-20241022_direct.csv") %>% mutate(Model = "Claude 3.5 Sonnet (Direct)", Color = "firebrick", Pattern = "none"),
  read_csv("Player_two_claude-3-5-sonnet-20241022_cot.csv") %>% mutate(Model = "Claude 3.5 Sonnet (CoT)", Color = "firebrick", Pattern = "stripe"),
  read_csv("Player_two.csv") %>% mutate(Model = "Human (p < 0.05)", Color = "steelblue", Pattern = "none")
)

# Preprocess and compute shift
summary_df <- combined_df %>%
  mutate(Shift = (Second_answer - First_answer) / (Social_info - First_answer)) %>%
  filter(Shift < 1.5, Shift > -0.5) %>%
  filter(Type %in% c("Advice", "Spied")) %>%
  group_by(Model, Type, Color, Pattern) %>%
  summarise(Average_Shift = mean(Shift, na.rm = TRUE), .groups = 'drop')

# Print average shifts
print(summary_df %>% arrange(factor(Model, levels = c(
  "Claude 3.5 Sonnet (CoT)",
  "Claude 3.5 Sonnet (Direct)",
  "GPT-4o (CoT)",
  "GPT-4o (Direct)",
  "Human (p < 0.05)"
))))

# Factor Model to control bar order to match legend
summary_df$Model <- factor(summary_df$Model, levels = c(
  "Claude 3.5 Sonnet (CoT)",
  "Claude 3.5 Sonnet (Direct)",
  "GPT-4o (CoT)",
  "GPT-4o (Direct)",
  "Human (p < 0.05)"
))

# Plot
ggplot(summary_df, aes(x = Type, y = Average_Shift, fill = Model, pattern = Pattern)) +
  geom_bar_pattern(
    stat = "identity",
    position = position_dodge(width = 0.9),
    width = 0.9,
    color = "black",
    pattern_fill = "black",
    pattern_angle = 45,
    pattern_density = 0.1,
    pattern_spacing = 0.02,
    pattern_key_scale_factor = 0.6
  ) +
  scale_fill_manual(
    values = c(
      "GPT-4o (Direct)" = "forestgreen",
      "GPT-4o (CoT)" = "forestgreen",
      "Claude 3.5 Sonnet (Direct)" = "firebrick",
      "Claude 3.5 Sonnet (CoT)" = "firebrick",
      "Human (p < 0.05)" = "steelblue"
    )
  ) +
  scale_pattern_manual(
    values = c("none" = "none", "stripe" = "stripe")
  ) +
  labs(
    x = "Information Type",
    y = "Average Proportion Shift",
    title = "Average Proportion Shift by Information Type and Model",
    fill = "Player 2 Entity"
  ) +
  guides(
    pattern = "none",
    fill = guide_legend(override.aes = list(
      pattern = c("stripe", "none", "stripe", "none", "none")
    ))
  ) +
  theme_minimal(base_size = 16) +
  theme(
    plot.title = element_text(size = 18),
    plot.margin = margin(t = 10, r = 50, b = 10, l = 10)
  )

ggsave("average_info_shift_graph.png", width = 10, height = 6, dpi = 300)
```
```{r}
library(boot)

# Combine all CSVs
combined_df <- bind_rows(
  read_csv("Player_two_gpt-4o_direct.csv") %>% mutate(Model = "GPT-4o (Direct)", Color = "forestgreen", Pattern = "none"),
  read_csv("Player_two_gpt-4o_cot.csv") %>% mutate(Model = "GPT-4o (CoT)", Color = "forestgreen", Pattern = "stripe"),
  read_csv("Player_two_claude-3-5-sonnet-20241022_direct.csv") %>% mutate(Model = "Claude 3.5 Sonnet (Direct)", Color = "firebrick", Pattern = "none"),
  read_csv("Player_two_claude-3-5-sonnet-20241022_cot.csv") %>% mutate(Model = "Claude 3.5 Sonnet (CoT)", Color = "firebrick", Pattern = "stripe"),
)

# Preprocess and compute shift
processed_df <- combined_df %>%
  mutate(Shift = (Second_answer - First_answer) / (Social_info - First_answer)) %>%
  filter(Shift < 1.5, Shift > -0.5) %>%
  filter(Type %in% c("Advice", "Spied"))


# Bootstrapping function
boot_ci <- function(data, indices) {
  mean(data[indices], na.rm = TRUE)
}

# Compute summary with bootstrapped CI
summary_df <- processed_df %>%
  group_by(Model, Type, Color, Pattern) %>%
  summarise(
    Average_Shift = mean(Shift, na.rm = TRUE),
    CI_Low = boot(Shift, boot_ci, R = 1000)$t0 - quantile(boot(Shift, boot_ci, R = 1000)$t, 0.975) + mean(Shift, na.rm = TRUE),
    CI_High = boot(Shift, boot_ci, R = 1000)$t0 - quantile(boot(Shift, boot_ci, R = 1000)$t, 0.025) + mean(Shift, na.rm = TRUE),
    .groups = 'drop'
  )

# Print average shifts
print(summary_df %>% arrange(factor(Model, levels = c(
  "Claude 3.5 Sonnet (CoT)",
  "Claude 3.5 Sonnet (Direct)",
  "GPT-4o (CoT)",
  "GPT-4o (Direct)"
))))

# Reorder for plotting
summary_df$Model <- factor(summary_df$Model, levels = c(
  "Claude 3.5 Sonnet (CoT)",
  "Claude 3.5 Sonnet (Direct)",
  "GPT-4o (CoT)",
  "GPT-4o (Direct)"
))

baseline_means <- processed_df %>%
  group_by(Type) %>%
  summarise(Baseline_Shift = mean(Shift, na.rm = TRUE)) %>%
  deframe() 

cat("Average Proportion Shift (All Models)\n")
cat("--------------------------------------\n")
cat(sprintf("Advice: %.4f\n", baseline_means["Advice"]))
cat(sprintf("Spied : %.4f\n", baseline_means["Spied"]))

baseline_df <- tibble(
  Type = c("Advice", "Spied"),
  y = c(baseline_means["Advice"], baseline_means["Spied"]),
  Label = c("Advice Baseline", "Spied Baseline")
)

# Plot
ggplot(summary_df, aes(x = Type, y = Average_Shift, fill = Model, pattern = Pattern)) +
  geom_bar_pattern(
    stat = "identity",
    position = position_dodge(width = 0.9),
    width = 0.9,
    color = "black",
    pattern_fill = "black",
    pattern_angle = 45,
    pattern_density = 0.1,
    pattern_spacing = 0.02,
    pattern_key_scale_factor = 0.6
  ) +
  geom_errorbar(aes(ymin = CI_Low, ymax = CI_High),
                position = position_dodge(width = 0.9),
                width = 0.2,
                color = "black") +
  geom_segment(
    data = baseline_df,
    aes(x = as.numeric(factor(Type)) - 0.45,
        xend = as.numeric(factor(Type)) + 0.45,
        y = y, yend = y,
        linetype = Label),
    inherit.aes = FALSE,
    color = "black",
    linewidth = 0.8
  ) +
  scale_fill_manual(
    values = c(
      "GPT-4o (Direct)" = "forestgreen",
      "GPT-4o (CoT)" = "forestgreen",
      "Claude 3.5 Sonnet (Direct)" = "firebrick",
      "Claude 3.5 Sonnet (CoT)" = "firebrick"
    )
  ) +
  scale_pattern_manual(
    values = c("none" = "none", "stripe" = "stripe")
  ) +
  scale_linetype_manual(
  name = "",
  values = c("Advice Baseline" = "dashed", "Spied Baseline" = "dotdash")
) +
  labs(
    x = "Information Type",
    y = "Average Proportion Shift",
    title = "Average Proportion Shift by Information Type and Model",
    fill = "Player 2 Entity"
  ) +
  guides(
    pattern = "none",
    fill = guide_legend(override.aes = list(
      pattern = c("stripe", "none", "stripe", "none")
    ))
  ) +
  theme_minimal(base_size = 16) +
  theme(
    plot.title = element_text(size = 18),
    plot.margin = margin(t = 10, r = 50, b = 10, l = 10)
  )

ggsave("average_info_shift_graph.png", width = 10, height = 6, dpi = 300)
```
```{r}
library(tidyverse)
library(boot)

# Combine all CSVs
combined_df <- bind_rows(
  read_csv("Player_two_gpt-4o_direct.csv") %>% mutate(Model = "GPT-4o (Direct)"),
  read_csv("Player_two_gpt-4o_cot.csv") %>% mutate(Model = "GPT-4o (CoT)"),
  read_csv("Player_two_claude-3-5-sonnet-20241022_direct.csv") %>% mutate(Model = "Claude 3.5 Sonnet (Direct)"),
  read_csv("Player_two_claude-3-5-sonnet-20241022_cot.csv") %>% mutate(Model = "Claude 3.5 Sonnet (CoT)"),
  read_csv("Player_two.csv") %>% mutate(Model = "Humans (p < .05)")
)

# Preprocess and compute shift
processed_df <- combined_df %>%
  mutate(Shift = (Second_answer - First_answer) / (Social_info - First_answer)) %>%
  filter(Shift < 1.5, Shift > -0.5) %>%
  filter(Type %in% c("Advice", "Spied"))

# Bootstrapping function
boot_ci <- function(data, indices) {
  mean(data[indices], na.rm = TRUE)
}

# Compute summary with bootstrapped CI
summary_df <- processed_df %>%
  group_by(Model, Type) %>%
  summarise(
    Average_Shift = mean(Shift, na.rm = TRUE),
    CI_Low = boot(Shift, boot_ci, R = 1000)$t0 - quantile(boot(Shift, boot_ci, R = 1000)$t, 0.975) + mean(Shift, na.rm = TRUE),
    CI_High = boot(Shift, boot_ci, R = 1000)$t0 - quantile(boot(Shift, boot_ci, R = 1000)$t, 0.025) + mean(Shift, na.rm = TRUE),
    .groups = 'drop'
  )

# Print average shifts
print(summary_df %>% arrange(factor(Model, levels = c(
  "Claude 3.5 Sonnet (CoT)",
  "Claude 3.5 Sonnet (Direct)",
  "GPT-4o (CoT)",
  "GPT-4o (Direct)",
  "Humans (p < .05)"
))))

# Reorder for plotting
summary_df$Model <- factor(summary_df$Model, levels = c(
  "Claude 3.5 Sonnet (CoT)",
  "Claude 3.5 Sonnet (Direct)",
  "GPT-4o (CoT)",
  "GPT-4o (Direct)",
  "Humans (p < .05)"
))


# Compute overall advice/spied average
baseline_means <- processed_df %>%
  group_by(Type) %>%
  summarise(Baseline_Shift = mean(Shift, na.rm = TRUE)) %>%
  deframe() 

cat("Average Proportion Shift (All Models)\n")
cat("--------------------------------------\n")
cat(sprintf("Advice: %.4f\n", baseline_means["Advice"]))
cat(sprintf("Spied : %.4f\n", baseline_means["Spied"]))

baseline_df <- tibble(
  Type = c("Advice", "Spied"),
  y = c(baseline_means["Advice"], baseline_means["Spied"]),
  Label = c("Advice Baseline", "Spied Baseline")
)

plot_df <- summary_df %>%
  filter(Model != "Humans (p < .05)")

# Plot
ggplot(plot_df, aes(x = Type, y = Average_Shift, fill = Model)) +
  geom_bar(
    stat = "identity",
    position = position_dodge(width = 0.9),
    width = 0.9,
    color = "black"
  ) +
  geom_errorbar(aes(ymin = CI_Low, ymax = CI_High),
                position = position_dodge(width = 0.9),
                width = 0.2,
                color = "black") +
  geom_segment(
    data = baseline_df,
    aes(x = as.numeric(factor(Type)) - 0.45,
        xend = as.numeric(factor(Type)) + 0.45,
        y = y, yend = y,
        linetype = Label),
    inherit.aes = FALSE,
    color = "black",
    linewidth = 0.8
  ) +
  scale_fill_manual(
  values = c(
    "GPT-4o (Direct)" = "#A6CEE3",
    "GPT-4o (CoT)"    = "#1F78B4",
    "Claude 3.5 Sonnet (Direct)" = "#FDD86E",
    "Claude 3.5 Sonnet (CoT)"    = "#FDBF11",
    "Humans (p < .05)" = "#33A02C"  # green
  )
) +
  scale_linetype_manual(
    name = "",
    values = c("Advice Baseline" = "dashed", "Spied Baseline" = "dotdash")
  ) +
  labs(
    x = "Information Type",
    y = "Average Proportion Shift",
    title = "Average Proportion Shift by Information Type and Model",
    fill = "Player 2 Entity"
  ) +
  coord_cartesian(ylim = c(0, max(summary_df$CI_High) + 0.1)) +
  theme_minimal(base_size = 16) +
  theme(
    plot.title = element_text(size = 18),
    plot.margin = margin(t = 10, r = 50, b = 10, l = 10)
  )

# Save
ggsave("average_info_shift_graph.png", width = 10, height = 4, dpi = 300)

```
```{r}
models <- c("GPT-4o (Direct)", "GPT-4o (CoT)", "Claude 3.5 Sonnet (Direct)", "Claude 3.5 Sonnet (CoT)")

for (m in models) {
  cat("\nModel:", m, "\n")
  
  advice_shifts <- processed_df %>% filter(Model == m, Type == "Advice") %>% pull(Shift)
  spied_shifts  <- processed_df %>% filter(Model == m, Type == "Spied")  %>% pull(Shift)
  
  t_result <- t.test(spied_shifts, advice_shifts, alternative = "greater")
  
  cat(sprintf("  t = %.3f, df = %.1f, p-value = %.4g\n",
              t_result$statistic,
              t_result$parameter,
              t_result$p.value))
}
```
```{r}
mean_shifts <- processed_df %>%
  filter(Model %in% c("GPT-4o (Direct)", "GPT-4o (CoT)", "Humans (p < .05)")) %>%
  group_by(Model, Type) %>%
  summarise(Mean = mean(Shift, na.rm = TRUE), .groups = "drop")

print(mean_shifts)

types <- c("Advice", "Spied")
variants <- c("GPT-4o (Direct)", "GPT-4o (CoT)")

for (variant in variants) {
  for (t in types) {
    gpt_shifts <- processed_df %>% filter(Model == variant, Type == t) %>% pull(Shift)
    human_shifts <- processed_df %>% filter(Model == "Humans (p < .05)", Type == t) %>% pull(Shift)
    
    t_result <- t.test(gpt_shifts, human_shifts)
    
    cat(sprintf("\n%s vs. Humans — %s:\n  t = %.3f, df = %.1f, p = %.4g\n",
                variant, t,
                t_result$statistic,
                t_result$parameter,
                t_result$p.value))
  }
}
```
```{r}
claude_models <- c("Claude 3.5 Sonnet (Direct)", "Claude 3.5 Sonnet (CoT)")

for (variant in claude_models) {
  for (t in c("Advice", "Spied")) {
    claude_shifts <- processed_df %>% filter(Model == variant, Type == t) %>% pull(Shift)
    human_shifts  <- processed_df %>% filter(Model == "Humans (p < .05)", Type == t) %>% pull(Shift)
    
    t_result <- t.test(claude_shifts, human_shifts)
    
    cat(sprintf("\n%s vs. Humans — %s:\n  t = %.3f, df = %.1f, p = %.4g\n",
                variant, t,
                t_result$statistic,
                t_result$parameter,
                t_result$p.value))
  }
}
```
