########################################
## Automated Scoring of Handwritten Mathematical Problem Posing Using Large Language Models
## Kar, Kisla, Kayici & Aydin
#####################################

# Read data
library(readxl)
mathas=data.frame(read_excel("MathASmainStudyScores.xlsx"))
# variables
## student id, task name, answer id, Human numeric score (see rubric in supp). 
## Original OCR reading (stdinput 1), revised OCR reading (expertrevision i2)
## AS is automated scoring, i is input 1 for original, 2 for revised.
## p is prompt, 1 for ascending scoring, 2 for descending scoring. 

# Descriptives 
## Load required packages
library(gtsummary)
library(dplyr)
library(tidyr)

## 1. Reshape the data so "Rater" becomes a grouping column
mathas_long <- mathas %>%
  select(Human, AS_i1_p1, AS_i1_p2, AS_i2_p1, AS_i2_p2) %>%
  pivot_longer(
    cols = everything(),
    names_to = "Rater_Type",
    values_to = "Score"
  ) %>%
  # Convert Score to a factor (0-6) so all score categories appear
  mutate(Score = factor(Score, levels = 0:6)) %>%
  # Convert Rater_Type to a factor to control the left-to-right column order and display labels
  mutate(Rater_Type = factor(Rater_Type, 
                             levels = c("Human", "AS_i1_p1", "AS_i1_p2", "AS_i2_p1", "AS_i2_p2"),
                             labels = c("Human Rater", "AS (i1, p1)", "AS (i1, p2)", "AS (i2, p1)", "AS (i2, p2)")))

## 2. Generate the side-by-side table
apa_table_side <- mathas_long %>%
  tbl_summary(
    by = Rater_Type, # This splits the raters into side-by-side columns
    statistic = all_categorical() ~ "{n} ({p}%)",
    missing_text = "Missing",
    label = list(Score ~ "Score Category")
  ) %>%
  # Format the column headers to be bold
  modify_header(all_stat_cols() ~ "**{level}**") %>%
  # Add an overarching header across all the rater columns
  modify_spanning_header(all_stat_cols() ~ "**Constructed Response Evaluation**") %>%
  bold_labels() %>%
  # Add the APA title
  modify_caption("**Table 1**\n*Descriptive Statistics for Human and Automated Scoring*")

## View the table
apa_table_side

## 1. Install the flextable package if you haven't already
## install.packages("flextable")
library(flextable)

## 2. Convert your existing gtsummary table to a flextable object
apa_table_word <- as_flex_table(apa_table_side)

## 3. Save it to your computer as a Word document
#save_as_docx(apa_table_word, path = "Table1_Descriptive_Statistics.docx")


# Descriptives done. Move to agreement
rm(list=setdiff(ls(),c("mathas")))

# Agreement metrics (Custom function named ASmetric)

## LWK, QWK manual, QWK by Metric, Pearson, Distance cor by Energy, Tertium Quid 1 and 2
library(irr)
library(Metrics) #does not work with integers
library(energy)

# Manual QWK function for two scores
qwktwo=function(a,b){
  vara=var(a)
  varb=var(b)
  meansqdiff=mean((a-b)^2)
  sqmeandiff=(mean(a)-mean(b))^2
  qwk=1-(meansqdiff)/(vara+varb+sqmeandiff)
  return(qwk)
}

ASmetric=function(a,b){
  lwk=kappa2(data.frame(a,b), weight = "equal")$value
  qwk=qwktwo(a,b)
  qwk2=Metrics::ScoreQuadraticWeightedKappa(a, 
                                   b,
                                   min.rating = 0, max.rating = 6)
  
  pear=as.numeric(cor(a,b))
  distr=as.numeric(energy::dcor(a,b))
  tq1p=sum(abs(a-b)>=1)/length(a)
  tq2p=sum(abs(a-b)>=2)/length(a)
  
  out=c(lwk,qwk,qwk2,pear,distr,tq1p,tq2p)
  return(out)
  }



# Load required packages
library(dplyr)
library(flextable)

