---
title: "descriptive_stat"
author: "Xuefeiyang Zhang"
date: "2025-08-12"
output: html_document
---

# Step 1: Event distribution analysis

# Set path
data_dir <- "../descriptive_stat"
output_dir <- "../descriptive_stat/output"

cat("=== DISTRIBUTION OF NUMBER OF EVENTS ===\n")

# Read data
data_file <- file.path(data_dir, "master_analysis_table.tsv")
data <- read.delim(data_file, sep = "\t", header = TRUE, stringsAsFactors = FALSE)

cat("Data loaded:", nrow(data), "observations\n")

# Create frequency table
label_counts <- table(data$Label)

# Convert to data frame and calculate percentages
event_dist <- data.frame(
  Label = names(label_counts),
  n = as.numeric(label_counts),
  stringsAsFactors = FALSE
)

# Calculate percentages
event_dist$percent <- paste0(round(event_dist$n / sum(event_dist$n) * 100, 1), "%")

# Sort by frequency
event_dist <- event_dist[order(event_dist$n, decreasing = TRUE), ]

# Add total row
total_row <- data.frame(
  Label = "Total",
  n = sum(event_dist$n),
  percent = "100%",
  stringsAsFactors = FALSE
)

# Combine table
table_event_distribution <- rbind(event_dist, total_row)

# Reset row names
rownames(table_event_distribution) <- NULL

# Print and save table
cat("\nEvent Distribution Table:\n")
print(table_event_distribution)
output_file <- file.path(output_dir, "table1_event_distribution.tsv")
write.table(table_event_distribution, output_file, sep = "\t", row.names = FALSE, quote = FALSE)

cat("\nTable saved to:", output_file, "\n")

# Create bar chart
cat("\nCreating bar chart...\n")

chart_data <- event_dist[event_dist$Label != "Total", ]

#save bar chart
png_file <- file.path(output_dir, "figure1_event_distribution.png")
png(png_file, width = 800, height = 600)

# Set bar chart
max_count <- max(chart_data$n)
y_limit <- max_count * 1.1 

barplot(chart_data$n, 
        names.arg = chart_data$Label,
        main = "Distribution of Paralinguistic Events",
        xlab = "Event Type",
        ylab = "Count",
        col = "steelblue",
        ylim = c(0,y_limit))
        
#Add text and labels
text(x = seq_along(chart_data$n) * 1.2 - 0.5, 
     y = chart_data$n + max(chart_data$n) * 0.02, 
     labels = chart_data$n, 
     pos = 3)

dev.off()

cat("Bar chart saved to:", png_file, "\n")

# Step 2: WER overview

# Set path
data_dir <- "../descriptive_stat"
output_dir <- "../descriptive_stat/output"

# Read data
data <- read.delim(file.path(data_dir, "master_analysis_table.tsv"), sep = "\t")
data$WER_Percent <- as.numeric(data$WER_Percent)
data <- data[!is.na(data$WER_Percent), ]

# Calculate WER statistics
event_types <- c("laughter", "sigh", "throat-clearing")
wer_stats <- data.frame()

for(label in event_types) {
  wer_values <- data[data$Label == label, "WER_Percent"]
  if(length(wer_values) > 0) {
    stats_row <- data.frame(
      Label = label,
      n = length(wer_values),
      Mean = round(mean(wer_values), 2),
      SD = round(sd(wer_values), 2),
      Median = round(median(wer_values), 2),
      Q1 = round(quantile(wer_values, 0.25), 2),
      Q3 = round(quantile(wer_values, 0.75), 2),
      IQR = round(quantile(wer_values, 0.75) - quantile(wer_values, 0.25), 2)
    )
    wer_stats <- rbind(wer_stats, stats_row)
  }
}

# Print and save table
print(wer_stats)
write.table(wer_stats, file.path(output_dir, "table2_wer_statistics.tsv"), 
            sep = "\t", row.names = FALSE, quote = FALSE)

# Create boxplot
boxplot_data <- list()
sample_sizes <- c()
for(label in event_types) {
  boxplot_data[[label]] <- data[data$Label == label, "WER_Percent"]
  sample_sizes <- c(sample_sizes, sum(data$Label == label))
}

label_with_n <- paste0(event_types, "\n(n=", sample_sizes, ")")

# Save boxplot
png(file.path(output_dir, "figure2_wer_boxplot.png"), width = 800, height = 600)

#set boxplot
boxplot(boxplot_data,
        names = label_with_n,
        main = "WER Distribution by Event Type",
        xlab = "Event Type",
        ylab = "WER (%)",
        col = "steelblue",
        ylim = c(0, 120))
