stats.Rmd

---
title: "stats"
output: html_document
---

```{r setup}
if (!requireNamespace("pacman", quietly = TRUE)) install.packages("pacman")
pacman::p_load(tidyverse, here, lme4, lmerTest, stats, kableExtra)
```


# Load Data
```{r, echo = FALSE, results = 'hide', message = FALSE, warning = FALSE} 
# (i'm surpressing the output of this code chunk)

version <- "3.0"
role_to_keep <- "assistant"

# Read text stats, filter for 'assistant' role
textstats_df <- read_csv(here("metrics", paste0("v", version, "_text_stats.csv"))) %>%
  filter(role == role_to_keep) %>%
  rename(level = group)

# Read textdescriptives, filter for 'assistant' role
textdescriptives_df <- read_csv(here("metrics", paste0("v", version, "_textdescriptives.csv"))) %>%
  filter(role == role_to_keep) %>%
  rename(text_length = doc_length) %>%
  rename(level = group)

# Read surprisal, filter for 'assistant' role
surprisal_df <- read_csv(here("metrics", paste0("v", version, "_surprisal.csv"))) %>%
  filter(role == role_to_keep) %>%
  rename(level = group)
```

```{r}
# check that they are the same
all(textstats_df$id == textdescriptives_df$id)
all(textstats_df$id == surprisal_df$id)
```

## Function DEF for running mixed effects
```{r}
# define the function to fit mixed effects models
# formula is metric ~ group + (1 | id)

fit_lmer_models <- function(data, models, metrics, output_dir) {
  # ensure the output directory exists
  dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
  
  # empty data frame to store results
  results_df <- data.frame(
    model = character(), 
    metric = character(), 
    term = character(), 
    estimate = numeric(),
    std_error = numeric(),
    t_value = numeric(),
    p_value = numeric(),
    stringsAsFactors = FALSE
  )
  
  # loop through models and metrics to fit the models and extract statistics
  for (model_name in models) {
    for (metric in metrics) {
      model_data <- data %>% filter(model == model_name)
      formula <- as.formula(paste(metric, "~ level + (1 | id)"))
      model <- lmer(formula, data = model_data)
      model_summary <- summary(model)
      
      # extract all coefficients and their statistics
      coefficients <- model_summary$coefficients
      
      # loop through terms to store results
      for (term in rownames(coefficients)) {
        results_df <- rbind(results_df, data.frame(
          model = model_name, 
          metric = metric, 
          term = term, 
          estimate = coefficients[term, "Estimate"],
          std_error = coefficients[term, "Std. Error"],
          t_value = coefficients[term, "t value"],
          p_value = round(coefficients[term, "Pr(>|t|)"], 4)  # round p-value
        ))
      }
      
      # save the model summary to a text file
      file_name <- paste0(model_name, "_", metric, "_summary.txt")
      file_path <- file.path(output_dir, file_name)
      capture.output(model_summary, file = file_path)
    }
  }
  
  # return the results data frame
  return(results_df)
}
```


## Run Mixed Effects
```{r}
models <- c("meta-llama--Llama-3.1-8B-Instruct", 
            "google--gemma-3-12b-it", 
            "mistralai--Mistral-7B-Instruct-v0.3", 
            "Qwen--Qwen2.5-7B-Instruct")
```


### Text Stats Models (Readability Metrics)
```{r}
output_dir <- here("results", paste0("v", version), "mixed_effects", "readability")  
metrics <- c("fernandez_huerta", "szigriszt_pazos", "gutierrez_polini")
results_textstats <- fit_lmer_models(textstats_df, models, metrics, output_dir)

# print the p-values
print(results_textstats)
```

### Text Descriptives Models (Structural Features)
```{r}
output_dir <- here("results", paste0("v", version), "mixed_effects", "structural")
metrics <- c("text_length", "dependency_distance_mean")
results_textdescriptives <- fit_lmer_models(textdescriptives_df, models, metrics, output_dir)

# print the p-values
print(results_textdescriptives)
```


