pdf(NULL) # suppress Rplots.pdf in batch mode
# Parallelization refactor (2026-03-21): QTE bootstrap calls now run in parallel
# via mclapply across available cores, reducing runtime from ~70min to ~34min.
# Load Packages
library(tidyverse)        # Collection of R packages for data science (includes ggplot2, dplyr, tidyr, etc.)
library(xtable)
library(qte)
library(combinat)
library(gtsummary)        # Summary tables for regression models and more
library(stringr)          # String manipulation
library(broom)            # Convert statistical analysis objects into tidy tibbles
library(broom.mixed)      # Tidy tools for mixed models
library(here)             # Better way to manage working directories
library(tidylog, warn.conflicts = FALSE) # Logging for tidyverse operations
library(parallel)         # Parallel computation for bootstrap-heavy QTE calls

# Clean and Set up workspace
rm(list = ls())
set.seed(1324)

# Import and Apply Utils.R to Data
source(here("shared", "utils.R"))
user_iter_file <- here("similarity_exp_replication", "inputs",
                       "main_df_promptLvl_2026-03-22_PROD_HASHED_TRIMMED.csv")

result <- generate_data_from_user_iter_file(user_iter_file)
user_iter_df <- result$user_iter_df
user_df <- result$user_df

skill_vars <- c("Reading Comprehension", "Active Listening", "Writing", "Speaking", "Critical Thinking",                       
                "Social Perceptiveness", "Coordination", "Instructing", "Programming",
                "Judgment and Decision Making", "Systems Evaluation", "Science", "Active Learning",
                "Learning Strategies", "Monitoring", "Complex Problem Solving", "Operations Analysis",
                "Technology Design", "Troubleshooting", "Quality Control Analysis", "Systems Analysis")
group_vars <- c('GPT_usage4', 'degree3', 'gender', 'age5', 'GenAI_outlook3',
                'programming_skill3', 'program_frequency3', 'imagery_skill3', 'write_image_use3',
                'instruction_skill3', 'write_instr_usage3')
group_vars <- c(group_vars, skill_vars)

# remove any row with NA
user_df_comp <- user_df[complete.cases(user_df[,group_vars]),]

dep_vars <- c('mean_cosine_sim_mean', 'mean_dream_sim_mean',
              'z_mean_cosine_sim_mean', 'z_mean_dream_sim_mean',
              'pct_rank_cosine_sim_mean', 'z_mean_dream_sim_mean')


dalle_version_z_score_calculation <- function(user_df, outcome, outcome_clean) {
  
  user_df_3v <- user_df %>% filter(dalle_version == '3-verbatim-prompt')
  user_df_3r <- user_df %>% filter(dalle_version == '3-revised-prompt')
  user_df_2 <- user_df %>% filter(dalle_version == '2')
  
  mean_2 <- mean(user_df_2[[outcome]], na.rm = T)
  mean_3v <- mean(user_df_3v[[outcome]], na.rm=T)
  mean_3r <- mean(user_df_3r[[outcome]], na.rm=T)
  
  sd_2 <- sd(user_df_2[[outcome]], na.rm = T)
  sd_3v <- sd(user_df_3v[[outcome]], na.rm=T)
  sd_3r <- sd(user_df_3r[[outcome]], na.rm=T)
  
  n_2 <- sum(!is.na(user_df_2[[outcome]]))
  n_3v <- sum(!is.na(user_df_3v[[outcome]]))
  n_3r <- sum(!is.na(user_df_3r[[outcome]]))
  
  effect_23r <- (mean_3r - mean_2)
  effect_23v <- (mean_3v - mean_2)
  effect_3r3v <- (mean_3v - mean_3r)
  
  effect_se_23r <- sqrt(sd_2^2 / n_2 + sd_3r^2 / n_3r)
  effect_se_23v <- sqrt(sd_2^2 / n_2 + sd_3v^2 / n_3v)
  effect_se_3r3v <- sqrt(sd_3r^2 / n_3r + sd_3v^2 / n_3v)
  
  z23r <- effect_23r / effect_se_23r
  z23v <- effect_23v / effect_se_23v
  z3r3v <- effect_3r3v / effect_se_3r3v
  
  p23r <- 2 * pnorm(abs(z23r), lower.tail = FALSE)
  p23v <- 2 * pnorm(abs(z23v), lower.tail = FALSE)
  p3r3v <- 2 * pnorm(abs(z3r3v), lower.tail = FALSE)
  
  results_df <- data.frame(outcome = rep(outcome_clean, 3),
                           contrast = c('DALL-E 3 (revised) - DALL-E 2',
                                        'DALL-E 3 (verbatim) - DALL-E 2',
                                        'DALL-E 3 (verbatim) - DALL-E 3 (revised)'),
                           effects = c(effect_23r, effect_23v, effect_3r3v),
                           sds = c(effect_se_23r, effect_se_23v, effect_se_3r3v),
                           z = c(z23r, z23v, z3r3v),
                           p = c(p23r, p23v, p3r3v))
  return(results_df)
}