dev.off()

cat("Files saved:\n- table2_wer_statistics.tsv\n- figure2_wer_boxplot.png\n")


# Step 3: error type composition
# Set paths
data_dir <- "../descriptive_stat"
output_dir <- "../descriptive_stat/output"

# Read data
data <- read.delim(file.path(data_dir, "master_analysis_table.tsv"), sep = "\t")
data$Deletions <- as.numeric(data$Deletions)
data$Insertions <- as.numeric(data$Insertions)
data$Substitutions <- as.numeric(data$Substitutions)
data <- data[!is.na(data$Deletions) & !is.na(data$Insertions) & !is.na(data$Substitutions), ]

# Calculate error composition
event_types <- c("laughter", "sigh", "throat-clearing")
error_composition <- data.frame()

for(label in event_types) {
  label_data <- data[data$Label == label, ]
  
  if(nrow(label_data) > 0) {
    total_deletions <- sum(label_data$Deletions, na.rm = TRUE)
    total_insertions <- sum(label_data$Insertions, na.rm = TRUE)
    total_substitutions <- sum(label_data$Substitutions, na.rm = TRUE)
    total_errors <- total_deletions + total_insertions + total_substitutions
    
    if(total_errors > 0) {
      d_percent <- round((total_deletions / total_errors) * 100, 1)
      i_percent <- round((total_insertions / total_errors) * 100, 1)
      s_percent <- round((total_substitutions / total_errors) * 100, 1)
    } else {
      d_percent <- i_percent <- s_percent <- 0
    }
    
    composition_row <- data.frame(
      Label = label,
      Deletions = total_deletions,
      Insertions = total_insertions,
      Substitutions = total_substitutions,
      TotalErrors = total_errors,
      D_Percent = d_percent,
      I_Percent = i_percent,
      S_Percent = s_percent
    )
    
    error_composition <- rbind(error_composition, composition_row)
  }
}

# Print and save table
print(error_composition)
write.table(error_composition, file.path(output_dir, "table3_error_composition.tsv"), 
            sep = "\t", row.names = FALSE, quote = FALSE)

# Create stacked bar chart
chart_data <- as.matrix(error_composition[, c("D_Percent", "I_Percent", "S_Percent")])
rownames(chart_data) <- error_composition$Label

# Save plot
png(file.path(output_dir, "figure3_error_composition.png"), width = 800, height = 600)
barplot(t(chart_data),
        main = "Error Type Composition by Event Type",
        xlab = "Event Type", ylab = "Percentage (%)",
        col = c("coral", "steelblue", "lightgreen"),
        legend.text = c("Deletions", "Insertions", "Substitutions"),
        args.legend = list(x = "topright", bty = "n"),
        ylim = c(0, 100), beside = FALSE)
dev.off()

cat("Files saved:\n- table3_error_composition.tsv\n- figure3_error_composition.png\n")

# Step 4: Temporal position of errors

# Set path
data_dir <- "../descriptive_stat"
output_dir <- "../descriptive_stat/output"

# Read and prepare data
data <- read.delim(file.path(data_dir, "master_analysis_table.tsv"), sep = "\t")
data$Feature_Before <- as.numeric(data$Feature_Before)
data$Feature_During <- as.numeric(data$Feature_During)
data$Feature_After <- as.numeric(data$Feature_After)
data$No_Feature <- as.numeric(data$No_Feature)

# Calculate by event type
event_types <- c("laughter", "sigh", "throat-clearing")
results <- data.frame()

for(label in event_types) {
  label_data <- data[data$Label == label, ]
  
  before <- sum(label_data$Feature_Before, na.rm = TRUE)
  during <- sum(label_data$Feature_During, na.rm = TRUE)
  after <- sum(label_data$Feature_After, na.rm = TRUE)
  no_feature <- sum(label_data$No_Feature, na.rm = TRUE)
  total <- before + during + after + no_feature
  
  row <- data.frame(
    Label = label,
    Before = before, During = during, After = after, No_Feature = no_feature,
    Row_Total = total,
    Before_Pct = round((before/total)*100, 1),
    During_Pct = round((during/total)*100, 1),
    After_Pct = round((after/total)*100, 1),
    No_Feature_Pct = round((no_feature/total)*100, 1)
  )
  results <- rbind(results, row)
}

# Save table
print(results)
write.table(results, file.path(output_dir, "table4_temporal_position.tsv"), 
            sep = "\t", row.names = FALSE, quote = FALSE)

