```{r setup, include=FALSE}
options(repos = c(CRAN = "https://cran.r-project.org")) 
if (!requireNamespace("groundhog", quietly = TRUE)) {
    install.packages("groundhog")
}
pkgs <- c("knitr", "stringr", "dplyr", "ggplot2", 
          "metafor", "readr", "forcats", "cowplot", "patchwork", 
          "kableExtra", "broom")
groundhog::groundhog.library(pkg = pkgs,
                           date = "2024-04-23")
rm(pkgs)

knitr::opts_chunk$set(
  echo = FALSE,
  warning = FALSE,
  message = FALSE,
  fig.width = 6,
  fig.height = 4,
  fig.align = "center"
)

# Load required packages
# library(dplyr)
# library(ggplot2)
# library(metafor)
# library(readr)
# library(forcats)
# library(cowplot)
# library(patchwork)
# library(kableExtra)
# library(knitr)
# library(broom)
# library(stringr)

# Set working directory and load data
setwd("/Users/meikezehnle/Library/CloudStorage/OneDrive-UniversitaetSt.Gallen/meta")
df <- read.csv("meta-analysis_data.csv", header = TRUE, sep = ",")

# Set output options
options(scipen = 999)
options(digits = 3)

# Compute effect size variances
df$variance_d <- ifelse(df$within_subject_design == "yes",
                   ((2 * (1 - 0.5)) / df$n) + ((df$cohens_d^2) / (2 * df$n)),
                   ((df$n1 + df$n2) / (df$n1 * df$n2)) + (df$cohens_d^2 / (2 * (df$n1 + df$n2))))

# Compute effect size standard errors
df$se <- sqrt(df$variance_d)

# Define article and study ids
df <- transform(df, article_id = match(authors, unique(authors)))
df$study_id <- paste(df$article_id, df$study_nr, sep = "-")

# Define AI development phases
df <- df %>%
  mutate(ai_development_phase = case_when(
    year >= 2002 & year <= 2009 ~ "before 2010",
    year >= 2010 & year <= 2014 ~ "2010-2014",
    year >= 2015 & year <= 2019 ~ "2015-2019",
    year >= 2020 & year <= 2025 ~ "2020-2025",
    TRUE ~ "< 2002"
  ))

df$ai_development_phase <- factor(df$ai_development_phase,
                          levels = c("before 2010", "2010-2014", "2015-2019", "2020-2025"))

# Specify column types
df$average_age <- as.numeric(as.character(df$average_age))
df$percentage_females <- as.numeric(as.character(df$percentage_females))
df$n1 <- as.numeric(as.character(df$n1))
df$n2 <- as.numeric(as.character(df$n2))
df$n <- as.numeric(as.character(df$n))

# Set plotting theme
theme <- theme(
  text = element_text(family = "Arial", size = 9),
  axis.text.x = element_text(angle = 45, hjust = 1, size = 9),
  axis.text.y = element_text(size = 9),
  axis.title.x = element_text(size = 9, margin = margin(t = 5)),
  axis.title.y = element_text(size = 9, margin = margin(r = 5))
)

# Define plotting colors
blue_light <- "#c6dbef"
blue_medium <- "#6baed6"
blue_dark <- "#08519c"
red_line <- "#a50f15"

# Create function to print model coefficients
create_coefficient_table <- function(model) {
  tidy_model <- broom::tidy(model) %>%
    select(term, estimate, std.error, statistic, p.value) %>%
    mutate(
      estimate = round(estimate, 2),
      std.error = round(std.error, 2),
      statistic = round(statistic, 2),
      p.value = case_when(
        p.value < 0.001 ~ "< .001 ***",
        p.value < 0.01 ~ sprintf("%.3f **", p.value),
        p.value < 0.05 ~ sprintf("%.3f *", p.value),
        TRUE ~ sprintf("%.2f", p.value)
      ),
      ref = str_extract(term, "ref = \".*?\"") %>% 
            str_remove("ref = \"") %>%
            str_remove("\""),
      term = case_when(
        term == "intrcpt" ~ "Intercept",
        TRUE ~ str_extract(term, "[^)]+$") %>%
          str_to_title() %>%
          str_replace("Ai ", "AI ")
      ),
      term = if_else(!is.na(ref) & term != "Intercept",
                    paste0(term, " (vs. ", str_to_title(ref) %>% str_replace("Ai ", "AI "), ")"),
                    term)
    ) %>%
    select(-ref)
  
  kable(tidy_model, 
        col.names = c("Predictor", "Estimate", "SE", "t-value", "p-value"),
        align = c('l', 'r', 'r', 'r', 'r'),
        format = "html") %>%
    kable_styling(
      bootstrap_options = c("striped", "hover"),
      full_width = TRUE,
      position = "left"
    ) %>%
    column_spec(1, bold = TRUE) %>%
    column_spec(2:5, width = "1.2in") %>%
    add_header_above(c(" " = 1, "Parameter Estimates" = 4)) %>%
    footnote(
      general = "Significance levels: *** p < .001; ** p < .01; * p < .05",
      general_title = ""
    )
}
```

# DESCRIPTIVE STATISTICS

## Sample Overview
```{r basic_stats}
k_publications <- length(unique(df$article_id))
k_studies <- length(unique(df$study_id))
k_observations <- length(unique(df$es_id))
n_total <- df %>%
  group_by(study_id) %>%
  slice(1) %>%
  ungroup() %>%
  summarize(n_total = sum(n))

stats_data <- data.frame(
  Metric = c("Number of Articles", "Number of Studies", "Number of Effect Sizes", "Number of Unique Study Subjects"),
  Count = c(k_publications, k_studies, k_observations, n_total$n_total)
)

knitr::kable(stats_data, caption = "Data Scope", align = c('l', 'r')) %>%
  kable_styling(full_width = TRUE) %>%
  column_spec(1, width = "4in") %>%
  column_spec(2, width = "1.5in")

# Disciplines
disc_dist <- df %>%
  count(discipline, .drop = FALSE) %>%
  mutate(percentage = round(n/sum(n) * 100, 1)) %>%
  arrange(desc(percentage)) %>%
  arrange(is.na(discipline), desc(percentage))

names(disc_dist) <- c("Discipline", "Count", "Percentage")
disc_dist$Discipline <- tools::toTitleCase(as.character(disc_dist$Discipline))
disc_dist$Discipline[is.na(disc_dist$Discipline)] <- "NA"

knitr::kable(disc_dist, align = c('l', 'r', 'r')) %>%
  kable_styling(full_width = TRUE) %>%
  column_spec(1, width = "4in") %>%
  column_spec(2, width = "1.5in") %>%
  column_spec(3, width = "1.5in")
```

## AI Label
```{r ai_label}
ai_labels <- df %>%
  count(ai_label, .drop = FALSE) %>%
  mutate(percentage = round(n/sum(n) * 100, 1)) %>%
  arrange(desc(percentage)) %>%
  arrange(is.na(ai_label), desc(percentage))

names(ai_labels) <- c("Label", "Count", "Percentage")
ai_labels$Label <- gsub("Ai", "AI", tools::toTitleCase(as.character(ai_labels$Label)))
ai_labels$Label[is.na(ai_labels$Label)] <- "NA"

knitr::kable(ai_labels, align = c('l', 'r', 'r')) %>%
  kable_styling(full_width = TRUE) %>%
  column_spec(1, width = "4in") %>%
  column_spec(2, width = "1.5in") %>%
  column_spec(3, width = "1.5in")
```

## Consumer Response
```{r consumer_response}
consumer_responses <- df %>%
  count(consumer_response, .drop = FALSE) %>%
  mutate(percentage = round(n/sum(n) * 100, 1)) %>%
  arrange(desc(percentage)) %>%
  arrange(is.na(consumer_response), desc(percentage))

names(consumer_responses) <- c("Response", "Count", "Percentage")
consumer_responses$Response <- tools::toTitleCase(as.character(consumer_responses$Response))
consumer_responses$Response[is.na(consumer_responses$Response)] <- "NA"

knitr::kable(consumer_responses, caption = "Consumer Response", align = c('l', 'r', 'r')) %>%
  kable_styling(full_width = TRUE) %>%
  column_spec(1, width = "4in") %>%
  column_spec(2, width = "1.5in") %>%
  column_spec(3, width = "1.5in")
```

## AI Characteristics
```{r ai_characteristics}
# AI Autonomy
autonomy <- df %>%
  count(autonomy, .drop = FALSE) %>%
  mutate(percentage = round(n/sum(n) * 100, 1)) %>%
  arrange(desc(percentage)) %>%
  arrange(is.na(autonomy), desc(percentage))

names(autonomy) <- c("Level", "Count", "Percentage")
autonomy$Level <- tools::toTitleCase(as.character(autonomy$Level))
autonomy$Level[is.na(autonomy$Level)] <- "NA"

knitr::kable(autonomy, caption = "Autonomy", align = c('l', 'r', 'r')) %>%
  kable_styling(full_width = TRUE) %>%
  column_spec(1, width = "4in") %>%
  column_spec(2, width = "1.5in") %>%
  column_spec(3, width = "1.5in")

# AI Performance
performance <- df %>%
  count(performance, .drop = FALSE) %>%
  mutate(percentage = round(n/sum(n) * 100, 1)) %>%
  arrange(desc(percentage)) %>%
  arrange(is.na(performance), desc(percentage))

names(performance) <- c("Level", "Count", "Percentage")
performance$Level <- tools::toTitleCase(as.character(performance$Level))
performance$Level[is.na(performance$Level)] <- "NA"

knitr::kable(performance, caption = "Performance", align = c('l', 'r', 'r')) %>%
  kable_styling(full_width = TRUE) %>%
  column_spec(1, width = "4in") %>%
  column_spec(2, width = "1.5in") %>%
  column_spec(3, width = "1.5in")

# Anthropomorphism
anthropomorphism <- df %>%
  count(anthropomorphism, .drop = FALSE) %>%
  mutate(percentage = round(n/sum(n) * 100, 1)) %>%
  arrange(desc(percentage)) %>%
  arrange(is.na(anthropomorphism), desc(percentage))

names(anthropomorphism) <- c("Level", "Count", "Percentage")
anthropomorphism$Level <- tools::toTitleCase(as.character(anthropomorphism$Level))
anthropomorphism$Level[is.na(anthropomorphism$Level)] <- "NA"

knitr::kable(anthropomorphism, caption = "Anthropomorphism", align = c('l', 'r', 'r')) %>%
  kable_styling(full_width = TRUE) %>%
  column_spec(1, width = "4in") %>%
  column_spec(2, width = "1.5in") %>%
  column_spec(3, width = "1.5in")

# Uniqueness Recognition
uniqueness <- df %>%
  count(uniqueness_recognition, .drop = FALSE) %>%
  mutate(percentage = round(n/sum(n) * 100, 1)) %>%
  arrange(desc(percentage)) %>%
  arrange(is.na(uniqueness_recognition), desc(percentage))

names(uniqueness) <- c("Level", "Count", "Percentage")
uniqueness$Level <- tools::toTitleCase(as.character(uniqueness$Level))
uniqueness$Level[is.na(uniqueness$Level)] <- "NA"

knitr::kable(uniqueness, caption = "Uniqueness Recognition", align = c('l', 'r', 'r')) %>%
  kable_styling(full_width = TRUE) %>%
  column_spec(1, width = "4in") %>%
  column_spec(2, width = "1.5in") %>%
  column_spec(3, width = "1.5in")
```

## Application Domain
```{r application_domain}
domains <- df %>%
  count(application_domain, .drop = FALSE) %>%
  mutate(percentage = round(n/sum(n) * 100, 1)) %>%
  arrange(desc(percentage)) %>%
  arrange(is.na(application_domain), desc(percentage))

names(domains) <- c("Domain", "Count", "Percentage")
domains$Domain <- tools::toTitleCase(as.character(domains$Domain))
domains$Domain[is.na(domains$Domain)] <- "NA"

knitr::kable(domains, align = c('l', 'r', 'r')) %>%
  kable_styling(full_width = TRUE) %>%
  column_spec(1, width = "4in") %>%
  column_spec(2, width = "1.5in") %>%
  column_spec(3, width = "1.5in")
```