z_score_table <- rbind(dalle_version_z_score_calculation(user_df, 'z_mean_dream_sim_mean', 'Within-task z-score of 1-expected DreamSim'),
                       dalle_version_z_score_calculation(user_df, 'embedding_var', 'Prompt embedding variance'),
                       dalle_version_z_score_calculation(user_df, 'n_total_topics_cosine_tasklvl_mean',
                                                         'Number of topical transitions (cosine similarity)'),
                       dalle_version_z_score_calculation(user_df, 'n_total_topics_tsr_tasklvl_mean',
                                                         'Number of topical transitions (token sort ratio)'),
                       dalle_version_z_score_calculation(user_df, 'mean_cosine_sim_prev_prompt',
                                                         'Mean cosine similarity to previous prompt'),
                       dalle_version_z_score_calculation(user_df, 'mean_tsr_prev_prompt',
                                                         'Mean token sort ratio with previous prompt'),
                       dalle_version_z_score_calculation(user_df, 'mean_contains_prev_prompt',
                                                         'Mean percentage of time containing previous prompt'),
                       dalle_version_z_score_calculation(user_df, 'mean_prompt_length_chars',
                                                         'Mean prompt length (characters)'),
                       dalle_version_z_score_calculation(user_df, 'mean_prompt_length_words',
                                                         'Mean prompt length (words)')
)

rounded_z_score_table_3r3v <- z_score_table %>%
  filter(contrast == 'DALL-E 3 (verbatim) - DALL-E 3 (revised)') %>%
  dplyr::select(-contrast) %>%
  mutate(across(c(effects, sds, z, p), ~round(., 4)),
         across(c(outcome), as.character), # Ensure these are characters for manipulation
         outcome = ifelse(p <= 0.05, paste0('\\textbf{', outcome, '}'), outcome)) %>%
  rename(
    Outcome = outcome,
    Effect = effects,
    SE = sds,
    Z = z,
    p = p
  )

rounded_z_score_table_23r <- z_score_table %>%
  filter(contrast == 'DALL-E 3 (revised) - DALL-E 2') %>%
  dplyr::select(-contrast) %>%
  mutate(across(c(effects, sds, z, p), ~round(., 4)),
         across(c(outcome), as.character), # Ensure these are characters for manipulation
         outcome = ifelse(p <= 0.05, paste0('\\textbf{', outcome, '}'), outcome)) %>%
  rename(
    Outcome = outcome,
    Effect = effects,
    SE = sds,
    Z = z,
    p = p
  )

rounded_z_score_table_23v <- z_score_table %>%
  filter(contrast == 'DALL-E 3 (verbatim) - DALL-E 2') %>%
  dplyr::select(-contrast) %>%
  mutate(across(c(effects, sds, z, p), ~round(., 4)),
         across(c(outcome), as.character), # Ensure these are characters for manipulation
         outcome = ifelse(p <= 0.05, paste0('\\textbf{', outcome, '}'), outcome)) %>%
  rename(
    Outcome = outcome,
    Effect = effects,
    SE = sds,
    Z = z,
    p = p
  )


latex_table_23v <- xtable(rounded_z_score_table_23v,
                          align = rep("l", ncol(rounded_z_score_table_23v) + 1),
                          caption = "Difference between performance and prompting behavior outcomes for DALL-E 3 (verbatim) vs. DALL-E 2.")

dir.create(here("similarity_exp_replication", "outputs", "h5-h6"), recursive = TRUE, showWarnings = FALSE)
dir.create(here("similarity_exp_replication", "outputs", "tables", "pnas-h6"), recursive = TRUE, showWarnings = FALSE)