# 1. Define the vectors to iterate over
as_columns <- c("AS_i1_p1", "AS_i1_p2", "AS_i2_p1", "AS_i2_p2")
task_types <- c("drive", "cake")

# Initialize an empty list to store the looped results
agreement_results <- list()

# 2. Iterate through each task type, and then through each AS scenario
for (t in task_types) {
  
  # Filter the data for the specific task
  task_data <- mathas %>% filter(Task == t)
  
  for (as_col in as_columns) {
    
    # Safely extract and convert factors back to numeric for the math operations
    human_scores <- as.numeric(as.character(task_data$Human))
    as_scores <- as.numeric(as.character(task_data[[as_col]]))
    
    # Run the custom metric function
    metrics <- ASmetric(human_scores, as_scores)
    
    # Structure the output and append to the list
    agreement_results[[paste(t, as_col)]] <- data.frame(
      # tools::toTitleCase capitalizes the first letter for a cleaner table output
      Task = tools::toTitleCase(t), 
      Scenario = as_col,
      LWK = metrics[1],
      QWK_Manual = metrics[2],
      QWK_Metrics = metrics[3],
      Pearson_r = metrics[4],
      Dist_Cor = metrics[5],
      TQ_1_Plus = metrics[6],
      TQ_2_Plus = metrics[7]
    )
  }
}

# 3. Bind the list into a single data frame and round numerics
agreement_df_grouped <- bind_rows(agreement_results) %>%
  mutate(across(where(is.numeric), ~round(.x, 3)))

# View the raw dataframe
print(agreement_df_grouped)


# 4. Convert the grouped dataframe to an APA-style Word table
grouped_table_word <- agreement_df_grouped %>%
  flextable() %>%
  # Merge identical consecutive values in the "Task" column vertically
  merge_v(j = "Task") %>% 
  # Align the merged Task labels to the top of their section
  valign(j = "Task", valign = "top") %>%
  # Add a horizontal line after the 4th row to visually separate Drive and Cake
  hline(i = 4, border = officer::fp_border(color = "black", width = 1)) %>%
  set_caption("**Table 2**\n*Agreement Metrics Between Human Raters and Automated Scoring by Task*") %>%
  autofit() %>%
  font(fontname = "Times New Roman", part = "all")

# 5. Export to Word
#save_as_docx(grouped_table_word, path = "Table2_Grouped_Agreement_Metrics.docx")



#analyze diff scores
rm(list=setdiff(ls(),c("mathas")))
head(mathas)
table(mathas$stdid)
length(unique(mathas$stdid))
#235 students

diffdata=mathas[,c(1,2,3,4,7:10)]
head(diffdata)
names(diffdata) <- c(
  "stuid", 
  "task", 
  "answer", 
  "human", 
  "inp_P1", 
  "inp_P2", 
  "rev_P1", 
  "rev_P2"
)
head(diffdata)

# wide to long
# Load required packages
library(tidyr)
library(dplyr)

# Transform from wide to long
diffdata_long <- diffdata %>%
  pivot_longer(
    cols = c(inp_P1, inp_P2, rev_P1, rev_P2),
    names_to = c("Input", "Prompt"),
    names_sep = "_",               # Splits "inp_P1" into "inp" and "P1"
    values_to = "ai_score"
  ) %>%
  mutate(
    # Make the Input names more descriptive
    Input = ifelse(Input == "inp", "Original", "Revised"),
    
    # Calculate our dependent variable (Human vs Automated difference)
    score_diff = abs(human - ai_score),
    
    # Convert variables to factors for the lmer model
    stuid = as.factor(stuid),
    task = as.factor(task),
    answer = as.factor(answer),    # Your constructed response items
    Input = as.factor(Input),
    Prompt = as.factor(Prompt)
  )


# Run the MLM
head(diffdata_long)
diffdata_long=data.frame(diffdata_long)
# Give each answer a uniqu id
diffdata_long$answer=paste0(diffdata_long$stuid,"_",diffdata_long$answer)

summary(diffdata_long$score_diff)
sd(diffdata_long$score_diff)