### Surprisal Models (Surprisal Metrics)
```{r}
output_dir <- here("results", paste0("v", version), "mixed_effects", "surprisal")
metrics <- c("surprisal_mean")
results_surprisal <- fit_lmer_models(surprisal_df, models, metrics, output_dir)

# print the p-values
print(results_surprisal)
```
## Adjusting P-values for Multiple Comparisons
```{r}
# combine all p-values data frames
all_results <- bind_rows(results_textstats, results_textdescriptives, results_surprisal)

# count total number of tests
total_tests <- nrow(all_results)

# apply bonferroni correction
all_results <- all_results %>%
  mutate(
    p_value_corrected = pmin(p_value * total_tests, 1),  # bonferroni correction
    p_value_corrected = round(p_value_corrected, 4),  # round to 4 decimals
    stars = case_when(  # add significance stars
      p_value_corrected < 0.001 ~ "***",
      p_value_corrected < 0.01  ~ "**",
      p_value_corrected < 0.05  ~ "*",
      TRUE ~ ""
    )
  )

# rename the uncorrected p-value column
all_results <- all_results %>%
  select(-p_value)

# print or save the corrected p-values
print(all_results)
```

## Order and Save results
Order the dataset in a meaningful way:
```{r}
# define the desired order for metrics and models
metric_order <- c("fernandez_huerta", "szigriszt_pazos", "gutierrez_polini", 
                  "text_length", "dependency_distance_mean", "surprisal_mean")

model_order <- c("meta-llama--Llama-3.1-8B-Instruct", 
                 "google--gemma-3-12b-it", 
                 "mistralai--Mistral-7B-Instruct-v0.3", 
                 "Qwen--Qwen2.5-7B-Instruct")

all_results$metric <- factor(all_results$metric, levels = metric_order)
all_results$model <- factor(all_results$model, levels = model_order)

# sort by metric, model
all_results <- all_results[order(all_results$metric, all_results$model), ]
```


```{r}
all_results
```


```{r}
# save the corrected p-values to a CSV file
output_dir <- here("results", paste0("v", version), "mixed_effects")
output_file <- file.path(output_dir, "results.csv")
write_csv(all_results, output_file)
```

## Make LATEX table
```{r}
output_dir <- here("results", paste0("v", version), "mixed_effects")
```


```{r}
# format data
all_results <- all_results %>%
  mutate(
    formatted_metric = case_when(
      metric == "fernandez_huerta" ~ "Fernandez Huerta",
      metric == "szigriszt_pazos" ~ "Szigriszt Pazos",
      metric == "gutierrez_polini" ~ "Gutierrez Polini",
      metric == "text_length" ~ "Text Length",
      metric == "dependency_distance_mean" ~ "Mean Dependency Distance",
      metric == "surprisal_mean" ~ "Message Surprisal",
      TRUE ~ as.character(metric)
    ),
    # Clean up model names for display
    display_model = case_when(
      model == "meta-llama--Llama-3.1-8B-Instruct" ~ "Llama 3.1 8B Instruct",
      model == "google--gemma-3-12b-it" ~ "Gemma 3 12B IT",
      model == "mistralai--Mistral-7B-Instruct-v0.3" ~ "Mistral 7B Instruct v0.3",
      model == "Qwen--Qwen2.5-7B-Instruct" ~ "Qwen 2.5 7B Instruct",
      TRUE ~ as.character(model)
    )
  )

# use the formatted metric names for display
all_results <- all_results %>%
  arrange(metric, model, term) %>%
  group_by(metric) %>%
  mutate(
    # flag the first occurrence of each model within a metric group
    first_model_occurrence = !duplicated(model)
  ) %>%
  ungroup()

# create a column that shows the model name only on first occurrence
all_results <- all_results %>%
  mutate(
    display_model_selective = ifelse(first_model_occurrence, display_model, "")
  )
```


