---
title: "logistic_models"
author: "Xuefeiyang Zhang"
date: "2025-08-15"
output: html_document
---
# Logistic Regression Analysis 

library(sandwich)
library(lmtest)

# Read data
data <- read.delim("../logistic_model/modeling_dataset.tsv", sep = "\t")
data$Speaker <- as.factor(data$Speaker)

print(paste("Sample size:", nrow(data)))
print(paste("Speakers:", length(unique(data$Speaker))))

# Fit models
m1 <- glm(AnyError ~ zHNR + zTilt + zZCR + Label, data = data, family = binomial)
m2 <- glm(DelOccur ~ zHNR + zTilt + zZCR + Label, data = data, family = binomial)
m3 <- glm(InsOccur ~ zHNR + zTilt + zZCR + Label, data = data, family = binomial)
m4 <- glm(SubOccur ~ zHNR + zTilt + zZCR + Label, data = data, family = binomial)

# Get robust results for each model
get_results <- function(model, data) {
  robust_se <- vcovCL(model, cluster = data$Speaker)
  robust_test <- coeftest(model, vcov = robust_se)
  
  coefs <- robust_test[, "Estimate"]
  ses <- robust_test[, "Std. Error"]
  pvals <- robust_test[, "Pr(>|z|)"]
  
  ors <- exp(coefs)
  ci_low <- exp(coefs - 1.96 * ses)
  ci_high <- exp(coefs + 1.96 * ses)
  
  results <- data.frame(
    Variable = names(coefs),
    OR = round(ors, 2),
    CI_Low = round(ci_low, 2),
    CI_High = round(ci_high, 2),
    p = round(pvals, 3)
  )
  
  results$p[results$p < 0.001] <- "<0.001"
  return(results)
}

# Get results
r1 <- get_results(m1, data)
r2 <- get_results(m2, data)
r3 <- get_results(m3, data)
r4 <- get_results(m4, data)

# Remove intercept
r1 <- r1[r1$Variable != "(Intercept)", ]
r2 <- r2[r2$Variable != "(Intercept)", ]
r3 <- r3[r3$Variable != "(Intercept)", ]
r4 <- r4[r4$Variable != "(Intercept)", ]

# Make table
variables <- r1$Variable
table_result <- data.frame(
  Variable = variables,
  M1_AnyError = paste0(r1$OR, " [", r1$CI_Low, ", ", r1$CI_High, "] (", r1$p, ")"),
  M2_DelOccur = paste0(r2$OR, " [", r2$CI_Low, ", ", r2$CI_High, "] (", r2$p, ")"),
  M3_InsOccur = paste0(r3$OR, " [", r3$CI_Low, ", ", r3$CI_High, "] (", r3$p, ")"),
  M4_SubOccur = paste0(r4$OR, " [", r4$CI_Low, ", ", r4$CI_High, "] (", r4$p, ")")
)

# Clean variable names
table_result$Variable <- gsub("zHNR", "z-HNR", table_result$Variable)
table_result$Variable <- gsub("zTilt", "z-Spectral Tilt", table_result$Variable)
table_result$Variable <- gsub("zZCR", "z-ZCR", table_result$Variable)
table_result$Variable <- gsub("Labelsigh", "Event: sigh", table_result$Variable)
table_result$Variable <- gsub("Labelthroat-clearing", "Event: throat-clearing", table_result$Variable)

# Show and save table
print(table_result)
write.table(table_result, "table6_logistic_main.tsv", sep = "\t", row.names = FALSE, quote = FALSE)

# Make plots
# Plot 1: Insertion by event type
new_data <- data.frame(
  zHNR = 0, zTilt = 0, zZCR = 0,
  Label = c("laughter", "sigh", "throat-clearing")
)
pred <- predict(m3, newdata = new_data, type = "response")
pred_pct <- pred * 100

png("figure8_insertion_by_label.png", width = 800, height = 600)
barplot(pred_pct, names.arg = new_data$Label,
        main = "Predicted Insertion Probability by Event Type",
        xlab = "Event Type", ylab = "Probability (%)",
        col = "steelblue")
dev.off()

# Plot 2: AnyError vs HNR
hnr_values <- seq(-2, 2, by = 0.1)
new_data2 <- data.frame(
  zHNR = hnr_values, zTilt = 0, zZCR = 0, Label = "laughter"
)
pred2 <- predict(m1, newdata = new_data2, type = "response")
pred2_pct <- pred2 * 100

png("figure9_anyerror_vs_hnr.png", width = 800, height = 600)
plot(hnr_values, pred2_pct, type = "l", lwd = 2, col = "steelblue",
     main = "Any Error Probability vs HNR",
     xlab = "z-HNR", ylab = "Probability (%)")
dev.off()

# Plot 3: AnyError vs ZCR
zcr_values <- seq(-2, 2, by = 0.1)
new_data3 <- data.frame(
  zHNR = 0, zTilt = 0, zZCR = zcr_values, Label = "laughter"
)
pred3 <- predict(m1, newdata = new_data3, type = "response")
pred3_pct <- pred3 * 100

png("figure10_anyerror_vs_zcr.png", width = 800, height = 600)
plot(zcr_values, pred3_pct, type = "l", lwd = 2, col = "steelblue",
     main = "Any Error Probability vs ZCR",
     xlab = "z-ZCR", ylab = "Probability (%)")
dev.off()

print("Files saved:")
print("- table6_logistic_main.tsv")
print("- figure8_insertion_by_label.png")
print("- figure9_anyerror_vs_hnr.png")
print("- figure10_anyerror_vs_zcr.png")
print("")
print("Note: Using cluster-robust SEs by Speaker")