sink(here("similarity_exp_replication", "outputs", "h5-h6", "23v_performance_and_prompting_comp.tex"))
print(latex_table_23v, include.rownames = FALSE, sanitize.text.function = function(x) {
  x <- gsub("\\\textbf\\{", "\\textbf{", x) # Ensure correct escaping for LaTeX
  x
}, hline.after = c(-1, 0, nrow(latex_table_23v)))
sink()

# Convert to an xtable object
latex_table_3r3v <- xtable(rounded_z_score_table_3r3v,
                           align = rep("l", ncol(rounded_z_score_table_3r3v) + 1),
                           caption = "Difference between performance and prompting behavior outcomes for DALL-E 3 (verbatim) vs. DALL-E 3 (revised).")

sink(here("similarity_exp_replication", "outputs", "h5-h6", "3r3v_performance_and_prompting_comp.tex"))
print(latex_table_3r3v, include.rownames = FALSE, sanitize.text.function = function(x) {
  x <- gsub("\\\textbf\\{", "\\textbf{", x) # Ensure correct escaping for LaTeX
  x
}, hline.after = c(-1, 0, nrow(latex_table_3r3v)))
sink()

# Convert to an xtable object
latex_table_23r <- xtable(rounded_z_score_table_23r,
                          align = rep("l", ncol(rounded_z_score_table_23r) + 1),
                          caption = "Difference between performance and prompting behavior outcomes for DALL-E 3 (revised) vs. DALL-E 2.")

sink(here("similarity_exp_replication", "outputs", "h5-h6", "23r_performance_and_prompting_comp.tex"))
print(latex_table_23r, include.rownames = FALSE, sanitize.text.function = function(x) {
  x <- gsub("\\\textbf\\{", "\\textbf{", x) # Ensure correct escaping for LaTeX
  x
}, hline.after = c(-1, 0, nrow(latex_table_23r)))
sink()

long_df <- user_df %>%
  dplyr::select(dalle_version, z_mean_dream_sim_mean, embedding_var,
                n_total_topics_cosine_tasklvl_mean, n_total_topics_tsr_tasklvl_mean,
                mean_cosine_sim_prev_prompt, mean_tsr_prev_prompt, mean_contains_prev_prompt,
                mean_prompt_length_chars, mean_prompt_length_words) %>%
  pivot_longer(cols = -dalle_version,
               names_to = "column_name",
               values_to = "value")

# Now, calculate mean and sd for each column for each dalle_version
prompt_behavior_plot <- long_df %>%
  filter(column_name != 'z_mean_dream_sim_mean') %>%
  group_by(dalle_version, column_name) %>%
  summarise(mean = mean(value, na.rm = TRUE),
            sd = sd(value, na.rm = TRUE)/sqrt(sum(!is.na(value)))) %>%
  ungroup() %>%
  mutate(column_name = case_when(
    column_name == 'embedding_var' ~ 'Prompt embedding variance',
    column_name == 'mean_contains_prev_prompt' ~ 'Percentage of prompts containing previous prompt',
    column_name == 'mean_cosine_sim_prev_prompt' ~ 'Mean cosine similarity to previous prompt',
    column_name == 'mean_tsr_prev_prompt' ~ 'Mean token sort ratio with previous prompt',
    column_name == 'n_total_topics_cosine_tasklvl_mean' ~ 'Number of topical transitions (cosine similarity)',
    column_name == 'n_total_topics_tsr_tasklvl_mean' ~ 'Number of topical transitions (token sort ratio)',
    column_name == 'mean_prompt_length_words' ~ 'Mean prompt length (words)',
    TRUE ~ 'Mean prompt length (characters)'),
    dalle_version = case_when(
      dalle_version == '2' ~ 'DALL-E 2',
      dalle_version == '3-revised-prompt' ~ 'DALL-E 3 (revised)',
      TRUE ~ 'DALL-E 3 (verbatim)')
  ) %>%
  arrange(column_name) %>%
  ggplot(., aes(y=dalle_version, x=mean, color=dalle_version)) + geom_point() + facet_wrap(~column_name, nrow=4, scales='free') +
  geom_pointrange(aes(xmin=mean - 1.96*sd, xmax = mean + 1.96*sd, color=dalle_version)) +
  theme_minimal() +
  theme(legend.position = 'none',
        text = element_text(size = 11),
        strip.text.x.top = element_text(size = 10)) +
  labs(
    x = '',
    y = ''
  )