# use effect coding for interpretive purposes
options(contrasts = c("contr.sum", "contr.poly"))
library(nlme)
model_lme <- nlme::lme(score_diff ~ task * Input * Prompt, 
                       random = ~ 1 | stuid/answer, 
                       data = diffdata_long, 
                       method = "REML")
# View the detailed coefficients, t-values, and p-values
summary(model_lme)
VarCorr(model_lme)
# View the ANOVA table for the overall significance of main effects and interactions
anova(model_lme)

# Treat score_diff ordinal
# install.packages("ordinal")
# 1. ConvertID columns to explicit factors
diffdata_long2=diffdata_long
diffdata_long2$stuid <- as.factor(diffdata_long2$stuid)
diffdata_long2$answer <- as.factor(diffdata_long2$answer)

# 2. Create the nested random effect term explicitly 
# The interaction() function safely pastes the two factors together 
diffdata_long2$stuid_answer <- interaction(diffdata_long2$stuid, diffdata_long2$answer)
diffdata_long2$score_diff_ord <- factor(diffdata_long2$score_diff, ordered = TRUE)

# 3. Re-run the CLMM using the newly created variable
library(ordinal)
model_clmm <- clmm(score_diff_ord ~ task * Input * Prompt + 
                     (1 | stuid) + 
                     (1 | stuid_answer), 
                   data = diffdata_long2,
                   link = "logit")

summary(model_clmm)


# Try one more alternative to traditional MLM
library(glmmTMB)

model_nb <- glmmTMB(score_diff ~ task * Input * Prompt + (1 | stuid/answer), 
                    data = diffdata_long, 
                    family = nbinom2) # nbinom2 is standard negative binomial parametrization

summary(model_nb)


#Create table
# Load required packages
library(broom.mixed)
library(dplyr)
library(flextable)
library(officer)
library(nlme)

# 1. Extract and format ONLY the Fixed Effects
fixed_results <- tidy(model_lme, effects = "fixed") %>%
  mutate(
    # Format p-values to APA style
    p.value = case_when(
      is.na(p.value) ~ "",
      p.value < 0.001 ~ "< .001",
      TRUE ~ sub("^0", "", sprintf("%.3f", p.value))
    ),
    # Round estimates
    estimate = sprintf("%.2f", estimate),
    std.error = sprintf("%.2f", std.error),
    statistic = sprintf("%.2f", statistic),
    
    # Rename terms for effect coding interpretation
    Parameter = case_when(
      term == "(Intercept)" ~ "Intercept (Grand Mean)",
      term == "task1" ~ "Task (Level 1)",
      term == "Input1" ~ "Input (Level 1)",
      term == "Prompt1" ~ "Prompt (Level 1)",
      term == "task1:Input1" ~ "Task × Input",
      term == "task1:Prompt1" ~ "Task × Prompt",
      term == "Input1:Prompt1" ~ "Input × Prompt",
      term == "task1:Input1:Prompt1" ~ "Task × Input × Prompt",
      TRUE ~ term
    )
  ) %>%
  select(Parameter, Estimate = estimate, SE = std.error, t = statistic, p = p.value)

# 2. Extract Random Effects manually using VarCorr
# VarCorr returns a character matrix. Coercing to numeric turns the headers into NAs.
vc_matrix <- suppressWarnings(as.numeric(VarCorr(model_lme)[, "StdDev"]))
ran_pars <- vc_matrix[!is.na(vc_matrix)] # Drop the NAs to leave exactly the 3 StdDev values

# Create the Random Effects dataframe
random_results <- data.frame(
  Parameter = c("Student Variance (Std. Dev.)", 
                "Answer w/in Student (Std. Dev.)", 
                "Residual (Std. Dev.)"),
  Estimate = sprintf("%.2f", ran_pars),
  SE = "", t = "", p = "" # Blank columns to align with fixed effects
)

# 3. Combine Fixed and Random effects into one table
formatted_table <- bind_rows(fixed_results, random_results)