## Temporal AI Advancement
```{r temporal_advancement}
# Publication Years Range
year_range <- data.frame(
  Years = paste(min(df$year, na.rm = TRUE), "—", max(df$year, na.rm = TRUE))
)
knitr::kable(year_range, col.names = NULL, caption = "Range of Publication Years") %>%
  kable_styling(full_width = TRUE) %>%
  column_spec(1, width = "4in")

# AI Development Phases
dev_phases <- df %>%
  count(ai_development_phase, .drop = FALSE) %>%
  mutate(percentage = round(n/sum(n) * 100, 1)) %>%
  arrange(desc(percentage)) %>%
  arrange(is.na(ai_development_phase), desc(percentage))

names(dev_phases) <- c("Development Phase", "Count", "Percentage")
dev_phases$`Development Phase`[is.na(dev_phases$`Development Phase`)] <- "NA"

knitr::kable(dev_phases, caption = "AI Development Phases", align = c('l', 'r', 'r')) %>%
  kable_styling(full_width = TRUE) %>%
  column_spec(1, width = "4in") %>%
  column_spec(2, width = "1.5in") %>%
  column_spec(3, width = "1.5in")
```

## Task Characteristics
```{r task_characteristics, fig.height = 8, fig.width = 8}
# Task Objectivity
task_obj <- df %>%
  count(task_objectivity, .drop = FALSE) %>%
  mutate(percentage = round(n/sum(n) * 100, 1)) %>%
  arrange(desc(percentage)) %>%
  arrange(is.na(task_objectivity), desc(percentage))

names(task_obj) <- c("Level", "Count", "Percentage")
task_obj$Level <- tools::toTitleCase(as.character(task_obj$Level))
task_obj$Level[is.na(task_obj$Level)] <- "NA"

knitr::kable(task_obj, caption = "Task Objectivity", align = c('l', 'r', 'r')) %>%
  kable_styling(full_width = TRUE) %>%
  column_spec(1, width = "4in") %>%
  column_spec(2, width = "1.5in") %>%
  column_spec(3, width = "1.5in")

# Task Consequentiality
task_cons <- df %>%
  count(task_consequentiality, .drop = FALSE) %>%
  mutate(percentage = round(n/sum(n) * 100, 1)) %>%
  arrange(desc(percentage)) %>%
  arrange(is.na(task_consequentiality), desc(percentage))

names(task_cons) <- c("Level", "Count", "Percentage")
task_cons$Level <- tools::toTitleCase(as.character(task_cons$Level))
task_cons$Level[is.na(task_cons$Level)] <- "NA"

knitr::kable(task_cons, caption = "Task Consequentiality", align = c('l', 'r', 'r')) %>%
  kable_styling(full_width = TRUE) %>%
  column_spec(1, width = "4in") %>%
  column_spec(2, width = "1.5in") %>%
  column_spec(3, width = "1.5in")
```

## Study Characteristics
```{r study_characteristics}
# Human Benchmark
benchmark <- df %>%
  count(human_benchmark, .drop = FALSE) %>%
  mutate(percentage = round(n/sum(n) * 100, 1)) %>%
  arrange(desc(percentage)) %>%
  arrange(is.na(human_benchmark), desc(percentage))

names(benchmark) <- c("Level", "Count", "Percentage")
benchmark$Level <- tools::toTitleCase(as.character(benchmark$Level))
benchmark$Level[is.na(benchmark$Level)] <- "NA"

knitr::kable(benchmark, caption = "Human Benchmark", align = c('l', 'r', 'r')) %>%
  kable_styling(full_width = TRUE) %>%
  column_spec(1, width = "4in") %>%
  column_spec(2, width = "1.5in") %>%
  column_spec(3, width = "1.5in")

# AI Explanation
explanation <- df %>%
  count(ai_explanation, .drop = FALSE) %>%
  mutate(percentage = round(n/sum(n) * 100, 1)) %>%
  arrange(desc(percentage)) %>%
  arrange(is.na(ai_explanation), desc(percentage))

names(explanation) <- c("Level", "Count", "Percentage")
explanation$Level <- tools::toTitleCase(as.character(explanation$Level))
explanation$Level[is.na(explanation$Level)] <- "NA"

knitr::kable(explanation, caption = "AI Explanation", align = c('l', 'r', 'r')) %>%
  kable_styling(full_width = TRUE) %>%
  column_spec(1, width = "4in") %>%
  column_spec(2, width = "1.5in") %>%
  column_spec(3, width = "1.5in")

# Field Setting
field <- df %>%
  count(field_setting, .drop = FALSE) %>%
  mutate(percentage = round(n/sum(n) * 100, 1)) %>%
  arrange(desc(percentage)) %>%
  arrange(is.na(field_setting), desc(percentage))

names(field) <- c("Level", "Count", "Percentage")
field$Level <- tools::toTitleCase(as.character(field$Level))
field$Level[is.na(field$Level)] <- "NA"

knitr::kable(field, caption = "Field Setting", align = c('l', 'r', 'r')) %>%
  kable_styling(full_width = TRUE) %>%
  column_spec(1, width = "4in") %>%
  column_spec(2, width = "1.5in") %>%
  column_spec(3, width = "1.5in")

# Stimulus Presentation
stimulus <- df %>%
  count(stimulus_presentation, .drop = FALSE) %>%
  mutate(percentage = round(n/sum(n) * 100, 1)) %>%
  arrange(desc(percentage)) %>%
  arrange(is.na(stimulus_presentation), desc(percentage))

names(stimulus) <- c("Level", "Count", "Percentage")
stimulus$Level <- tools::toTitleCase(as.character(stimulus$Level))
stimulus$Level[is.na(stimulus$Level)] <- "NA"

knitr::kable(stimulus, caption = "Stimulus Presentation") %>%
  kable_styling(full_width = TRUE) %>%
  column_spec(1, width = "4in") %>%
  column_spec(2, width = "1.5in") %>%
  column_spec(3, width = "1.5in")

# Incentive Compatibility
incentive <- df %>%
  count(incentive_compatibility, .drop = FALSE) %>%
  mutate(percentage = round(n/sum(n) * 100, 1)) %>%
  arrange(desc(percentage)) %>%
  arrange(is.na(incentive_compatibility), desc(percentage))

names(incentive) <- c("Level", "Count", "Percentage")
incentive$Level <- tools::toTitleCase(as.character(incentive$Level))
incentive$Level[is.na(incentive$Level)] <- "NA"

knitr::kable(incentive, caption = "Incentive Compatibility") %>%
  kable_styling(full_width = TRUE) %>%
  column_spec(1, width = "4in") %>%
  column_spec(2, width = "1.5in") %>%
  column_spec(3, width = "1.5in")

# Behavioral DV
behavioral <- df %>%
  count(behavioral_dv, .drop = FALSE) %>%
  mutate(percentage = round(n/sum(n) * 100, 1)) %>%
  arrange(desc(percentage)) %>%
  arrange(is.na(behavioral_dv), desc(percentage))

names(behavioral) <- c("Level", "Count", "Percentage")
behavioral$Level <- tools::toTitleCase(as.character(behavioral$Level))
behavioral$Level[is.na(behavioral$Level)] <- "NA"

knitr::kable(behavioral, caption = "Behavioral DV") %>%
  kable_styling(full_width = TRUE) %>%
  column_spec(1, width = "4in") %>%
  column_spec(2, width = "1.5in") %>%
  column_spec(3, width = "1.5in")

# Within-Subject Design
within_subj <- df %>%
  count(within_subject_design, .drop = FALSE) %>%
  mutate(percentage = round(n/sum(n) * 100, 1)) %>%
  arrange(desc(percentage)) %>%
  arrange(is.na(within_subject_design), desc(percentage))

names(within_subj) <- c("Level", "Count", "Percentage")
within_subj$Level <- tools::toTitleCase(as.character(within_subj$Level))
within_subj$Level[is.na(within_subj$Level)] <- "NA"

knitr::kable(within_subj, caption = "Within-Subject Design") %>%
  kable_styling(full_width = TRUE) %>%
  column_spec(1, width = "4in") %>%
  column_spec(2, width = "1.5in") %>%
  column_spec(3, width = "1.5in")

# Online Sample
online <- df %>%
  count(online_sample, .drop = FALSE) %>%
  mutate(percentage = round(n/sum(n) * 100, 1)) %>%
  arrange(desc(percentage)) %>%
  arrange(is.na(online_sample), desc(percentage))

names(online) <- c("Level", "Count", "Percentage")
online$Level <- tools::toTitleCase(as.character(online$Level))
online$Level[is.na(online$Level)] <- "NA"

knitr::kable(online, caption = "Online Sample") %>%
  kable_styling(full_width = TRUE) %>%
  column_spec(1, width = "4in") %>%
  column_spec(2, width = "1.5in") %>%
  column_spec(3, width = "1.5in")
```

## Sample Characteristics
```{r sample_characteristics}
# Age and Gender
age_gender_avail <- data.frame(
  Characteristic = c(
    "Age Data Available",
    "Age Data Missing",
    "Gender Data Available",
    "Gender Data Missing"
  ),
  Count = c(
    length(na.omit(df$average_age)),
    sum(is.na(df$average_age)),
    length(na.omit(df$percentage_females)),
    sum(is.na(df$percentage_females))
  )
) %>%
  mutate(Percentage = round(Count/sum(Count[1:2]) * 100, 1))

names(age_gender_avail) <- c("Availability", "Count", "Percentage")

knitr::kable(age_gender_avail, caption = "Age and Gender") %>%
  kable_styling(full_width = TRUE) %>%
  column_spec(1, width = "4in") %>%
  column_spec(2, width = "1.5in") %>%
  column_spec(3, width = "1.5in")

demo_means <- data.frame(
  Metric = c("Mean Age", "Mean % Female"),
  Value = c(
    round(mean(df$average_age, na.rm = TRUE), 1),
    round(mean(df$percentage_females, na.rm = TRUE), 1)
  )
)
names(demo_means) <- c("Variable", "Percentage")
knitr::kable(demo_means, caption = "") %>% 
  kable_styling(full_width = TRUE) %>% 
  column_spec(1, width = "4in") %>% 
  column_spec(2, width = "1.5in")

# Geographical Background
geo_dist <- df %>%
  count(geographical_background, .drop = FALSE) %>%
  mutate(percentage = round(n/sum(n) * 100, 1)) %>%
  arrange(desc(percentage)) %>%
  arrange(is.na(geographical_background), desc(percentage))

names(geo_dist) <- c("Region", "Count", "Percentage")
geo_dist$Region <- tools::toTitleCase(as.character(geo_dist$Region))
geo_dist$Region[is.na(geo_dist$Region)] <- "NA"

knitr::kable(geo_dist, caption = "Geographic Distribution of Samples") %>%
  kable_styling(full_width = TRUE) %>%
  column_spec(1, width = "4in") %>%
  column_spec(2, width = "1.5in") %>%
  column_spec(3, width = "1.5in")
```

# PRIMARY META-ANALYSIS

## Overall Effect Size Estimate
```{r meta_analysis}
# Overall effect size.
overall <- metafor::rma.mv(cohens_d, variance_d, data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, slab = authors)

# Robust test and confidence intervals accounting for dependence in sampling errors (i.e., clusters within studies due to overlap of participants)
overall_robust <- robust(overall, cluster = study_id)
summary(overall_robust)
```

## Prediction Interval
```{r prediction_interval}
predict_overall <- predict(overall_robust)
predict_overall
```

## Forest Plot
```{r forest_plot, fig.height = 8, fig.width = 9}
forest(overall_robust,
       addfit = FALSE,
       xlim=c(-3,3),
       order="obs",
       slab=NA, annotate=FALSE,
       efac=0,
       pch=19,
       col="black",
       psize=2,
       cex.lab=1, cex.axis=1,
       xlab = "",
       lwd = 0.5)
addpoly.rma(overall_robust, addpred=TRUE, mlab="", efac = 4, lty = "solid")
text(-3, -7, expression("Model estimate with prediction\ninterval: d = -0.21*** [-1.20, 0.78]"), pos=4, offset=0, cex=1)
mtext("Effect size (Cohen's d) with 95% CI", side=1, line=3, cex=1)
```