ggsave(prompt_behavior_plot, file = here("similarity_exp_replication", "outputs", "h5-h6", "prompt_behavior_plot.pdf"), height=11, width=12, units='in')

ate_bands <- long_df %>%
  filter(column_name == 'z_mean_dream_sim_mean') %>%
  group_by(dalle_version) %>%
  summarise(mean = mean(value, na.rm = TRUE),
            sd = sd(value, na.rm = TRUE)/sqrt(sum(!is.na(value))),
            lb = mean - 1.96*sd,
            ub = mean + 1.96*sd) %>%
  ungroup()

default_colors <- scales::hue_pal()(3)

ate_dalle_version <- user_df %>%
  ggplot(., aes(x=z_mean_dream_sim_mean, color=dalle_version)) + stat_ecdf() +
  geom_vline(data=ate_bands, aes(xintercept=mean, color=dalle_version)) +
  geom_rect(data=ate_bands, aes(xmin = lb, xmax = ub, ymin = -Inf, ymax = Inf, fill=dalle_version),
            alpha = 0.2, inherit.aes = FALSE) +
  scale_color_manual(values = c("2" = default_colors[1],
                                "3-revised-prompt" = default_colors[2],
                                "3-verbatim-prompt" = default_colors[3]),
                     labels = c("2" = "DALL-E 2",
                                "3-revised-prompt" = "DALL-E 3 (Revised)",
                                "3-verbatim-prompt" = "DALL-E 3 (Verbatim)")) +
  scale_fill_manual(values = c("2" = default_colors[1],
                               "3-revised-prompt" = default_colors[2],
                               "3-verbatim-prompt" = default_colors[3]),
                    labels = c("2" = "DALL-E 2",
                               "3-revised-prompt" = "DALL-E 3 (Revised)",
                               "3-verbatim-prompt" = "DALL-E 3 (Verbatim)")) +
  guides(color = guide_legend(title = "Model"),
         fill = guide_legend(title = "Model")) +
  theme_minimal() +
  labs(x='Within-task z-score of 1-expected DreamSim',
       y='Empirical CDF') +
  theme(legend.position = 'bottom',
        text = element_text(size = 11),
        strip.text.x.top = element_text(size = 10),
        legend.text = element_text(size = 9),
        legend.title = element_text(size = 9))

ggsave(ate_dalle_version, file = here("similarity_exp_replication", "outputs", "h5-h6", "ate_dalle_version.pdf"), height=5, width=8.5, units='in')

### Hypothesis 6

# Pre-compute filtered data frames once (used by all calculate_qte calls)
user_df_for_qte <- user_df %>%
  mutate(dalle_version_2 = ifelse(dalle_version == '2', 1, 0),
         dalle_version_3r = ifelse(dalle_version == '3-revised-prompt', 1, 0),
         dalle_version_3v = ifelse(dalle_version == '3-verbatim-prompt', 1, 0))
user_df_for_qte_no_2 <- user_df_for_qte %>% filter(dalle_version_2 == 0)
user_df_for_qte_no_3r <- user_df_for_qte %>% filter(dalle_version_3r == 0)
user_df_for_qte_no_3v <- user_df_for_qte %>% filter(dalle_version_3v == 0)

calculate_qte <- function(data_no_2, data_no_3r, data_no_3v, outcome_var, outcome_var_pretty) {
  prob_vector <- c(.1, .2, .3, .4, .5, .6, .7, .8, .9)
  n_iters <- 10000

  # Running QTE calculations
  qte_3rv <- ci.qte(reformulate("dalle_version_3v", response = outcome_var), x=NULL, data=data_no_2,
                    se=TRUE, probs = prob_vector, iters=n_iters)

  qte_23v <- ci.qte(reformulate("dalle_version_3v", response = outcome_var), x=NULL, data=data_no_3r,
                    se=TRUE, probs = prob_vector, iters=n_iters)

  qte_23r <- ci.qte(reformulate("dalle_version_3r", response = outcome_var), x=NULL, data=data_no_3v,
                    se=TRUE, probs = prob_vector, iters=n_iters)

  # Combining results
  qte_df <- data.frame(
    qte = c(qte_3rv$qte, qte_23v$qte, qte_23r$qte),
    qte_lb = c(qte_3rv$qte.lower, qte_23v$qte.lower, qte_23r$qte.lower),
    qte_ub = c(qte_3rv$qte.upper, qte_23v$qte.upper, qte_23r$qte.upper),
    prob = c(qte_3rv$probs, qte_23v$probs, qte_23r$probs),
    contrast = c(rep('DALL-E 3 (verbatim) - DALL-E 3 (revised)', length(qte_3rv$qte)),
                 rep('DALL-E 3 (verbatim) - DALL-E 2', length(qte_23v$qte)),
                 rep('DALL-E 3 (revised) - DALL-E 2', length(qte_23r$qte))
    ),
    outcome = rep(outcome_var_pretty, length(qte_3rv$qte*3)),
    ate = c(rep(qte_3rv$ate, length(qte_3rv$qte)),
            rep(qte_23v$ate, length(qte_23v$qte)),
            rep(qte_23r$ate, length(qte_23r$qte))
    )
  )

  return(qte_df)
}