# 4. Create the APA-style flextable
lme_flextable <- formatted_table %>%
  flextable() %>%
  # Update caption to note the outcome variable
  set_caption("**Table 3**\n*Linear Mixed-Effects Model Predicting Absolute Difference Scores for Constructed Response Items*") %>%
  autofit() %>%
  font(fontname = "Times New Roman", part = "all") %>%
  align(j = 2:5, align = "center", part = "all") %>%
  # Horizontal line separating Fixed and Random effects
  hline(i = 8, border = fp_border(color = "black", width = 1)) %>%
  add_header_row(values = c("Fixed Effects", "", "", "", ""), top = FALSE) %>%
  compose(i = 10, j = 1, part = "body", value = as_paragraph(as_b("Random Effects")))

# 5. Export to Word
#save_as_docx(lme_flextable, path = "Table3_LME_Results_EffectCoding.docx")




# Remove interactions
options(contrasts = c("contr.treatment", "contr.poly"))
model_lme_simpler <- nlme::lme(score_diff ~ task + Input , 
                       random = ~ 1 | stuid/answer, 
                       data = diffdata_long, 
                       method = "REML")

summary(model_lme_simpler)
VarCorr(model_lme_simpler)

# Load required packages
library(broom.mixed)
library(dplyr)
library(flextable)
library(officer)
library(nlme)

# 1. Extract and format ONLY the Fixed Effects
fixed_results <- tidy(model_lme_simpler, effects = "fixed") %>%
  mutate(
    # Format p-values to APA style
    p.value = case_when(
      is.na(p.value) ~ "",
      p.value < 0.001 ~ "< .001",
      TRUE ~ sub("^0", "", sprintf("%.3f", p.value))
    ),
    # Round estimates
    estimate = sprintf("%.2f", estimate),
    std.error = sprintf("%.2f", std.error),
    statistic = sprintf("%.2f", statistic),
    
    # Rename terms for dummy coding interpretation
    Parameter = case_when(
      term == "(Intercept)" ~ "Intercept (Task: Cake, Input: Orig)",
      term == "taskdrive" ~ "Task: Drive",
      term == "InputRevised" ~ "Input: Revised",
      TRUE ~ term
    )
  ) %>%
  select(Parameter, Estimate = estimate, SE = std.error, t = statistic, p = p.value)

# 2. Extract Random Effects manually using VarCorr
vc_matrix <- suppressWarnings(as.numeric(VarCorr(model_lme_simpler)[, "StdDev"]))
ran_pars <- vc_matrix[!is.na(vc_matrix)]

# Create the Random Effects dataframe with a dedicated header row
random_results <- data.frame(
  Parameter = c("Random Effects", 
                "  Student Variance (Std. Dev.)", 
                "  Answer w/in Student (Std. Dev.)", 
                "  Residual (Std. Dev.)"),
  Estimate = c("", sprintf("%.2f", ran_pars)),
  SE = "", t = "", p = ""
)

# 3. Combine Fixed and Random effects
formatted_table <- bind_rows(fixed_results, random_results)

# 4. Create the APA-style flextable
lme_flextable_simpler <- formatted_table %>%
  flextable() %>%
  set_caption("**Table 4**\n*Main Effects Linear Mixed Model Predicting Absolute Difference Scores*") %>%
  autofit() %>%
  font(fontname = "Times New Roman", part = "all") %>%
  align(j = 2:5, align = "center", part = "all") %>%
  # Add a horizontal line after row 4 (the last fixed effect)
  hline(i = 4, border = fp_border(color = "black", width = 1)) %>%
  # Add the 'Fixed Effects' overarching column header
  add_header_row(values = c("Fixed Effects", "", "", "", ""), top = FALSE) %>%
  # Bold the "Random Effects" separator row in the body
  bold(i = 5, j = 1, part = "body")

# 5. Export to Word
#save_as_docx(lme_flextable_simpler, path = "Table4_Main_Effects_LME.docx")