## Sensitivity Analysis
We set `eval = FALSE` for this part of the analyses to keep the document concise. The full analyses can be found in the source script.
```{r sensitivity_analysis, eval = FALSE}
# Standardized residuals.
resid <- rstandard(overall, progbar = TRUE, reestimate = TRUE)
outliers_resid <- data.frame(
  es_id = overall$data$es_id[resid$z > 3 | resid$z < -3],
  author = overall$data$author[resid$z > 3 | resid$z < -3],
  z_score = resid$z[resid$z > 3 | resid$z < -3]
) %>%
  arrange(desc(abs(z_score)))
print(outliers_resid)

# Cook's distance.
cooks_distance <- cooks.distance.rma.mv(overall, progbar = TRUE, reestimate = TRUE)
cutoff <- 4 / (440 - 1 - 1)
outliers_cooks <- data.frame(
  es_id = overall$data$es_id[cooks_distance > cutoff],
  author = overall$data$author[cooks_distance > cutoff],
  cooks_d = cooks_distance[cooks_distance > cutoff]
) %>%
  arrange(desc(cooks_d))
print(outliers_cooks)

unique_outliers <- unique(c(outliers_resid$es_id, outliers_cooks$es_id))
unique_outliers <- sort(unique_outliers)
unique_outliers
outliers_cohens_d <- df$cohens_d[df$es_id %in% unique_outliers]
range(outliers_cohens_d)

# overall meta-analysis w/out influential outliers.
overall_sens <- metafor::rma.mv(cohens_d, variance_d, data = df, subset = c(-67, -166, -170, -175, -176, -189, -190, -209, -251, -253, -288, -341, -342), random = ~ 1 | article_id/es_id, tdist = TRUE, slab = authors)
summary(overall_sens)

# Robust test and confidence intervals accounting for dependence in sampling errors.
overall_sens_robust <- robust(overall_sens, cluster = study_id)
summary(overall_sens_robust)

# Leave-one-out analysis at effect size level.
df_loo <- matrix(0, ncol = 12, nrow = nrow(df))
df_loo <- data.frame(df_loo)
colnames(df_loo) <- c("estimate", "se", "tval", "pval", "ci.lb", "ci.ub", "sigma.total", "sigma2.1",
                      "sigma2.2", "i2total", "i2between", "i2within")
for(i in 1:nrow(df))
{
  df1 <- df[c(-i),]
  loo <- metafor::rma.mv(cohens_d, variance_d, data = df1, random = ~ 1 | article_id/es_id, tdist = TRUE)
  loo_robust <- robust(loo, cluster = df1$study_id, adjust = TRUE)
  df_loo[i, 1] <- loo_robust$b
  df_loo[i, 2] <- loo_robust$se
  df_loo[i, 4] <- loo_robust$pval
  df_loo[i, 5] <- loo_robust$ci.lb
  df_loo[i, 6] <- loo_robust$ci.ub
  df_loo[i, 7] <- loo_robust$sigma2[1] + loo_robust$sigma2[2]
  df_loo[i, 8] <- loo_robust$sigma2[1]
  df_loo[i, 9] <- loo_robust$sigma2[2]
  w <- diag(1/df1$variance_d)
  x <- model.matrix(loo_robust)
  p <- w - w %*% x %*% solve(t(x) %*% w %*% x) %*% t(x) %*% w
  i2.ml <- 100 * loo_robust$sigma2 / (sum(loo_robust$sigma2) + (loo_robust$k - loo_robust$p)/sum(diag(p)))
  df_loo[i, 10] <- i2.ml[1] + i2.ml[2]
  df_loo[i, 11] <- i2.ml[1]
  df_loo[i, 12] <- i2.ml[2]
}
range(df_loo$estimate)
range(df_loo$pval)

# Leave-one-out analysis at publication size level.
df_loopub <- matrix(0, ncol = 13, nrow = 72)
df_loopub <- data.frame(df_loopub)
colnames(df_loopub) <- c("numES", "estimate", "se", "tval", "pval", "ci.lb", "ci.ub", "sigma.total",
                         "sigma2.1", "sigma2.2", "i2total", "i2between", "i2within")

for(i in 1:72)
{
  df2 <- subset(df, article_id != i)
  loopub <- metafor::rma.mv(cohens_d, variance_d, data = df2, random = ~ 1 | article_id/es_id, tdist = TRUE)
  loopub_robust <- robust(loopub, cluster = df2$study_id, adjust = TRUE)
  df_loopub[i, 1] <- nrow(df) - nrow(df2)
  df_loopub[i, 2] <- loopub_robust$b
  df_loopub[i, 3] <- loopub_robust$se
  df_loopub[i, 5] <- loopub_robust$pval
  df_loopub[i, 6] <- loopub_robust$ci.lb
  df_loopub[i, 7] <- loopub_robust$ci.ub
  df_loopub[i, 8] <- loopub_robust$sigma2[1] + loopub_robust$sigma2[2]
  df_loopub[i, 9] <- loopub_robust$sigma2[1]
  df_loopub[i, 10] <- loopub_robust$sigma2[2]
  w <- diag(1/df2$variance_d)
  x <- model.matrix(loopub_robust)
  p <- w - w %*% x %*% solve(t(x) %*% w %*% x) %*% t(x) %*% w
  i2.ml <- 100 * loopub_robust$sigma2 / (sum(loopub_robust$sigma2) + (loopub_robust$k - loopub_robust$p)/sum(diag(p)))
  df_loopub[i,11] <- i2.ml[1]+i2.ml[2]
  df_loopub[i,12] <- i2.ml[1]
  df_loopub[i,13] <- i2.ml[2]
}
range(df_loopub$estimate)
range(df_loopub$pval)
```

## Publication Bias
We set `eval = FALSE` for this part of the analyses to keep the document concise. The full analyses can be found in the source script.
```{r publication_bias, eval = FALSE}
# Funnel plot.
funnel(overall,
       yaxis = "sei",
       main = "",
       xlab = "Effect Size (Cohen's d)",
       ylab = "Standard Error",
       cex = 0.5,
       cex.axis = 0.6,
       cex.lab = 0.6,
       back = "white")

# Rank correlation test.
ranktest(overall)

# Egger's test of asymmetry.
overall_eggers <- metafor::rma.mv(cohens_d, variance_d, mod = ~ sqrt(variance_d), random = ~ 1 | article_id/es_id, data = df)
overall_eggers_robust <- robust(overall_eggers, cluster = study_id, adjust = TRUE)
summary(overall_eggers_robust)

# Trim and fill analysis.
overall_rma <- rma.uni(cohens_d, variance_d, data = df)
summary(overall_rma)
trimfill(overall_rma, "left")
trimfill(overall_rma, "right")
taf <- trimfill(overall_rma, "left")
funnel(taf,
       yaxis = "sei",
       main = "",
       xlab = "Effect Size (Cohen's d)",
       ylab = "Standard Error",
       cex = 0.5,
       cex.axis = 0.6,
       cex.lab = 0.6,
       back = "white")

# PET (Precision Effect Test).
pet_model <- rma.mv(cohens_d, variance_d,
                    mods = ~ sqrt(variance_d),
                    random = ~ 1 | article_id/es_id,
                    data = df)
pet_robust <- robust(pet_model, cluster = study_id)
pet_robust

# PEESE (Precision Effect Estimate with Standard Error)
peese_model <- rma.mv(cohens_d, variance_d,
                      mods = ~ variance_d,
                      random = ~ 1 | article_id/es_id,
                      data = df)
peese_robust <- robust(peese_model, cluster = study_id)
peese_robust

# Selection Models.
tab <- data.frame(
steps = c(0.005, 0.01, 0.05, 0.10, 0.25, 0.35, 0.50, 0.65, 0.75, 0.90, 0.95, 0.99, 0.995, 1),
delta.mod.1 = c(1, 0.99, 0.95, 0.80, 0.75, 0.65, 0.60, 0.55, 0.50, 0.50, 0.50, 0.50, 0.50, 0.50),
delta.sev.1 = c(1, 0.99, 0.90, 0.75, 0.60, 0.50, 0.40, 0.35, 0.30, 0.25, 0.10, 0.10, 0.10, 0.10))

overall_rma <- rma.uni(cohens_d, variance_d, data = df)
summary(overall_rma)

# Moderate one-tailed selection bias scenario.
overall_mod1 <- metafor::selmodel(overall_rma, type = 'stepfun', alternative = 'greater', delta = tab$delta.mod.1, steps = tab$steps)
summary(overall_mod1)

# Severe one-tailed selection bias scenario.
overall_sev1 <- metafor::selmodel(overall_rma, type = 'stepfun', alternative = 'greater', delta = tab$delta.sev.1, steps = tab$steps)
summary(overall_sev1)

# Fail-safe N.
(5*440)+10 # Threshold
fsn(cohens_d, vi=variance_d, data = df, type = "Rosenthal")

# Ruecker's Limit.
rucker_limit <- function(rma_model, side = "left", maxiter = 100) {
    prev_est <- rma_model$b[1]
    current_model <- rma_model
    cat("\nStarting Ruckers Limit analysis...\n")
    for(i in 1:maxiter) {
        progress <- round(i/maxiter * 100, 1)
        cat(sprintf("\rIteration %d/%d (%g%%)", i, maxiter, progress))
        utils::flush.console()  # Ensure output is displayed immediately
        tf <- trimfill(current_model, side = side)
        current_est <- tf$b[1]
        if(abs(current_est - prev_est) < 0.0001) {
            cat("\nConverged successfully at iteration", i, "!\n")
            return(tf)
        }
        prev_est <- current_est
        current_model <- tf
    }
    warning("Maximum iterations reached without convergence")
    return(current_model)
}
rucker_result <- rucker_limit(overall_rma, side = "left")

P-Curve analysis.

# Step 1: Clean the p_value column
df_coded <- df %>%
  mutate(
    # Remove any text in parentheses
    p_value_cleaned = gsub("\\s*\\(.*\\)", "", p_value),

    # Convert "<" and ">" indicators to approximate values
    p_value_cleaned = case_when(
      grepl("< 0.001", p_value_cleaned) ~ "0.001",
      grepl("< 0.01", p_value_cleaned) ~ "0.01",
      grepl("< 0.05", p_value_cleaned) ~ "0.05",
      grepl("< 0.1", p_value_cleaned) ~ "0.1",
      grepl("> 0.1", p_value_cleaned) ~ "0.1",
      grepl("n\\.s\\.", p_value_cleaned, ignore.case = TRUE) ~ NA_character_,
      TRUE ~ p_value_cleaned
    ),

    # Convert cleaned p-values to numeric
    p_value_numeric = as.numeric(p_value_cleaned)
  ) %>%
  filter(!is.na(p_value_numeric) & p_value_numeric < 0.05) %>%
  mutate(p_value_capped = ifelse(p_value_numeric < 1e-8, 1e-8, p_value_numeric))

# Step 2: Create bins for p-values and summarize proportions
p_curve_data_coded <- df_coded %>%
  mutate(p_bin = cut(p_value_capped,
                     breaks = c(0, 0.01, 0.02, 0.03, 0.04, 0.05),
                     labels = c(".01", ".02", ".03", ".04", ".05"),
                     include.lowest = TRUE)) %>%
  group_by(p_bin) %>%
  summarise(n = n()) %>%
  mutate(prop = n / sum(n), method = "Coded p-values")

p_curve <- ggplot(p_curve_data_coded,
                  aes(x = p_bin, y = prop)) +
  geom_bar(stat = "identity", width = 0.7) +
  labs(x = expression(italic("p")*"-value"),
       y = "Percentage of Test Results") +
  theme_minimal() +
  theme(
  text = element_text(family = "Arial", size = 9),
  axis.text.x = element_text(hjust = 1, size = 9),
  axis.text.y = element_text(size = 9),
  axis.title.x = element_text(size = 9, margin = margin(t = 5)),
  axis.title.y = element_text(size = 9, margin = margin(r = 5))
  ) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1))
p_curve
```