qte_outcomes <- list(
  list(var = 'embedding_var',                       pretty = 'Prompt embedding variance'),
  list(var = 'mean_contains_prev_prompt',           pretty = 'Percentage of prompts containing previous prompt'),
  list(var = 'mean_cosine_sim_prev_prompt',         pretty = 'Mean cosine similarity to previous prompt'),
  list(var = 'mean_tsr_prev_prompt',                pretty = 'Mean token sort ratio with previous prompt'),
  list(var = 'n_total_topics_cosine_tasklvl_mean',  pretty = 'Number of topical transitions (cosine similarity)'),
  list(var = 'n_total_topics_tsr_tasklvl_mean',     pretty = 'Number of topical transitions (token sort ratio)'),
  list(var = 'mean_prompt_length_words',            pretty = 'Mean prompt length (words)'),
  list(var = 'z_mean_dream_sim_mean',               pretty = 'Task-Iter z-score of DreamSim'),
  list(var = 'z_mean_cosine_sim_mean',              pretty = 'Task-Iter z-score of CosineSim'),
  list(var = 'pct_rank_mean_cosine_sim_mean',       pretty = 'Avg Pctile Rank of CosineSim'),
  list(var = 'pct_rank_mean_dream_sim_mean',        pretty = 'Avg Pctile Rank of DreamSim'),
  list(var = 'mean_dream_sim_mean',                 pretty = 'Avg Raw DreamSim'),
  list(var = 'mean_cosine_sim_mean',                pretty = 'Avg Raw CosineSim')
)

n_cores <- min(detectCores() - 2, length(qte_outcomes))
cat(sprintf("Running %d QTE outcome vars in parallel across %d cores\n", length(qte_outcomes), n_cores))

qte_results <- mclapply(qte_outcomes, function(o) {
  calculate_qte(user_df_for_qte_no_2, user_df_for_qte_no_3r, user_df_for_qte_no_3v, o$var, o$pretty)
}, mc.cores = n_cores, mc.set.seed = TRUE)

qte_table <- do.call(rbind, qte_results)


hline_data <- qte_table %>%
  filter(outcome != 'Mean percentage of time containing previous prompt') %>%
  filter(outcome != 'Number of topical transitions (token sort ratio)') %>%
  filter(outcome != 'Number of topical transitions (cosine similarity)') %>%
  group_by(outcome, contrast) %>%
  summarise(ate = mean(ate, na.rm=T)) %>%
  ungroup()

qte_figure <- qte_table %>%
  filter(outcome != 'Mean percentage of time containing previous prompt') %>%
  filter(outcome != 'Number of topical transitions (token sort ratio)') %>%
  filter(outcome != 'Number of topical transitions (cosine similarity)') %>%
  mutate(prob = case_when(contrast == 'DALL-E 3 (verbatim) - DALL-E 3 (revised)' ~ prob + .02,
                          contrast == 'DALL-E 3 (verbatim) - DALL-E 2' ~ prob - .02,
                          TRUE ~ prob)) %>%
  ggplot(., aes(x=prob, y=qte, color=contrast)) + geom_point() +
  geom_pointrange(aes(ymin=qte_lb, ymax=qte_ub, color=contrast)) +
  theme_minimal() +
  theme(legend.position = 'bottom',
        text = element_text(size = 11),
        strip.text.x.top = element_text(size = 10),
        legend.text = element_text(size = 9),
        legend.title = element_text(size = 9)) +
  labs(
    color = 'Contrast',
    x = 'Quantile',
    y = 'Treatment Effect'
  ) +
  facet_wrap(~outcome, nrow=6, scales='free') +
  geom_hline(data = hline_data, aes(yintercept = ate, color=contrast), alpha=.3, linetype = "dashed") +
  geom_hline(aes(yintercept=0)) +
  scale_x_continuous(breaks = seq(0.1, 0.9, by = 0.1)) # This line adjusts the x-axis ticks