# Create plot
chart_data <- as.matrix(results[, c("Before_Pct", "During_Pct", "After_Pct", "No_Feature_Pct")])
rownames(chart_data) <- results$Label

png(file.path(output_dir, "figure4_temporal_position.png"), width = 800, height = 600)
barplot(t(chart_data),
        main = "Temporal Position of Errors by Event Type",
        xlab = "Event Type", ylab = "Percentage (%)",
        col = c("lightgreen", "orange", "coral", "lightgray"),
        legend.text = c("Before", "During", "After", "No Feature"),
        args.legend = list(x = "topright", bty = "n"),
        ylim = c(0, 100), beside = FALSE)
dev.off()

cat("Files saved: table4_temporal_position.tsv, figure4_temporal_position.png\n")

# Step 5: Acoustic features vs. error occurrence

# Set path
data_dir <- "../descriptive_stat"
output_dir <- "../descriptive_stat/output"

# Read data
data <- read.delim(file.path(data_dir, "master_analysis_table.tsv"), sep = "\t")

#Prepare data
data$AnyError <- ifelse(data$Total_WER_Errors > 0, 1, 0)

# Filter undefined labels
data <- data[data$Mean_HNR != "--undefined--" & 
             data$Spectral_Tilt != "--undefined--" & 
             data$ZCR != "--undefined--", ]

# Convert to numeric
data$Mean_HNR <- as.numeric(data$Mean_HNR)
data$Spectral_Tilt <- as.numeric(data$Spectral_Tilt)
data$ZCR <- as.numeric(data$ZCR)
data <- data[!is.na(data$Mean_HNR) & !is.na(data$Spectral_Tilt) & !is.na(data$ZCR), ]

# Create plots
features <- list(
  list(data = data$Mean_HNR, name = "HNR", ylabel = "HNR (dB)", file = "figure5_hnr_anyerror.png"),
  list(data = data$Spectral_Tilt, name = "Spectral Tilt", ylabel = "Spectral Tilt (dB)", file = "figure6_tilt_anyerror.png"),
  list(data = data$ZCR, name = "ZCR", ylabel = "ZCR (Hz)", file = "figure7_zcr_anyerror.png")
)

for(feature in features) {
  no_error <- feature$data[data$AnyError == 0]
  with_error <- feature$data[data$AnyError == 1]
  plot_data <- list("0" = no_error, "1" = with_error)
  
  png(file.path(output_dir, feature$file), width = 800, height = 600)
  boxplot(plot_data,
          main = paste("Distribution of", feature$name, "by Error Occurrence"),
          xlab = "Any Error (0 = no error; 1 = at least one error)",
          ylab = feature$ylabel,
          col = c("steelblue", "coral"))
  dev.off()
}

cat("Files saved:\n- figure5_hnr_anyerror.png\n- figure6_tilt_anyerror.png\n- figure7_zcr_anyerror.png\n")

# Step 6: spearman correlation

# Set path
data_dir <- "../descriptive_stat"
output_dir <- "../descriptive_stat/output"

# Read data
data <- read.delim(file.path(data_dir, "master_analysis_table.tsv"), sep = "\t")

# Prepare data
data <- data[data$Mean_HNR != "--undefined--" & 
             data$Spectral_Tilt != "--undefined--" & 
             data$ZCR != "--undefined--", ]

data$Mean_HNR <- as.numeric(data$Mean_HNR)
data$Spectral_Tilt <- as.numeric(data$Spectral_Tilt)
data$ZCR <- as.numeric(data$ZCR)
data <- data[!is.na(data$Mean_HNR) & !is.na(data$Spectral_Tilt) & !is.na(data$ZCR), ]

# Calculate correlation
cor1 <- cor.test(data$Mean_HNR, data$Spectral_Tilt, method = "spearman")
cor2 <- cor.test(data$Mean_HNR, data$ZCR, method = "spearman")
cor3 <- cor.test(data$Spectral_Tilt, data$ZCR, method = "spearman")

# Create table
results <- data.frame(
  Pair = c("HNR — Spectral Tilt", "HNR — ZCR", "Spectral Tilt — ZCR"),
  rho = round(c(cor1$estimate, cor2$estimate, cor3$estimate), 3),
  p = round(c(cor1$p.value, cor2$p.value, cor3$p.value), 3)
)

# Print and save table
print(results)
write.table(results, file.path(output_dir, "table5_spearman_correlations.tsv"), 
            sep = "\t", row.names = FALSE, quote = FALSE)

cat("File saved: table5_spearman_correlations.tsv\n")