# Remove interactions and standardize dependent var.
options(contrasts = c("contr.treatment", "contr.poly"))
diffdata_long$std_diff=scale(diffdata_long$score_diff)
model_lme_simplerstd <- nlme::lme(std_diff ~ task + Input , 
                               random = ~ 1 | stuid/answer, 
                               data = diffdata_long, 
                               method = "REML")

summary(model_lme_simplerstd)
VarCorr(model_lme_simplerstd)


# Load required packages
library(emmeans)
library(ggplot2)
library(dplyr)

# 1. Extract the predicted means and Confidence Intervals as a data frame
# Removed 'Prompt' from the formula so we just get Input | task
plot_data <- emmip(model_lme_simpler, ~ Input | task, 
                   CIs = TRUE, plotit = FALSE)

# 2. Build the customized APA-style plot
publication_plot <- ggplot(plot_data, aes(x = Input, y = yvar, group = 1)) +
  # Add lines and points (removed dodging and color/shape mappings)
  geom_line(linewidth = 0.8, color = "black") +
  geom_point(size = 3, color = "black") +
  # Add the confidence intervals (removed dodging)
  geom_errorbar(aes(ymin = LCL, ymax = UCL), 
                width = 0.15, 
                linewidth = 0.8, color = "black") +
  # Create panels for the Tasks and clean up the facet labels
  facet_wrap(~ task, labeller = as_labeller(c("cake" = "Task: Cake", "drive" = "Task: Drive"))) +
  # Add standard APA labels (removed legend titles)
  labs(
    title = "Figure X",
    subtitle = "Estimated Marginal Means of Absolute Difference Scores by Condition",
    x = "Input Type",
    y = "Absolute Difference Score"
  ) +
  # Apply a clean, minimal theme
  theme_classic() +
  theme(
    # Clean up the facet headers (strip)
    strip.background = element_rect(fill = "grey95", color = "black", linewidth = 1),
    strip.text = element_text(size = 11, face = "bold"),
    # Standardize font sizes for readability
    axis.text = element_text(size = 10, color = "black"),
    axis.title = element_text(size = 12, face = "bold"),
    # APA requires the word "Figure X" to be bolded, and the title to be italicized
    plot.title = element_text(size = 12, face = "bold", margin = margin(b = 5)),
    plot.subtitle = element_text(size = 12, face = "italic", margin = margin(b = 15))
  )

# 3. View the plot
print(publication_plot)

# 4. Save the plot in high resolution
# ggsave("Figure1_MainEffects.tiff", plot = publication_plot, 
#        width = 8, height = 5, units = "in", dpi = 300, compression = "lzw")


# additional descriptives
rm(list=setdiff(ls(),c("mathas","diffdata_long")))
head(diffdata_long)

# Load required package
library(dplyr)
mean(diffdata_long$score_diff)
# Replace 'your_data' with the actual name of your dataframe
descriptives_summary <- diffdata_long %>%
  # Group by your specified variables
  group_by(task, Input) %>%
  # Calculate the descriptives for score_diff
  summarise(
    n = n(),
    mean_diff = mean(score_diff, na.rm = TRUE),
    sd_diff = sd(score_diff, na.rm = TRUE),
    median_diff = median(score_diff, na.rm = TRUE),
    min_diff = min(score_diff, na.rm = TRUE),
    max_diff = max(score_diff, na.rm = TRUE),
    .groups = "drop" # Drops the grouping structure afterwards for a clean tibble
  )

# View the results
print(descriptives_summary)


descriptives_summary2 <- diffdata_long %>%
  # Group by your specified variables
  group_by( Input) %>%
  # Calculate the descriptives for score_diff
  summarise(
    n = n(),
    mean_diff = mean(score_diff, na.rm = TRUE),
    sd_diff = sd(score_diff, na.rm = TRUE),
    median_diff = median(score_diff, na.rm = TRUE),
    min_diff = min(score_diff, na.rm = TRUE),
    max_diff = max(score_diff, na.rm = TRUE),
    .groups = "drop" # Drops the grouping structure afterwards for a clean tibble
  )

# View the results
print(descriptives_summary2)