## Heterogeneity
```{r heterogeneity}
# Range of cohen's d values.
cohens_d_range <- range(df$cohens_d)

# Total heterogeneity across studies (Tau squared).
tau_squared <- sum(overall_robust$sigma2)

# Ratio of true to total variance across studies (I squared).
w <- diag(1 / df$variance_d)
x <- model.matrix(overall_robust)
p <- w - w %*% x %*% solve(t(x) %*% w %*% x) %*% t(x) %*% w
i_square <- 100 * sum(overall_robust$sigma2) / (sum(overall_robust$sigma2) + (overall_robust$k-overall_robust$p)/sum(diag(p)))

# Display results with clear formatting
cat("Range of Effect Sizes (Cohen's d):", 
    sprintf("[%.2f, %.2f]", cohens_d_range[1], cohens_d_range[2]), "\n\n")
cat("Total Heterogeneity (τ²):", 
    sprintf("%.2f", tau_squared), "\n\n")
cat("Ratio of True to Total Variance (I²):", 
    sprintf("%.2f%%", i_square), "\n")
```

# META-REGRESSIONS

## AI Labels & Consumer Responses
```{r}
############## MODERATOR ANALYSIS AI LABEL ##############
# Effect size estimates for each ai label.
ai_label <- metafor::rma.mv(cohens_d, variance_d, mods = ~ factor(ai_label) - 1, data = df, random = ~ 1 | article_id/es_id, tdist = TRUE)
ai_label_robust <- robust(ai_label, cluster = study_id, adjust=TRUE)

# Comparison of effect sizes across ai label categories (reference category: ai algorithms). 
ai_algorithms <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(ai_label), ref = 'ai algorithms'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:4)
ai_algorithms_robust <- robust(ai_algorithms, cluster = study_id)

# Comparison of effect sizes across ai label categories (reference category: ai systems). 
ai_systems <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(ai_label), ref = 'ai systems'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:4)
ai_systems_robust <- robust(ai_systems, cluster = study_id)

# Comparison of effect sizes across ai label categories (reference category: ai robots). 
ai_robots <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(ai_label), ref = 'ai robots'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:4)
ai_robots_robust <- robust(ai_robots, cluster = study_id)

# AI label plot.
coef_order <- coef(ai_label_robust)
ai_label_plot <- df %>%
  filter(!is.na(ai_label)) %>%
  mutate(ai_label = fct_recode(fct_reorder(ai_label, coef_order[paste0("factor(ai_label)", ai_label)], .desc = TRUE),
                               "AI Algorithms" = "ai algorithms",
                               "AI Systems" = "ai systems",
                               "AI Assistants" = "ai assistants",
                               "AI Robots" = "ai robots")) %>%
  ggplot(aes(x = ai_label, y = cohens_d)) +
  geom_violin(trim = FALSE, fill = blue_light, color = blue_dark) +
  geom_boxplot(width = 0.1, fill = blue_medium, color = blue_dark, 
               outlier.shape = NA) +
  stat_summary(fun = median, geom = "point", shape = 21, size = 2,
               fill = "white", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", size = 1, 
             color = red_line) +
  scale_y_continuous(limits = c(-3, 2), breaks = seq(-3, 2, 1)) +
  labs(x = "AI Label", y = "Cohen's d") +
  theme_classic() +
  theme

############## MODERATOR ANALYSIS CONSUMER RESPONSE ##############
# Effect size estimates for each consumer response type.
consumer_response <- metafor::rma.mv(cohens_d, variance_d, mods = ~ factor(consumer_response) - 1, data = df, random = ~ 1 | article_id/es_id, tdist = TRUE)
consumer_response_robust <- robust(consumer_response, cluster = study_id, adjust=TRUE)

# Comparison of effect sizes across consumer responses (reference category: affective). 
affective <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(consumer_response), ref = 'affective'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:3)
affective_robust <- robust(affective, cluster = study_id)

# Comparison of effect sizes across consumer responses (reference category: behavioral). 
behavioral <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(consumer_response), ref = 'behavioral'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:3)
behavioral_robust <- robust(behavioral, cluster = study_id)

# Consumer response plot.
coef_order <- coef(consumer_response_robust)
consumer_response_plot <- df %>%
  filter(!is.na(consumer_response)) %>%
  mutate(consumer_response = fct_recode(fct_reorder(consumer_response, 
                                                    coef_order[paste0("factor(consumer_response)", consumer_response)], 
                                                    .desc = TRUE),
                                        "Affective" = "affective",
                                        "Behavioral" = "behavioral",
                                        "Cognitive" = "cognitive")) %>%
  ggplot(aes(x = consumer_response, y = cohens_d)) +
  geom_violin(trim = FALSE, fill = blue_light, color = blue_dark) +
  geom_boxplot(width = 0.1, fill = blue_medium, color = blue_dark, 
               outlier.shape = NA) +
  stat_summary(fun = median, geom = "point", shape = 21, size = 2,
               fill = "white", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", size = 1, 
             color = red_line) +
  scale_y_continuous(limits = c(-3, 2), breaks = seq(-3, 2, 1)) +
  labs(x = "Consumer Response", y = "Cohen's d") +
  theme_classic() +
  theme

# Combine ai label and consumer response plots.
label_response_plot <- plot_grid(
  ai_label_plot, consumer_response_plot,
  labels = c("A", "B"),
  label_size = 10,
  ncol = 2
)
```

### Variation in Effect Sizes across AI Labels & Consumer Responses
```{r}
label_response_plot
```

### Effect Size Estimates for Each AI Label
```{r}
# Coefficient table
create_coefficient_table(ai_label_robust)
```

### Comparison of Effect Sizes Across AI Label Categories (Reference Category: AI Robots)
```{r}
# Coefficient table
create_coefficient_table(ai_robots)
```

### Effect Size Estimates for Each Consumer Response
```{r}
# Coefficient table
create_coefficient_table(consumer_response_robust)
```

### Comparison of Effect Sizes Across Consumer Responses (Reference Category: Affective)
```{r}
# Coefficient table
create_coefficient_table(affective_robust)
```

## AI Characteristics
```{r}
############## MODERATOR ANALYSIS AUTONOMY ##############
# Effect size estimates for each level of autonomy.
autonomy <- metafor::rma.mv(cohens_d, variance_d, mods = ~ factor(autonomy) - 1, data = df, random = ~ 1 | article_id/es_id, tdist = TRUE)
autonomy_robust <- robust(autonomy, cluster = study_id, adjust=TRUE)

# Comparison of effect sizes across levels of autonomy (reference category: extensive). 
extensive <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(autonomy), ref = 'extensive'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:2)
extensive_robust <- robust(extensive, cluster = study_id)

# Autonomy plot.
coef_order <- coef(autonomy_robust)
autonomy_plot <- df %>%
  filter(!is.na(autonomy)) %>%
  mutate(autonomy = fct_recode(fct_reorder(autonomy, 
                                           coef_order[paste0("factor(autonomy)", autonomy)], 
                                           .desc = TRUE),
                               "Limited" = "limited",
                               "Extensive" = "extensive")) %>%
  ggplot(aes(x = autonomy, y = cohens_d)) +
  geom_violin(trim = FALSE, fill = blue_light, color = blue_dark) +
  geom_boxplot(width = 0.1, fill = blue_medium, color = blue_dark, 
               outlier.shape = NA) +
  stat_summary(fun = median, geom = "point", shape = 21, size = 2,
               fill = "white", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", size = 1, 
             color = red_line) +
  scale_y_continuous(limits = c(-3, 2), breaks = seq(-3, 2, 1)) +
  labs(x = "Autonomy", y = "Cohen's d") +
  theme_classic() +
  theme

############## MODERATOR ANALYSIS PERFORMANCE ##############
# Effect size estimates for each level of performance.
performance <- metafor::rma.mv(cohens_d, variance_d, mods = ~ factor(performance) - 1, data = df, random = ~ 1 | article_id/es_id, tdist = TRUE)
performance_robust <- robust(performance, cluster = study_id, adjust=TRUE)

# Comparison of effect sizes across level of performance (reference category: inferior). 
inferior <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(performance), ref = 'inferior'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:3)
inferior_robust <- robust(inferior, cluster = study_id)

# Comparison of effect sizes across level of performance (reference category: equivalent). 
equivalent <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(performance), ref = 'equivalent'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:3)
equivalent_robust <- robust(equivalent, cluster = study_id)

# Performance plot.
coef_order <- coef(performance_robust)
performance_plot <- df %>%
  filter(!is.na(performance)) %>%
  mutate(performance = fct_recode(fct_reorder(performance, 
                                              coef_order[paste0("factor(performance)", performance)], 
                                              .desc = TRUE),
                                  "Superior" = "superior",
                                  "Equivalent" = "equivalent",
                                  "Inferior" = "inferior")) %>%
  ggplot(aes(x = performance, y = cohens_d)) +
  geom_violin(trim = FALSE, fill = blue_light, color = blue_dark) +
  geom_boxplot(width = 0.1, fill = blue_medium, color = blue_dark, 
               outlier.shape = NA) +
  stat_summary(fun = median, geom = "point", shape = 21, size = 2,
               fill = "white", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", size = 1, 
             color = red_line) +
  scale_y_continuous(limits = c(-3, 2), breaks = seq(-3, 2, 1)) +
  labs(x = "Performance", y = "Cohen's d") +
  theme_classic() +
  theme

############## MODERATOR ANALYSIS ANTHROPOMORPHISM ##############
# Effect size estimates for each level of anthropomorphism.
anthropomorphism <- metafor::rma.mv(cohens_d, variance_d, mods = ~ factor(anthropomorphism) - 1, data = df, random = ~ 1 | article_id/es_id, tdist = TRUE)
anthropomorphism_robust <- robust(anthropomorphism, cluster = study_id, adjust=TRUE)

# Comparison of effect sizes across levels of anthropomorphism (reference category: present). 
anthropomorphism_present <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(anthropomorphism), ref = 'present'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:2)
anthropomorphism_present_robust <- robust(anthropomorphism_present, cluster = study_id)

# Anthropomorphism plot.
coef_order <- coef(anthropomorphism_robust)
anthropomorphism_plot <- df %>%
  filter(!is.na(anthropomorphism)) %>%
  mutate(anthropomorphism = fct_recode(fct_reorder(anthropomorphism, 
                                                   coef_order[paste0("factor(anthropomorphism)", anthropomorphism)], 
                                                   .desc = TRUE),
                                       "Present" = "present",
                                       "Absent" = "absent")) %>%
  ggplot(aes(x = anthropomorphism, y = cohens_d)) +
  geom_violin(trim = FALSE, fill = blue_light, color = blue_dark) +
  geom_boxplot(width = 0.1, fill = blue_medium, color = blue_dark, 
               outlier.shape = NA) +
  stat_summary(fun = median, geom = "point", shape = 21, size = 2,
               fill = "white", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", size = 1, 
             color = red_line) +
  scale_y_continuous(limits = c(-3, 2), breaks = seq(-3, 2, 1)) +
  labs(x = "Anthropomorphism", y = "Cohen's d") +
  theme_classic() +
  theme

############## MODERATOR ANALYSIS UNIQUENESS RECOGNITION ##############
# Effect size estimates for each level of uniqueness recognition
uniqueness_recognition <- metafor::rma.mv(cohens_d, variance_d, mods = ~ factor(uniqueness_recognition) - 1, data = df, random = ~ 1 | article_id/es_id, tdist = TRUE)
uniqueness_recognition_robust <- robust(uniqueness_recognition, cluster = study_id, adjust=TRUE)

# Comparison of effect sizes across levels of uniqueness recognition (reference category: present). 
uniqueness_recognition_present <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(uniqueness_recognition), ref = 'present'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:2)
uniqueness_recognition_present_robust <- robust(uniqueness_recognition_present, cluster = study_id)

# Uniqueness recognition plot.
coef_order <- coef(uniqueness_recognition_robust)
uniqueness_recognition_plot <- df %>%
  filter(!is.na(uniqueness_recognition)) %>%
  mutate(uniqueness_recognition = fct_recode(fct_reorder(uniqueness_recognition, 
                                                         coef_order[paste0("factor(uniqueness_recognition)", uniqueness_recognition)], 
                                                         .desc = TRUE),
                                             "Present" = "present",
                                             "Absent" = "absent")) %>%
  ggplot(aes(x = uniqueness_recognition, y = cohens_d)) +
  geom_violin(trim = FALSE, fill = blue_light, color = blue_dark) +
  geom_boxplot(width = 0.1, fill = blue_medium, color = blue_dark, 
               outlier.shape = NA) +
  stat_summary(fun = median, geom = "point", shape = 21, size = 2,
               fill = "white", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", size = 1, 
             color = red_line) +
  scale_y_continuous(limits = c(-3, 2), breaks = seq(-3, 2, 1)) +
  labs(x = "Uniqueness\nRecognition", y = "Cohen's d") +
  theme_classic() +
  theme

# Combine AI characteristics plots.
ai_characteristics_plot <- plot_grid(
  autonomy_plot, performance_plot, 
  anthropomorphism_plot, uniqueness_recognition_plot,
  labels = c("A", "B", "C", "D"),
  label_size = 10,
  ncol = 4
)
```