```{r}
# create base table
html_table <- all_results %>%
  # sort by metric first
  arrange(metric, model, term) %>%
  # select relevant columns
  select(display_model_selective, term, estimate, std_error, t_value, p_value_corrected, stars) %>%
  # create kable
  kable(
    format = "html", 
    booktabs = TRUE, 
    digits = 4,
    col.names = c("", "Term", "Estimate", "Std. Error", "t-value", "Corrected p-value", "Significance"), # model is first column, but is empty as it does not make sense with metric groupings
    align = "lcccccc", 
    escape = FALSE
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )

# group by metrics
metric_rows <- all_results %>%
  group_by(metric) %>%
  summarize(
    count = n(),
    formatted_name = first(formatted_metric)  # get the formatted name for each metric
  ) %>%
  ungroup()

# apply the grouping
row_index <- 1
for(i in 1:nrow(metric_rows)) {
  html_table <- html_table %>%
    pack_rows(
      group_label = metric_rows$formatted_name[i],
      start_row = row_index,
      end_row = row_index + metric_rows$count[i] - 1
    )
  
  row_index <- row_index + metric_rows$count[i]
}

# display the table
html_table
```


```{r}
html_table %>% save_kable(file = file.path(output_dir, "results_table.html"), self_contained = T)
```


```{r}
# First, create separate dataframes for the two tables
first_three_metrics <- c("fernandez_huerta", "szigriszt_pazos", "gutierrez_polini")
last_three_metrics <- c("text_length", "dependency_distance_mean", "surprisal_mean")
# Filter results for first table
first_table_results <- all_results %>%
  filter(metric %in% first_three_metrics)
# Filter results for second table
second_table_results <- all_results %>%
  filter(metric %in% last_three_metrics)
# Function to create formatted table with better model grouping
create_formatted_table <- function(data) {
  # Create base table
  formatted_table <- data %>%
    arrange(metric, model, term) %>%
    select(display_model_selective, term, estimate, std_error, t_value, p_value_corrected, stars) %>%
    kable(
      format = "latex", 
      booktabs = TRUE, 
      digits = 4,
      col.names = c("", "Term", "Est.", "SE", "t", "p (Adj.)", "Sig."),
      align = "llrrrrc",
      escape = FALSE
    ) %>%
    kable_styling(
      latex_options = c("scale_down", "hold_position"),
      font_size = 8,
    ) %>%
    column_spec(1, width = "12em") %>%
    column_spec(2, width = "6em") %>%
    column_spec(3:6, width = "5em")

  # Get unique metrics for grouping
  metrics <- unique(data$metric)

  # For each metric, apply spacing and grouping
  current_row <- 1
  for (current_metric in metrics) {
    # Filter data for the current metric
    metric_data <- data %>% filter(metric == current_metric)

    # Add the metric header first
    metric_name <- metric_data$formatted_metric[1]
    metric_count <- nrow(metric_data)

    formatted_table <- formatted_table %>%
      pack_rows(
        group_label = metric_name,
        start_row = current_row,
        end_row = current_row + metric_count - 1,
        bold = TRUE,
        italic = FALSE,
        hline_after = FALSE
      )

    # Now handle model grouping within each metric
    models <- unique(metric_data$model)
    model_row <- current_row

    for (current_model in models) {
      model_data <- metric_data %>% filter(model == current_model)
      model_name <- model_data$display_model_selective[1]
      model_count <- nrow(model_data)

      # If model appears multiple times with different terms, add subtle grouping
      if (model_count > 1) {
        formatted_table <- formatted_table %>%
          pack_rows(
            group_label = "", 
            start_row = model_row,
            end_row = model_row + model_count - 1,
            indent = TRUE,    # Indent the rows
            hline_after = FALSE
          )
      }

      model_row <- model_row + model_count
    }

    # Add line after this metric group
    formatted_table <- formatted_table %>%
      row_spec(current_row + metric_count - 1, extra_latex_after = "\\addlinespace")

    current_row <- current_row + metric_count
  }

  return(formatted_table)
}
# Create and save first table (readability metrics)
first_table <- create_formatted_table(first_table_results)
first_table %>% save_kable(file = file.path(output_dir, "readability_table.tex"), 
                          self_contained = TRUE)
# Create and save second table (structural metrics)
second_table <- create_formatted_table(second_table_results)
second_table %>% save_kable(file = file.path(output_dir, "structural_surprisal_table.tex"), 
                           self_contained = TRUE)
```