ggsave(qte_figure, file = here("similarity_exp_replication", "outputs", "tables", "pnas-h6", "qte_figure.pdf"), height=11, width=8.5, units='in')




perform_var_tests_df <- function(data, outcome_vars, group_var) {
  # Find all unique pairs of group_var (dalle_version)
  version_pairs <- combn(unique(data[[group_var]]), 2, simplify = FALSE)
  
  # Initialize an empty data frame to store results
  results_df <- data.frame()
  
  # Iterate over each outcome variable
  for(outcome_var in outcome_vars) {
    # Iterate over each pair of group_var (dalle_version) values
    for(pair in version_pairs) {
      subset_1 <- data %>% filter(.[[group_var]] == pair[1])
      subset_2 <- data %>% filter(.[[group_var]] == pair[2])
      
      # Combine subsets for var.test
      combined_data <- rbind(subset_1, subset_2)
      
      # Perform var.test
      test_result <- var.test(combined_data[[outcome_var]] ~ combined_data[[group_var]])
      
      # Extract relevant values from test_result
      temp_df <- data.frame(
        outcome_var = outcome_var,
        dalle_version_pair = paste(pair[1], pair[2], sep = " vs. "),
        f_statistic = test_result$statistic,
        variance_ratio = test_result$estimate,
        conf_lower = test_result$conf.int[1],
        conf_upper = test_result$conf.int[2],
        p_value = test_result$p.value
      )
      
      # Append temp_df to the master results_df
      results_df <- rbind(results_df, temp_df)
    }
  }
  
  return(results_df)
}

outcome_vars <- c('z_mean_cosine_sim_mean', 'pct_rank_mean_cosine_sim_mean', 'pct_rank_mean_dream_sim_mean',
                  'mean_dream_sim_mean', 'mean_cosine_sim_mean',
                  'z_mean_dream_sim_mean','embedding_var','n_total_topics_cosine_tasklvl_mean',
                  'n_total_topics_tsr_tasklvl_mean','mean_cosine_sim_prev_prompt','mean_tsr_prev_prompt',
                  'mean_contains_prev_prompt','mean_prompt_length_words')
f_tests <- perform_var_tests_df(user_df, outcome_vars, 'dalle_version')

vr_23v <- f_tests %>% mutate(outcome_var = case_when(
  outcome_var == 'embedding_var' ~ 'Prompt embedding variance',
  outcome_var == 'mean_contains_prev_prompt' ~ 'Percentage of prompts containing previous prompt',
  outcome_var == 'mean_cosine_sim_prev_prompt' ~ 'Mean cosine similarity to previous prompt',
  outcome_var == 'mean_tsr_prev_prompt' ~ 'Mean token sort ratio with previous prompt',
  outcome_var == 'n_total_topics_cosine_tasklvl_mean' ~ 'Number of topical transitions (cosine similarity)',
  outcome_var == 'n_total_topics_tsr_tasklvl_mean' ~ 'Number of topical transitions (token sort ratio)',
  outcome_var == 'mean_prompt_length_words' ~ 'Mean prompt length (words)',
  outcome_var == 'z_mean_dream_sim_mean' ~ 'Task-Iter z-score of DreamSim',
  outcome_var == 'z_mean_cosine_sim_mean' ~ 'Task-Iter z-score of CosineSim', 
  outcome_var == 'pct_rank_mean_cosine_sim_mean' ~ 'Avg Pctile Rank of CosineSim', 
  outcome_var == 'pct_rank_mean_dream_sim_mean' ~ 'Avg Pctile Rank of DreamSim',
  outcome_var == 'mean_dream_sim_mean' ~ 'Avg Raw DreamSim', 
  outcome_var == 'mean_cosine_sim_mean' ~ 'Avg Raw CosineSim',
  TRUE ~ 'Mean prompt length (characters)')) %>%
  filter(dalle_version_pair == '3-verbatim-prompt vs. 2') %>%
  dplyr::select(-dalle_version_pair) %>%
  dplyr::select(outcome_var, variance_ratio, conf_lower, conf_upper, f_statistic, p_value) %>%
  arrange(variance_ratio) %>%
  mutate(across(c(f_statistic, variance_ratio, conf_lower, conf_upper, p_value), ~round(., 4)),
         across(c(outcome_var), as.character), # Ensure these are characters for manipulation
         outcome_var = ifelse(p_value <= 0.05, paste0('\\textbf{', outcome_var, '}'), outcome_var)) %>%
  mutate(across(c(f_statistic, variance_ratio, conf_lower, conf_upper, p_value),
                ~if_else(.x < 0.0001, "$\\leq 10^{-4}$", as.character(.x)))) %>%
  rename(
    Outcome = outcome_var,
    `Variance ratio` = variance_ratio,
    `VR lower bound` = conf_lower,
    `VR upper bound` = conf_upper,
    `F stat` = f_statistic,
    `p-value` = p_value
  )