### Variation in Effect Sizes across AI Characteristics
```{r}
ai_characteristics_plot
```

### Effect Size Estimates for Each Level of Autonomy
```{r}
# Coefficient table
create_coefficient_table(autonomy_robust)
```

### Comparison of Effect Sizes Across Levels of Autonomy (Reference Category: Extensive)
```{r}
# Coefficient table
create_coefficient_table(extensive_robust)
```

### Effect Size Estimates for Each Level of Performance
```{r}
# Coefficient table
create_coefficient_table(performance_robust)
```

### Comparison of Effect Sizes Across Levels of Performance (Reference Category: Inferior)
```{r}
# Coefficient table
create_coefficient_table(inferior_robust)
```

### Effect Size Estimates for Each Level of Anthropomorphism
```{r}
# Coefficient table
create_coefficient_table(anthropomorphism_robust)
```

### Comparison of Effect Sizes Across Levels of Anthropomorphism (Reference Category: Present)
```{r}
# Coefficient table
create_coefficient_table(anthropomorphism_present_robust)
```

### Effect Size Estimates for Each Level of Uniqueness Recognition
```{r}
# Coefficient table
create_coefficient_table(uniqueness_recognition_robust)
```

### Comparison of Effect Sizes Across Levels of Uniqueness Recognition (Reference Category: Present)
```{r}
# Coefficient table
create_coefficient_table(uniqueness_recognition_present_robust)
```

## Application Domain
```{r}
############## MODERATOR ANALYSIS APPLICATION DOMAIN ##############
# Effect size estimates for each application domain.
application_domain <- metafor::rma.mv(cohens_d, variance_d, mods = ~ factor(application_domain) - 1, data = df, random = ~ 1 | article_id/es_id, tdist = TRUE)
application_domain_robust <- robust(application_domain, cluster = study_id, adjust=TRUE)

# Comparison of effect sizes across application domains (reference category: entertainment & lifestyle). 
entertainment_lifestyle <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(application_domain), ref = 'entertainment & lifestyle'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:8)
entertainment_lifestyle_robust <- robust(entertainment_lifestyle, cluster = study_id)

# Comparison of effect sizes across application domains (reference category: healthcare). 
healthcare <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(application_domain), ref = 'healthcare'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:8)
healthcare_robust <- robust(healthcare, cluster = study_id)

# Comparison of effect sizes across application domains (reference category: investing & finance). 
investing_finance <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(application_domain), ref = 'investing & finance'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:8)
investing_finance_robust <- robust(investing_finance, cluster = study_id)

# Comparison of effect sizes across application domains (reference category: legal & public safety). 
legal_public_safety <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(application_domain), ref = 'legal & public safety'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:8)
legal_public_safety_robust <- robust(legal_public_safety, cluster = study_id)

# Comparison of effect sizes across application domains (reference category: operations & management). 
operations_management <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(application_domain), ref = 'operations & management'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:8)
operations_management_robust <- robust(operations_management, cluster = study_id)

# Comparison of effect sizes across application domains (reference category: other). 
other <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(application_domain), ref = 'other'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:8)
other_robust <- robust(other, cluster = study_id)

# Comparison of effect sizes across application domains (reference category: social welfare). 
social_welfare <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(application_domain), ref = 'social welfare'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:8)
social_welfare_robust <- robust(social_welfare, cluster = study_id)

# Application domain plot.
coef_order <- coef(application_domain_robust)
application_domain_plot <- df %>%
  filter(!is.na(application_domain)) %>%
  mutate(application_domain = fct_recode(fct_reorder(application_domain, 
                                                     coef_order[paste0("factor(application_domain)", application_domain)], 
                                                     .desc = TRUE),
                                         "Entertainment &\nLifestyle" = "entertainment & lifestyle",
                                         "Healthcare" = "healthcare",
                                         "Investing &\nFinance" = "investing & finance",
                                         "Legal &\nPublic Safety" = "legal & public safety",
                                         "Operations &\nManagement" = "operations & management",
                                         "Other" = "other",
                                         "Social Welfare" = "social welfare",
                                         "Transportation" = "transportation")) %>%
  ggplot(aes(x = application_domain, y = cohens_d)) +
  geom_violin(trim = FALSE, fill = blue_light, color = blue_dark) +
  geom_boxplot(width = 0.1, fill = blue_medium, color = blue_dark, 
               outlier.shape = NA) +
  stat_summary(fun = median, geom = "point", shape = 21, size = 2,
               fill = "white", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", size = 1, 
             color = red_line) +
  scale_y_continuous(limits = c(-3, 2), breaks = seq(-3, 2, 1)) +
  labs(x = "Application Domain", y = "Cohen's d") +
  theme_classic() +
  theme +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
```

### Variation in Effect Sizes across AI Application Domains
```{r}
application_domain_plot
```

### Effect Size Estimates for Each AI Application Domain
```{r}
# Coefficient table
create_coefficient_table(application_domain_robust)
```

### Comparison of Effect Sizes Across AI Application Domains (Reference Category: Operations & Management)
```{r}
# Coefficient table
create_coefficient_table(operations_management_robust)
```

## Temporal AI Advancement
```{r}
temporal_evolution <- metafor::rma.mv(cohens_d, variance_d, mods = ~ scale(year) - 1, data = df, random = ~ 1 | article_id/es_id, tdist = TRUE)
temporal_evolution_robust <- robust(temporal_evolution, cluster = study_id, adjust=TRUE)

# Temporal evolution plot.
events <- data.frame(
  year = c(2011, 2015, 2018, 2020, 2021, 2022),
  event = c("IBM Watson wins Jeopardy",
            "AlphaGo defeats World Champion", 
            "BERT Release",
            "GPT-3 Launch", 
            "DALL-E Launch", 
            "ChatGPT Launch"))

events <- events %>%
  arrange(year) %>%
  group_by(year) %>%
  mutate(y_pos = if(n() > 1) c(0.75, 0.5) else 0.75) %>%
  ungroup()

my_theme <- theme(
  text = element_text(family = "Arial", size = 9),
  axis.text.x = element_text(angle = 45, hjust = 1, size = 9),
  axis.text.y = element_text(size = 9),
  axis.title.x = element_text(size = 9, margin = margin(t = 5)),
  axis.title.y = element_text(size = 9, margin = margin(r = 5)))

temporal_evolution_plot <- ggplot(df, aes(x = year, y = cohens_d)) +
  annotate("rect", 
           xmin = min(df$year), xmax = 2016, 
           ymin = -Inf, ymax = Inf, 
           fill = "#666666", alpha = 0.15) +
  annotate("rect", 
           xmin = 2016, xmax = max(df$year), 
           ymin = -Inf, ymax = Inf, 
           fill = "#4A90E2", alpha = 0.15) +
  annotate("text",
           x = mean(c(min(df$year), 2016)),
           y = 0.95,
           label = "Era of Rule-Based\nAlgorithmic Forms of AI",
           size = 3.5,
           color = "#666666",
           fontface = "bold",
           hjust = 0.5) +
  annotate("segment",
           x = min(df$year), xend = 2016,
           y = 0.80, yend = 0.80,
           color = "#666666",
           size = 0.5,
           arrow = arrow(ends = "both", length = unit(0.1, "cm"), type = "closed")) +
  annotate("text",
           x = mean(c(2016, max(df$year))),
           y = 0.95,
           label = "Era of LLMs and\nGenerative Forms of AI",
           size = 3.5,
           color = "#4A90E2",
           fontface = "bold",
           hjust = 0.5) +
  annotate("segment",
           x = 2016, xend = max(df$year),
           y = 0.80, yend = 0.80,
           color = "#4A90E2",
           size = 0.5,
           arrow = arrow(ends = "both", length = unit(0.1, "cm"), type = "closed")) +
  geom_segment(data = events,
               aes(x = year, xend = year,
                   y = -1, yend = 0.75),
               color = "gray20",
               linetype = "dashed",
               linewidth = 0.3) +
  geom_smooth(method = "loess", se = FALSE,
              color = "#0072B2") +
  geom_text(data = events,
            aes(x = year, y = y_pos, label = event),
            angle = 90,
            hjust = 1,
            vjust = -0.25,
            size = 3,
            color = "gray20",
            fontface = "bold") +
  labs(x = "Year of Publication",
       y = "Cohen's d") +
  theme_minimal() + 
  my_theme +
  theme(panel.grid.minor = element_blank()) +
  scale_x_continuous(
    breaks = seq(min(df$year), max(df$year), by = 2),
    expand = c(0.05, 0.05)
  ) +
  coord_cartesian(ylim = c(-1, 1))

############## MODERATOR ANALYSIS AI DEVELOPMENT PHASES ##############
# Effect size estimates for each ai development phase.
ai_development_phase <- metafor::rma.mv(cohens_d, variance_d, mods = ~ factor(ai_development_phase) - 1, data = df, random = ~ 1 | article_id/es_id, tdist = TRUE)
ai_development_phase_robust <- robust(ai_development_phase, cluster = study_id, adjust=TRUE)

# Comparison of effect sizes across ai development phases (reference category: before 2010). 
before_2010 <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(ai_development_phase), ref = 'before 2010'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:4)
before_2010_robust <- robust(before_2010, cluster = study_id)

# Comparison of effect sizes across ai development phases (reference category: 2010-2014). 
from_2010_2014 <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(ai_development_phase), ref = '2010-2014'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:4)
from_2010_2014_robust <- robust(from_2010_2014, cluster = study_id)

# Comparison of effect sizes across ai development phases (reference category: 2015-2019). 
from_2015_2019 <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(ai_development_phase), ref = '2015-2019'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:4)
from_2015_2019_robust <- robust(from_2015_2019, cluster = study_id)

# Comparison of effect sizes across ai development phases (reference category: 2020-2025). 
from_2020_2025 <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(ai_development_phase), ref = '2020-2025'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:4)
from_2020_2025_robust <- robust(from_2020_2025, cluster = study_id)

############## MODERATOR ANALYSIS AI DEVELOPMENT PHASES FOR COGNITIVE CONSUMER RESPONSES ##############
cognitive_data <- df %>% 
  filter(consumer_response == "cognitive")

ai_development_phase_cognitive <- metafor::rma.mv(cohens_d, variance_d, mods = ~ factor(ai_development_phase) - 1, data = cognitive_data, random = ~ 1 | article_id/es_id, tdist = TRUE)
ai_development_phase_cognitive_robust <- robust(ai_development_phase_cognitive, cluster = study_id, adjust=TRUE)

before_2010_cognitive <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(ai_development_phase), ref = 'before 2010'), data = cognitive_data, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:4)
before_2010_cognitive_robust <- robust(before_2010_cognitive, cluster = study_id)

from_2010_2014_cognitive <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(ai_development_phase), ref = "2010-2014"), data = cognitive_data, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:4)
from_2010_2014_cognitive_robust <- robust(from_2010_2014_cognitive, cluster = study_id)

from_2015_2019_cognitive <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(ai_development_phase), ref = '2015-2019'), data = cognitive_data, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:4)
from_2015_2019_cognitive_robust <- robust(from_2015_2019_cognitive, cluster = study_id)

from_2020_2025_cognitive <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(ai_development_phase), ref = '2020-2025'), data = cognitive_data, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:4)
from_2020_2025_cognitive_robust <- robust(from_2020_2025_cognitive, cluster = study_id)

############## MODERATOR ANALYSIS AI DEVELOPMENT PHASES FOR AFFECTIVE CONSUMER RESPONSES ##############
affective_data <- df %>% 
  filter(consumer_response == "affective")

ai_development_phase_affective <- metafor::rma.mv(cohens_d, variance_d, mods = ~ factor(ai_development_phase) - 1, data = affective_data, random = ~ 1 | article_id/es_id, tdist = TRUE)
ai_development_phase_affective_robust <- robust(ai_development_phase_affective, cluster = study_id, adjust=TRUE)

before_2010_affective <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(ai_development_phase), ref = 'before 2010'), data = affective_data, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:4)
before_2010_affective_robust <- robust(before_2010_affective, cluster = study_id)

from_2010_2014_affective <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(ai_development_phase), ref = "2010-2014"), data = affective_data, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:4)
from_2010_2014_affective_robust <- robust(from_2010_2014_affective, cluster = study_id)

from_2015_2019_affective <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(ai_development_phase), ref = '2015-2019'), data = affective_data, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:4)
from_2015_2019_affective_robust <- robust(from_2015_2019_affective, cluster = study_id)

from_2020_2025_affective <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(ai_development_phase), ref = '2020-2025'), data = affective_data, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:4)
from_2020_2025_affective_robust <- robust(from_2020_2025_affective, cluster = study_id)

############## MODERATOR ANALYSIS AI DEVELOPMENT PHASES FOR BEHAVIORAL CONSUMER RESPONSES ##############
behavioral_data <- df %>% 
  filter(consumer_response == "behavioral")

ai_development_phase_behavioral <- metafor::rma.mv(cohens_d, variance_d, mods = ~ factor(ai_development_phase) - 1, data = behavioral_data, random = ~ 1 | article_id/es_id, tdist = TRUE)
ai_development_phase_behavioral_robust <- robust(ai_development_phase_behavioral, cluster = study_id, adjust=TRUE)

before_2010_behavioral <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(ai_development_phase), ref = 'before 2010'), data = behavioral_data, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:4)
before_2010_behavioral_robust <- robust(before_2010_behavioral, cluster = study_id)

from_2015_2019_behavioral <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(ai_development_phase), ref = '2015-2019'), data = behavioral_data, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:4)
from_2015_2019_behavioral_robust <- robust(from_2015_2019_behavioral, cluster = study_id)

from_2020_2025_behavioral <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(ai_development_phase), ref = '2020-2025'), data = behavioral_data, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:4)
from_2020_2025_behavioral_robust <- robust(from_2020_2025_behavioral, cluster = study_id)

# Consumer responses over time plot.
prepare_data <- function(data, response_type, response_label) {
  data %>%
    filter(consumer_response == response_type) %>%
    mutate(ai_development_phase = case_when(
      ai_development_phase %in% c("before 2010", "2010-2014") ~ "before 2015",
      TRUE ~ as.character(ai_development_phase)
    )) %>%
    group_by(ai_development_phase) %>%
    summarise(
      mean_d = mean(cohens_d, na.rm = TRUE),
      se = sd(cohens_d, na.rm = TRUE) / sqrt(n()),
      response_label = response_label
    ) %>%
    mutate(ai_development_phase = factor(ai_development_phase, 
                                         levels = c("before 2015", "2015-2019", "2020-2025")))
}

cognitive_data <- prepare_data(df, "cognitive", "Cognitive")
affective_data <- prepare_data(df, "affective", "Affective")
behavioral_data <- prepare_data(df, "behavioral", "Behavioral")

create_response_plot <- function(data, response_name, show_y_axis = TRUE) {
  ggplot(data, aes(x = ai_development_phase, y = mean_d, group = response_label)) +
    geom_hline(yintercept = 0, linetype = "dashed", size = 1, color = "#a50f15") +
    geom_errorbar(aes(ymin = mean_d - se, ymax = mean_d + se), 
                  width = 0.2, 
                  color = "#08519c") +
    geom_line(color = "#08519c", linewidth = 1) +
    geom_point(color = "#08519c", fill = "#6baed6", shape = 21, size = 3) +
    theme +
    theme(
      axis.text.y = if (show_y_axis) element_text(size = 9) else element_blank(),
      axis.ticks.y = if (show_y_axis) element_line() else element_blank(),
      legend.position = "none",
      panel.grid.minor = element_blank(),
      panel.grid.major.y = element_line(color = "gray90"),
      panel.grid.major.x = element_blank(),
      strip.text = element_text(size = 9, face = "bold", color = "#08519c"),
      plot.title = element_text(size = 9, face = "bold", hjust = 0.5)
    ) +
    scale_y_continuous(limits = c(-0.9, 0.20), breaks = seq(-0.9, 0.20, by = 0.1)) +
    labs(x = "AI Development Phase", 
         y = if (show_y_axis) "Cohen's d" else NULL,
         title = response_name)
}

cognitive_plot <- create_response_plot(cognitive_data, "Cognitive", show_y_axis = TRUE)
affective_plot <- create_response_plot(affective_data, "Affective", show_y_axis = FALSE)
behavioral_plot <- create_response_plot(behavioral_data, "Behavioral", show_y_axis = FALSE)

responses_over_time_plot <- cognitive_plot | affective_plot | behavioral_plot
```

### Variation in Effect Sizes Over Time
```{r}
temporal_evolution_plot
```

### Effect Size Estimate for Temporal Evolution
```{r}
# Coefficient table
create_coefficient_table(temporal_evolution_robust)
```

### Variation in Effect Sizes Across Consumer Responses Over Time
```{r}
responses_over_time_plot
```

### Comparison of Effect Sizes Across AI Development Phases
```{r}
# Coefficient table
create_coefficient_table(ai_development_phase_robust)
```

### Comparison of Effect Sizes Across Consumer Responses and AI Development Phases

#### Cognitive Responses
```{r}
# Coefficient table
create_coefficient_table(ai_development_phase_cognitive_robust)
```

#### Affective Responses
```{r}
# Coefficient table
create_coefficient_table(ai_development_phase_affective_robust)
```

#### Behavioral Responses
```{r}
# Coefficient table
create_coefficient_table(ai_development_phase_behavioral_robust)
```