latex_table <- xtable(vr_23v,
                      align = c("l", "l", "r", "r", "r", "r", "r"),
                      caption = "Variance ratio and F stat for distributions of different outcomes under Dall-E 2 and DALL-E 3 (verbatim).")

# Print table without 'add.to.row'
sink(here("similarity_exp_replication", "outputs", "tables", "pnas-h6", "vr_23v.tex"))
print(latex_table, include.rownames = FALSE, sanitize.text.function = function(x) {
  x <- gsub("< 10^{-4}", "$< 10^{-4}$", x)
  x <- gsub("\\\textbf\\{", "\\textbf{", x) # Ensure correct escaping for LaTeX
  x
}, hline.after = c(-1, 0, nrow(vr_23v)),
size = "\\scriptsize")
sink()

vr_23r <- f_tests %>% mutate(outcome_var = case_when(
  outcome_var == 'embedding_var' ~ 'Prompt embedding variance',
  outcome_var == 'mean_contains_prev_prompt' ~ 'Percentage of prompts containing previous prompt',
  outcome_var == 'mean_cosine_sim_prev_prompt' ~ 'Mean cosine similarity to previous prompt',
  outcome_var == 'mean_tsr_prev_prompt' ~ 'Mean token sort ratio with previous prompt',
  outcome_var == 'n_total_topics_cosine_tasklvl_mean' ~ 'Number of topical transitions (cosine similarity)',
  outcome_var == 'n_total_topics_tsr_tasklvl_mean' ~ 'Number of topical transitions (token sort ratio)',
  outcome_var == 'mean_prompt_length_words' ~ 'Mean prompt length (words)',
  outcome_var == 'z_mean_dream_sim_mean' ~ 'Task-Iter z-score of DreamSim',
  outcome_var == 'z_mean_cosine_sim_mean' ~ 'Task-Iter z-score of CosineSim', 
  outcome_var == 'pct_rank_mean_cosine_sim_mean' ~ 'Avg Pctile Rank of CosineSim', 
  outcome_var == 'pct_rank_mean_dream_sim_mean' ~ 'Avg Pctile Rank of DreamSim',
  outcome_var == 'mean_dream_sim_mean' ~ 'Avg Raw DreamSim', 
  outcome_var == 'mean_cosine_sim_mean' ~ 'Avg Raw CosineSim',
  TRUE ~ 'Mean prompt length (characters)')) %>%
  filter(dalle_version_pair == '3-revised-prompt vs. 2') %>%
  dplyr::select(-dalle_version_pair) %>%
  dplyr::select(outcome_var, variance_ratio, conf_lower, conf_upper, f_statistic, p_value) %>%
  arrange(variance_ratio) %>%
  mutate(across(c(f_statistic, variance_ratio, conf_lower, conf_upper, p_value), ~round(., 4)),
         across(c(outcome_var), as.character), # Ensure these are characters for manipulation
         outcome_var = ifelse(p_value <= 0.05, paste0('\\textbf{', outcome_var, '}'), outcome_var)) %>%
  mutate(across(c(f_statistic, variance_ratio, conf_lower, conf_upper, p_value),
                ~if_else(.x <= 0.0001, "$\\leq 10^{-4}$", as.character(.x)))) %>%
  rename(
    Outcome = outcome_var,
    `Variance ratio` = variance_ratio,
    `VR lower bound` = conf_lower,
    `VR upper bound` = conf_upper,
    `F stat` = f_statistic,
    `p-value` = p_value
  )