## Task and Study Characteristics
```{r}
############## MODERATOR ANALYSIS TASK OBJECTIVITY ##############
# Effect size estimates for each level of task objectivity.
task_objectivity <- metafor::rma.mv(cohens_d, variance_d, mods = ~ factor(task_objectivity) - 1, data = df, random = ~ 1 | article_id/es_id, tdist = TRUE)
task_objectivity_robust <- robust(task_objectivity, cluster = study_id, adjust=TRUE)

# Comparison of effect sizes across levels of task objectivity (reference category: objective). 
objective <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(task_objectivity), ref = 'objective'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:2)
objective_robust <- robust(objective, cluster = study_id)

# Task objectivity plot.
coef_order <- coef(task_objectivity_robust)
task_objectivity_plot <- df %>%
  filter(!is.na(task_objectivity)) %>%
  mutate(task_objectivity = fct_recode(fct_reorder(task_objectivity, 
                                                   coef_order[paste0("factor(task_objectivity)", task_objectivity)], 
                                                   .desc = TRUE),
                                       "Objective" = "objective",
                                       "Subjective" = "subjective")) %>%
  ggplot(aes(x = task_objectivity, y = cohens_d)) +
  geom_violin(trim = FALSE, fill = blue_light, color = blue_dark) +
  geom_boxplot(width = 0.1, fill = blue_medium, color = blue_dark, 
               outlier.shape = NA) +
  stat_summary(fun = median, geom = "point", shape = 21, size = 2,
               fill = "white", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", size = 1, 
             color = red_line) +
  scale_y_continuous(limits = c(-3, 2), breaks = seq(-3, 2, 1)) +
  labs(x = "Task\nObjectivity", y = "Cohen's d") +
  theme_classic() +
  theme

############## MODERATOR ANALYSIS TASK CONSEQUENTIALITY ##############
# Effect size estimates for each level of task consequentiality.
task_consequentiality <- metafor::rma.mv(cohens_d, variance_d, mods = ~ factor(task_consequentiality) - 1, data = df, random = ~ 1 | article_id/es_id, tdist = TRUE)
task_consequentiality_robust <- robust(task_consequentiality, cluster = study_id, adjust=TRUE)

# Comparison of effect sizes across levels of task consequentiality (reference category: low). 
low <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(task_consequentiality), ref = 'low'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:2)
low_robust <- robust(low, cluster = study_id)

# Task consequentiality plot.
coef_order <- coef(task_consequentiality_robust)
task_consequentiality_plot <- df %>%
  filter(!is.na(task_consequentiality)) %>%
  mutate(task_consequentiality = fct_recode(fct_reorder(task_consequentiality, 
                                                        coef_order[paste0("factor(task_consequentiality)", task_consequentiality)], 
                                                        .desc = TRUE),
                                            "High" = "high",
                                            "Low" = "low")) %>%
  ggplot(aes(x = task_consequentiality, y = cohens_d)) +
  geom_violin(trim = FALSE, fill = blue_light, color = blue_dark) +
  geom_boxplot(width = 0.1, fill = blue_medium, color = blue_dark, 
               outlier.shape = NA) +
  stat_summary(fun = median, geom = "point", shape = 21, size = 2,
               fill = "white", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", size = 1, 
             color = red_line) +
  scale_y_continuous(limits = c(-3, 2), breaks = seq(-3, 2, 1)) +
  labs(x = "Task\nConsequentiality", y = "Cohen's d") +
  theme_classic() +
  theme

############## MODERATOR ANALYSIS HUMAN BENCHMARK ##############
# Effect size estimates for each level of human expert.
human_benchmark <- metafor::rma.mv(cohens_d, variance_d, mods = ~ factor(human_benchmark) - 1, data = df, random = ~ 1 | article_id/es_id, tdist = TRUE)
human_benchmark_robust <- robust(human_benchmark, cluster = study_id, adjust=TRUE)

# Comparison of effect sizes across levels of human expert (reference category: non-expert). 
non_expert <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(human_benchmark), ref = 'non-expert'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:2)
non_expert_robust <- robust(non_expert, cluster = study_id)

# Human benchmark plot.
coef_order <- coef(human_benchmark_robust)
human_benchmark_plot <- df %>%
  filter(!is.na(human_benchmark)) %>%
  mutate(human_benchmark = fct_recode(fct_reorder(human_benchmark, 
                                                  coef_order[paste0("factor(human_benchmark)", human_benchmark)], 
                                                  .desc = TRUE),
                                      "Non-\nExpert" = "non-expert",
                                      "Expert" = "expert")) %>%
  ggplot(aes(x = human_benchmark, y = cohens_d)) +
  geom_violin(trim = FALSE, fill = blue_light, color = blue_dark) +
  geom_boxplot(width = 0.1, fill = blue_medium, color = blue_dark, 
               outlier.shape = NA) +
  stat_summary(fun = median, geom = "point", shape = 21, size = 2,
               fill = "white", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", size = 1, 
             color = red_line) +
  scale_y_continuous(limits = c(-3, 2), breaks = seq(-3, 2, 1)) +
  labs(x = "Human\nBenchmark", y = "Cohen's d") +
  theme_classic() +
  theme

############## MODERATOR ANALYSIS AI EXPLANATION ##############
# Effect size estimates for each level of AI explanation.
ai_explanation <- metafor::rma.mv(cohens_d, variance_d, mods = ~ factor(ai_explanation) - 1, data = df, random = ~ 1 | article_id/es_id, tdist = TRUE)
ai_explanation_robust <- robust(ai_explanation, cluster = study_id, adjust=TRUE)

# Comparison of effect sizes across levels of AI explanation (reference category: absent). 
absent_explanation <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(ai_explanation), ref = 'absent'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:2)
absent_explanation_robust <- robust(absent_explanation, cluster = study_id)

# AI explanation plot.
coef_order <- coef(ai_explanation_robust)
ai_explanation_plot <- df %>%
  filter(!is.na(ai_explanation)) %>%
  mutate(ai_explanation = fct_recode(fct_reorder(ai_explanation, 
                                                 coef_order[paste0("factor(ai_explanation)", ai_explanation)], 
                                                 .desc = TRUE),
                                     "Present" = "present",
                                     "Absent" = "absent")) %>%
  ggplot(aes(x = ai_explanation, y = cohens_d)) +
  geom_violin(trim = FALSE, fill = blue_light, color = blue_dark) +
  geom_boxplot(width = 0.1, fill = blue_medium, color = blue_dark, 
               outlier.shape = NA) +
  stat_summary(fun = median, geom = "point", shape = 21, size = 2,
               fill = "white", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", size = 1, 
             color = red_line) +
  scale_y_continuous(limits = c(-3, 2), breaks = seq(-3, 2, 1)) +
  labs(x = "AI\nExplanation", y = "Cohen's d") +
  theme_classic() +
  theme

############## MODERATOR ANALYSIS FIELD SETTING ##############
# Effect size estimates for each level of field setting.
field_setting <- metafor::rma.mv(cohens_d, variance_d, mods = ~ factor(field_setting) - 1, data = df, random = ~ 1 | article_id/es_id, tdist = TRUE)
field_setting_robust <- robust(field_setting, cluster = study_id, adjust=TRUE)

# Comparison of effect sizes across levels of field setting (reference category: no). 
no_field <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(field_setting), ref = 'no'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:2)
no_field_robust <- robust(no_field, cluster = study_id)

# Field setting plot.
coef_order <- coef(field_setting_robust)
field_setting_plot <- df %>%
  filter(!is.na(field_setting)) %>%
  mutate(field_setting = fct_recode(fct_reorder(field_setting, 
                                                coef_order[paste0("factor(field_setting)", field_setting)], 
                                                .desc = TRUE),
                                    "Yes" = "yes",
                                    "No" = "no")) %>%
  ggplot(aes(x = field_setting, y = cohens_d)) +
  geom_violin(trim = FALSE, fill = blue_light, color = blue_dark) +
  geom_boxplot(width = 0.1, fill = blue_medium, color = blue_dark, 
               outlier.shape = NA) +
  stat_summary(fun = median, geom = "point", shape = 21, size = 2,
               fill = "white", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", size = 1, 
             color = red_line) +
  scale_y_continuous(limits = c(-3, 2), breaks = seq(-3, 2, 1)) +
  labs(x = "Field\nSetting", y = "Cohen's d") +
  theme_classic() +
  theme

############## MODERATOR ANALYSIS STIMULUS PRESENTATION ##############
# Effect size estimates for each level of stimulus presentation.
stimulus_presentation <- metafor::rma.mv(cohens_d, variance_d, mods = ~ factor(stimulus_presentation) - 1, data = df, random = ~ 1 | article_id/es_id, tdist = TRUE)
stimulus_presentation_robust <- robust(stimulus_presentation, cluster = study_id, adjust=TRUE)

# Comparison of effect sizes across levels of stimulus presentation (reference category: text only). 
text_only <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(stimulus_presentation), ref = 'text only'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:2)
text_only_robust <- robust(text_only, cluster = study_id)

# Stimulus presentation plot.
coef_order <- coef(stimulus_presentation_robust)
stimulus_presentation_plot <- df %>%
  filter(!is.na(stimulus_presentation)) %>%
  mutate(stimulus_presentation = fct_recode(fct_reorder(stimulus_presentation, 
                                                        coef_order[paste0("factor(stimulus_presentation)", stimulus_presentation)], 
                                                        .desc = TRUE),
                                            "Percept.-\nRich" = "perceptually rich",
                                            "Text-\nOnly" = "text only")) %>%
  ggplot(aes(x = stimulus_presentation, y = cohens_d)) +
  geom_violin(trim = FALSE, fill = blue_light, color = blue_dark) +
  geom_boxplot(width = 0.1, fill = blue_medium, color = blue_dark, 
               outlier.shape = NA) +
  stat_summary(fun = median, geom = "point", shape = 21, size = 2,
               fill = "white", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", size = 1, 
             color = red_line) +
  scale_y_continuous(limits = c(-3, 2), breaks = seq(-3, 2, 1)) +
  labs(x = "Stimulus\nPresentation", y = "Cohen's d") +
  theme_classic() +
  theme

############## MODERATOR ANALYSIS INCENTIVE COMPATIBILITY ##############
# Effect size estimates for each level of incentive compatibility.
incentive_compatibility <- metafor::rma.mv(cohens_d, variance_d, mods = ~ factor(incentive_compatibility) - 1, data = df, random = ~ 1 | article_id/es_id, tdist = TRUE)
incentive_compatibility_robust <- robust(incentive_compatibility, cluster = study_id, adjust=TRUE)

# Comparison of effect sizes across levels of incentive compatibility  (reference category: no). 
no_ic <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(incentive_compatibility), ref = 'no'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:2)
no_ic_robust <- robust(no_ic, cluster = study_id)

# Incentive compatibility plot.
coef_order <- coef(incentive_compatibility_robust)
incentive_compatibility_plot <- df %>%
  filter(!is.na(incentive_compatibility)) %>%
  mutate(incentive_compatibility = fct_recode(fct_reorder(incentive_compatibility, 
                                                          coef_order[paste0("factor(incentive_compatibility)", incentive_compatibility)], 
                                                          .desc = TRUE),
                                              "Yes" = "yes",
                                              "No" = "no")) %>%
  ggplot(aes(x = incentive_compatibility, y = cohens_d)) +
  geom_violin(trim = FALSE, fill = blue_light, color = blue_dark) +
  geom_boxplot(width = 0.1, fill = blue_medium, color = blue_dark, 
               outlier.shape = NA) +
  stat_summary(fun = median, geom = "point", shape = 21, size = 2,
               fill = "white", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", size = 1, 
             color = red_line) +
  scale_y_continuous(limits = c(-3, 2), breaks = seq(-3, 2, 1)) +
  labs(x = "Incentive\nCompatibility", y = "Cohen's d") +
  theme_classic() +
  theme

############## MODERATOR ANALYSIS BEHAVIORAL DV ##############
# Effect size estimates for each level of behavioral dv.
behavioral_dv <- metafor::rma.mv(cohens_d, variance_d, mods = ~ factor(behavioral_dv) - 1, data = df, random = ~ 1 | article_id/es_id, tdist = TRUE)
behavioral_dv_robust <- robust(behavioral_dv, cluster = study_id, adjust=TRUE)

# Comparison of effect sizes across levels of behavioral dv  (reference category: no). 
no_bdv <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(behavioral_dv), ref = 'no'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:2)
no_bdv_robust <- robust(no_bdv, cluster = study_id)

# Behavioral DV plot.
coef_order <- coef(behavioral_dv_robust)
behavioral_dv_plot <- df %>%
  filter(!is.na(behavioral_dv)) %>%
  mutate(behavioral_dv = fct_recode(fct_reorder(behavioral_dv, 
                                                coef_order[paste0("factor(behavioral_dv)", behavioral_dv)], 
                                                .desc = TRUE),
                                    "Yes" = "yes",
                                    "No" = "no")) %>%
  ggplot(aes(x = behavioral_dv, y = cohens_d)) +
  geom_violin(trim = FALSE, fill = blue_light, color = blue_dark) +
  geom_boxplot(width = 0.1, fill = blue_medium, color = blue_dark, 
               outlier.shape = NA) +
  stat_summary(fun = median, geom = "point", shape = 21, size = 2,
               fill = "white", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", size = 1, 
             color = red_line) +
  scale_y_continuous(limits = c(-3, 2), breaks = seq(-3, 2, 1)) +
  labs(x = "Behavioral\nDV", y = "Cohen's d") +
  theme_classic() +
  theme

############## MODERATOR ANALYSIS WITHIN-SUBJECT DESIGN ##############
# Effect size estimates for each level of within-subject design.
within_subject_design <- metafor::rma.mv(cohens_d, variance_d, mods = ~ factor(within_subject_design) - 1, data = df, random = ~ 1 | article_id/es_id, tdist = TRUE)
within_subject_design_robust <- robust(within_subject_design, cluster = study_id, adjust=TRUE)

# Comparison of effect sizes across levels of within-subject design  (reference category: no). 
no_wsd <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(within_subject_design), ref = 'no'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:2)
no_wsd_robust <- robust(no_wsd, cluster = study_id)

# Within-subject design plot.
coef_order <- coef(within_subject_design_robust)
within_subject_design_plot <- df %>%
  filter(!is.na(within_subject_design)) %>%
  mutate(within_subject_design = fct_recode(fct_reorder(within_subject_design, 
                                                        coef_order[paste0("factor(within_subject_design)", within_subject_design)], 
                                                        .desc = TRUE),
                                            "Yes" = "yes",
                                            "No" = "no")) %>%
  ggplot(aes(x = within_subject_design, y = cohens_d)) +
  geom_violin(trim = FALSE, fill = blue_light, color = blue_dark) +
  geom_boxplot(width = 0.1, fill = blue_medium, color = blue_dark, 
               outlier.shape = NA) +
  stat_summary(fun = median, geom = "point", shape = 21, size = 2,
               fill = "white", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", size = 1, 
             color = red_line) +
  scale_y_continuous(limits = c(-3, 2), breaks = seq(-3, 2, 1)) +
  labs(x = "Within-Subject\nDesign", y = "Cohen's d") +
  theme_classic() +
  theme

############## MODERATOR ANALYSIS ONLINE SAMPLE ##############
# Effect size estimates for each level of online sample.
online_sample <- metafor::rma.mv(cohens_d, variance_d, mods = ~ factor(online_sample) - 1, data = df, random = ~ 1 | article_id/es_id, tdist = TRUE)
online_sample_robust <- robust(online_sample, cluster = study_id, adjust=TRUE)

# Comparison of effect sizes across levels of online sample  (reference category: no). 
no_os <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(online_sample), ref = 'no'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:2)
no_os_robust <- robust(no_os, cluster = study_id)

# Sample Type plot.
coef_order <- coef(online_sample_robust)
sample_type_plot <- df %>%
  filter(!is.na(online_sample)) %>%
  mutate(online_sample = fct_recode(fct_reorder(online_sample, 
                                                coef_order[paste0("factor(online_sample)", online_sample)], 
                                                .desc = TRUE),
                                    "Online" = "yes",
                                    "Other" = "no")) %>%
  ggplot(aes(x = online_sample, y = cohens_d)) +
  geom_violin(trim = FALSE, fill = blue_light, color = blue_dark) +
  geom_boxplot(width = 0.1, fill = blue_medium, color = blue_dark, 
               outlier.shape = NA) +
  stat_summary(fun = median, geom = "point", shape = 21, size = 2,
               fill = "white", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", size = 1, 
             color = red_line) +
  scale_y_continuous(limits = c(-3, 2), breaks = seq(-3, 2, 1)) +
  labs(x = "Sample\nType", y = "Cohen's d") +
  theme_classic() +
  theme

# Create plot grid
plots <- list(
  task_objectivity_plot,
  task_consequentiality_plot,
  human_benchmark_plot,
  ai_explanation_plot,
  field_setting_plot,
  stimulus_presentation_plot,
  incentive_compatibility_plot,
  behavioral_dv_plot,
  within_subject_design_plot,
  sample_type_plot
)

task_study_characteristics_grid <- plot_grid(plotlist = plots, nrow = 2, ncol = 5, align = "h",
                  labels = "AUTO", label_size = 10)
```

### Variation in Effect Sizes across Task & Study Characteristics
```{r, fig.height = 8, fig.width = 9}
task_study_characteristics_grid
```

### Effect Size Estimates for Each Level of Task Objectivity
```{r}
# Coefficient table
create_coefficient_table(task_objectivity_robust)
```

### Comparison of Effect Sizes Across Levels of Task Objectivity (Reference Category: Objective)
```{r}
# Coefficient table
create_coefficient_table(objective_robust)
```

### Effect Size Estimates for Each Level of Task Consequentiality
```{r}
# Coefficient table
create_coefficient_table(task_consequentiality_robust)
```

### Comparison of Effect Sizes Across Levels of Task Consequentiality (Reference Category: Low)
```{r}
# Coefficient table
create_coefficient_table(low_robust)
```

### Effect Size Estimates for Each Level of Human Benchmark
```{r}
# Coefficient table
create_coefficient_table(human_benchmark_robust)
```

### Comparison of Effect Sizes Across Levels of Human Benchmark (Reference Category: Non-Expert)
```{r}
# Coefficient table
create_coefficient_table(non_expert_robust)
```

### Effect Size Estimates for Each Level of AI Explanation
```{r}
# Coefficient table
create_coefficient_table(ai_explanation_robust)
```

### Comparison of Effect Sizes Across Levels of AI Explanation (Reference Category: Absent)
```{r}
# Coefficient table
create_coefficient_table(absent_explanation_robust)
```

### Effect Size Estimates for Each Level of Field Setting
```{r}
# Coefficient table
create_coefficient_table(field_setting_robust)
```

### Comparison of Effect Sizes Across Levels of Field Setting (Reference Category: No)
```{r}
# Coefficient table
create_coefficient_table(no_field_robust)
```

### Effect Size Estimates for Each Level of Stimulus Presentation
```{r}
# Coefficient table
create_coefficient_table(stimulus_presentation_robust)
```

### Comparison of Effect Sizes Across Levels of Stimulus Presentation (Reference Category: Text Only)
```{r}
# Coefficient table
create_coefficient_table(text_only_robust)
```

### Effect Size Estimates for Each Level of Incentive Compatibility
```{r}
# Coefficient table
create_coefficient_table(incentive_compatibility_robust)
```

### Comparison of Effect Sizes Across Levels of Incentive Compatibility (Reference Category: No)
```{r}
# Coefficient table
create_coefficient_table(no_ic_robust)
```

### Effect Size Estimates for Each Level of Behavioral DV
```{r}
# Coefficient table
create_coefficient_table(behavioral_dv_robust)
```

### Comparison of Effect Sizes Across Levels of Behavioral DV (Reference Category: No)
```{r}
# Coefficient table
create_coefficient_table(no_bdv_robust)
```

### Effect Size Estimates for Each Level of Within-Subject Design
```{r}
# Coefficient table
create_coefficient_table(within_subject_design_robust)
```

### Comparison of Effect Sizes Across Levels of Within-Subject Design (Reference Category: No)
```{r}
# Coefficient table
create_coefficient_table(no_wsd_robust)
```

### Effect Size Estimates for Each Level of Online Sample
```{r}
# Coefficient table
create_coefficient_table(online_sample_robust)
```

### Comparison of Effect Sizes Across Levels of Online Sample (Reference Category: No)
```{r}
# Coefficient table
create_coefficient_table(no_os_robust)
```

## Sample Characteristics
```{r}
############## MODERATOR ANALYSIS AVERAGE AGE ##############
average_age <- metafor::rma.mv(cohens_d, variance_d, mods = ~ scale(average_age) - 1, data = df, random = ~ 1 | article_id/es_id, tdist = TRUE)
average_age_robust <- robust(average_age, cluster = study_id, adjust=TRUE)

# Plot for average age.
average_age_plot <- ggplot(df, aes(x = average_age, y = cohens_d)) +
  geom_point(alpha = 0.2, color = "#08519c") +
  geom_smooth(method = "loess", color = "#08519c", fill = "#c6dbef") +
  theme_classic() +
  theme +
  geom_hline(yintercept = 0, linetype = "dashed", size = 1, color = "#a50f15") +
  scale_y_continuous(limits = c(-3, 2), breaks = seq(-3, 2, by = 1)) +
  labs(x = "Sample Age", y = "Cohen's d")

############## MODERATOR ANALYSIS PERCENTAGE FEMALES ##############
percentage_females <- metafor::rma.mv(cohens_d, variance_d, mods = ~ scale(percentage_females) - 1, data = df, random = ~ 1 | article_id/es_id, tdist = TRUE)
percentage_females_robust <- robust(percentage_females, cluster = study_id, adjust=TRUE)

# Plot for percentage of female participants.
percentage_females_plot <- ggplot(df, aes(x = percentage_females, y = cohens_d)) +
  geom_point(alpha = 0.2, color = "#08519c") +
  geom_smooth(method = "loess", color = "#08519c", fill = "#c6dbef") +
  theme_classic() +
  theme +
  geom_hline(yintercept = 0, linetype = "dashed", size = 1, color = "#a50f15") +
  scale_y_continuous(limits = c(-3, 2), breaks = seq(-3, 2, by = 1)) +
  scale_x_continuous(
    limits = c(0, 100),
    breaks = seq(0, 100, by = 20),
    labels = function(x) paste0(x, "%")
  ) +
  labs(x = "Female Participants", y = "Cohen's d")

############## MODERATOR ANALYSIS GEOGRAPHICAL BACKGROUND ##############
# Effect size estimates for each level of geographical background.
geographical_background <- metafor::rma.mv(cohens_d, variance_d, mods = ~ factor(geographical_background) - 1, data = df, random = ~ 1 | article_id/es_id, tdist = TRUE)
geographical_background_robust <- robust(geographical_background, cluster = study_id, adjust=TRUE)

# Comparison of effect sizes across levels of geographical background  (reference category: europe). 
europe <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(geographical_background), ref = 'europe'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:3)
europe_robust <- robust(europe, cluster = study_id)

# Comparison of effect sizes across levels of geographical background  (reference category: asia). 
asia <- metafor::rma.mv(cohens_d, variance_d, mods = ~ relevel(factor(geographical_background), ref = 'asia'), data = df, random = ~ 1 | article_id/es_id, tdist = TRUE, btt = 2:3)
asia_robust <- robust(asia, cluster = study_id)

# Plot for geographical background
coef_order <- coef(geographical_background_robust)
geographical_background_plot <- df %>%
  filter(!is.na(geographical_background)) %>%
  mutate(geographical_background = fct_recode(
    fct_reorder(geographical_background, 
                coef_order[paste0("factor(geographical_background)", geographical_background)], 
                .desc = TRUE),
    "America" = "america",
    "Asia" = "asia",
    "Europe" = "europe")) %>%
  ggplot(aes(x = geographical_background, y = cohens_d)) +
  geom_violin(trim = FALSE, fill = blue_light, color = blue_dark) +
  geom_boxplot(width = 0.1, fill = blue_medium, color = blue_dark, outlier.shape = NA) +
  stat_summary(fun = median, geom = "point", shape = 21, size = 2, fill = "white", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", size = 1, color = red_line) +
  scale_y_continuous(limits = c(-3, 2), breaks = seq(-3, 2, 1)) +
  labs(x = "Continent", y = "Cohen's d") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
geographical_background_plot

sample_grid <- plot_grid(average_age_plot, percentage_females_plot, geographical_background_plot, nrow = 1, align = "h", labels = "AUTO", label_size = 10)
```

### Variation in Effect Sizes across Sample Characteristics
```{r}
sample_grid
```

### Effect Size Estimate for Age
```{r}
# Coefficient table
create_coefficient_table(average_age_robust)
```

### Effect Size Estimate for Share of Female Participants
```{r}
# Coefficient table
create_coefficient_table(percentage_females_robust)
```

### Effect Size Estimates for Each Geographical Region
```{r}
# Coefficient table
create_coefficient_table(online_sample_robust)
```

### Comparison of Effect Sizes Across Geographical Regions (Reference Category: Europe)
```{r}
# Coefficient table
create_coefficient_table(europe_robust)
```

# AI LABELS X APPLICATION DOMAINS
```{r, fig.height = 6, fig.width = 8}
ai_label_levels <- c("ai systems", "ai assistants", "ai algorithms", "ai robots")
ai_label_labels <- c("AI Systems", "AI Assistants", "AI Algorithms", "AI Robots")

domain_totals <- df %>%
  filter(!is.na(application_domain)) %>%
  group_by(application_domain) %>%
  summarise(
    total_count = n(),
    avg_d = mean(cohens_d, na.rm = TRUE)
  ) %>%
  arrange(total_count) %>%
  mutate(domain_label = sprintf("%s\n(d = %.2f)", tools::toTitleCase(application_domain), avg_d))

ai_label_totals <- df %>%
  filter(!is.na(ai_label)) %>%
  group_by(ai_label) %>%
  summarise(
    total_count = n(),
    avg_d = mean(cohens_d, na.rm = TRUE)
  ) %>%
  mutate(
    ai_label = factor(ai_label, levels = ai_label_levels, labels = ai_label_labels),
    ai_label_label = sprintf("%s\n(d = %.2f)", ai_label, avg_d)
  )

plot_data <- df %>%
  filter(!is.na(application_domain), !is.na(ai_label)) %>%
  group_by(application_domain, ai_label) %>%
  summarise(
    count = n(),
    mean_cohens_d = mean(cohens_d, na.rm = TRUE),
    .groups = 'drop'
  ) %>%
  mutate(ai_label = factor(ai_label, levels = ai_label_levels, labels = ai_label_labels))

labelxdomain_matrix <- ggplot(plot_data, aes(
  x = ai_label,
  y = factor(application_domain, levels = domain_totals$application_domain)
)) +
  geom_point(aes(size = count, color = mean_cohens_d), alpha = 0.8) +
  geom_text(aes(label = sprintf("%.2f", mean_cohens_d)), color = "black", size = 2.5) +
  scale_size_continuous(range = c(3, 15), name = "Number of Effect Sizes") +
  scale_color_gradient2(
    low = "red", mid = "grey", high = "grey", midpoint = 0,
    limits = c(-1.6, 0.2), breaks = seq(-1.6, 0.2, by = 0.2),
    name = "Mean Cohen's d"
  ) +
  scale_x_discrete(labels = setNames(ai_label_totals$ai_label_label, ai_label_totals$ai_label)) +
  scale_y_discrete(labels = setNames(domain_totals$domain_label, domain_totals$application_domain)) +
  theme_minimal() +
  theme(
    text = element_text(family = "Arial", size = 10),
    axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
    axis.text.y = element_text(size = 10),
    axis.title.x = element_text(size = 10, margin = margin(t = 5)),
    axis.title.y = element_text(size = 10, margin = margin(r = 5)),
    legend.position = "right",
    legend.box = "vertical",
    legend.margin = margin(0, 0, 0, 0),
    legend.spacing = unit(0.3, "cm"),
    legend.key.size = unit(0.8, "cm"),
    legend.title = element_text(hjust = 0.5, size = 10),
    legend.text = element_text(size = 9),
    legend.box.just = "center",
    legend.justification = "center",
    plot.margin = unit(c(0.2, 0.2, 0.2, 0.2), "cm")
  ) +
  labs(x = "AI Label", y = "Application Domain") +
  guides(
    size = guide_legend(
      override.aes = list(color = "grey50"),
      order = 2,
      title.position = "top",
      label.position = "right",
      title.hjust = 0.5,
      label.hjust = 0.5,
      keywidth = 0.8,
      ncol = 1
    ),
    color = guide_colorbar(
      order = 1,
      barwidth = 0.8,
      barheight = 10,
      title.position = "top",
      label.position = "right",
      title.hjust = 0.5,
      label.hjust = 1
    )
  )
labelxdomain_matrix
```

# SESSION INFO {.appendix}

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