latex_table <- xtable(vr_23r,
                      align = c("l", "l", "r", "r", "r", "r", "r"),
                      caption = "Variance ratio and F stat for distributions of different outcomes under Dall-E 2 and DALL-E 3 (revised).")

# Print table without 'add.to.row'
sink(here("similarity_exp_replication", "outputs", "tables", "pnas-h6", "vr_23r.tex"))
print(latex_table, include.rownames = FALSE, sanitize.text.function = function(x) {
  x <- gsub("< 10^{-4}", "$< 10^{-4}$", x)
  x <- gsub("\\\textbf\\{", "\\textbf{", x) # Ensure correct escaping for LaTeX
  x
}, hline.after = c(-1, 0, nrow(vr_23r)),
size = "\\scriptsize")
sink()

vr_3r3v <- f_tests %>% mutate(outcome_var = case_when(
  outcome_var == 'embedding_var' ~ 'Prompt embedding variance',
  outcome_var == 'mean_contains_prev_prompt' ~ 'Percentage of prompts containing previous prompt',
  outcome_var == 'mean_cosine_sim_prev_prompt' ~ 'Mean cosine similarity to previous prompt',
  outcome_var == 'mean_tsr_prev_prompt' ~ 'Mean token sort ratio with previous prompt',
  outcome_var == 'n_total_topics_cosine_tasklvl_mean' ~ 'Number of topical transitions (cosine similarity)',
  outcome_var == 'n_total_topics_tsr_tasklvl_mean' ~ 'Number of topical transitions (token sort ratio)',
  outcome_var == 'mean_prompt_length_words' ~ 'Mean prompt length (words)',
  outcome_var == 'z_mean_dream_sim_mean' ~ 'Task-Iter z-score of DreamSim',
  outcome_var == 'z_mean_cosine_sim_mean' ~ 'Task-Iter z-score of CosineSim', 
  outcome_var == 'pct_rank_mean_cosine_sim_mean' ~ 'Avg Pctile Rank of CosineSim', 
  outcome_var == 'pct_rank_mean_dream_sim_mean' ~ 'Avg Pctile Rank of DreamSim',
  outcome_var == 'mean_dream_sim_mean' ~ 'Avg Raw DreamSim', 
  outcome_var == 'mean_cosine_sim_mean' ~ 'Avg Raw CosineSim',
  TRUE ~ 'Mean prompt length (characters)')) %>%
  filter(dalle_version_pair == '3-verbatim-prompt vs. 3-revised-prompt') %>%
  dplyr::select(-dalle_version_pair) %>%
  dplyr::select(outcome_var, variance_ratio, conf_lower, conf_upper, f_statistic, p_value) %>%
  arrange(variance_ratio) %>%
  mutate(across(c(f_statistic, variance_ratio, conf_lower, conf_upper, p_value), ~round(., 4)),
         across(c(outcome_var), as.character), # Ensure these are characters for manipulation
         outcome_var = ifelse(p_value <= 0.05, paste0('\\textbf{', outcome_var, '}'), outcome_var)) %>%
  mutate(across(c(f_statistic, variance_ratio, conf_lower, conf_upper, p_value),
                ~if_else(.x <= 0.0001, "$\\leq 10^{-4}$", sprintf("%.4f", .x)))) %>%
  rename(
    Outcome = outcome_var,
    `Variance ratio` = variance_ratio,
    `VR lower bound` = conf_lower,
    `VR upper bound` = conf_upper,
    `F stat` = f_statistic,
    `p-value` = p_value
  )

latex_table <- xtable(vr_3r3v,
                      align = c("l", "l", "r", "r", "r", "r", "r"),
                      caption = "Variance ratio and F stat for distributions of different outcomes under Dall-E 3 (verbatim) and DALL-E 3 (revised).")

# Print table without 'add.to.row'
sink(here("similarity_exp_replication", "outputs", "tables", "pnas-h6", "vr_3r3v.tex"))
print(latex_table, include.rownames = FALSE, sanitize.text.function = function(x) {
  x <- gsub("< 10^{-4}", "$< 10^{-4}$", x)
  x <- gsub("\\\textbf\\{", "\\textbf{", x) # Ensure correct escaping for LaTeX
  x
}, hline.after = c(-1, 0, nrow(vr_3r3v)),
size = "\\scriptsize")
sink()

