Kendall's Tau (τ)

# Kendall's Tau-b: Ordinal correlation with tied ranks
# Sleep quality and academic performance in college students

library(tidyverse)
library(DescTools)   # For KendallTauB with CI
library(psych)       # For corr.test

# Simulate realistic data (or load: data <- read.csv("sleep_gpa.csv"))
set.seed(2025)
n <- 150

# Sleep quality (1-5 Likert) - clustered at 3 (average)
sleep_quality <- sample(1:5, n, replace=TRUE, 
                        prob=c(0.10, 0.20, 0.40, 0.20, 0.10))

# GPA categories (1-6) correlated with sleep (monotonic relationship)
true_latent_gpa <- 0.5 * sleep_quality + rnorm(n, 0, 0.8)
gpa_category <- cut(true_latent_gpa, 
                    breaks=c(-Inf, 1.5, 2.0, 2.5, 3.0, 3.5, Inf),
                    labels=c("<2.0", "2.0-2.5", "2.5-3.0", "3.0-3.5", "3.5-4.0", "4.0"))
gpa_numeric <- as.numeric(gpa_category)

data <- data.frame(
  student_id = 1:n,
  sleep_quality = factor(sleep_quality, levels=1:5, 
                         labels=c("Very poor", "Poor", "Fair", "Good", "Excellent")),
  sleep_numeric = sleep_quality,
  gpa_category = gpa_category,
  gpa_numeric = gpa_numeric
)

head(data, 10)

# === STEP 1: Check Assumptions ===

# 1. Check for tied ranks
cat("=== Frequency of Sleep Quality Ratings ===\n")
table(data$sleep_numeric)
cat("\nProportion of tied values in sleep:", 
    1 - length(unique(data$sleep_numeric))/n, "\n")

cat("\n=== Frequency of GPA Categories ===\n")
table(data$gpa_numeric)

# 2. Visual check for monotonicity
ggplot(data, aes(x=sleep_numeric, y=gpa_numeric)) +
  geom_jitter(width=0.2, height=0.2, alpha=0.4, size=2) +
  geom_smooth(method="loess", color="red", se=TRUE) +
  labs(title="Sleep Quality vs Academic Performance",
       subtitle="Red curve shows monotonic trend",
       x="Sleep Quality(1=Very Poor to 5=Excellent)",
       y="GPA Category(1=<2.0 to 6=4.0)") +
  theme_classic()

# Contingency table heatmap
tab <- table(data$sleep_numeric, data$gpa_numeric)
pheatmap::pheatmap(tab, 
                   main="Frequency Heatmap: Sleep × GPA",
                   display_numbers=TRUE, 
                   cluster_rows=FALSE, 
                   cluster_cols=FALSE)

# === STEP 2: Compute Kendall's Tau-b ===

# Method 1: Base R cor.test
result <- cor.test(data$sleep_numeric, data$gpa_numeric, 
                   method="kendall", alternative="two.sided")
print(result)

cat("\n=== Kendall's Tau-b Results ===\n")
cat(sprintf("τb = %.3f\n", result$estimate))
cat(sprintf("z = %.2f, p = %.4f\n", result$statistic, result$p.value))

# Method 2: DescTools (provides CI)
tau_ci <- KendallTauB(data$sleep_numeric, data$gpa_numeric, conf.level=0.95)
cat(sprintf("τb = %.3f, 95%% CI [%.3f, %.3f]\n", 
            tau_ci[1], tau_ci[2], tau_ci[3]))

# === STEP 3: Interpretation ===

cat("\n=== Interpretation Guidelines ===\n")
cat("Tau magnitude: 0.1=small, 0.3=medium, 0.5=large(Cohen, 1988)\n")
cat("Note: Tau < Pearson r for same data(tau is more conservative)\n\n")

tau_val <- result$estimate
if (abs(tau_val) < 0.1) {
  strength <- "negligible"
} else if (abs(tau_val) < 0.3) {
  strength <- "small"
} else if (abs(tau_val) < 0.5) {
  strength <- "moderate"
} else {
  strength <- "large"
}

cat(sprintf("Effect size: %s(%s association)\n", 
            strength, ifelse(tau_val > 0, "positive", "negative")))

# === STEP 4: Concordant vs Discordant Pairs ===

# Calculate manually
n_pairs <- n * (n - 1) / 2
concordant <- sum(outer(data$sleep_numeric, data$sleep_numeric, "<") & 
                  outer(data$gpa_numeric, data$gpa_numeric, "<"))
discordant <- sum(outer(data$sleep_numeric, data$sleep_numeric, "<") & 
                  outer(data$gpa_numeric, data$gpa_numeric, ">"))

cat(sprintf("\nTotal pairs: %.0f\n", n_pairs))
cat(sprintf("Concordant pairs: %d(%.1f%%)\n", concordant, 100*concordant/n_pairs))
cat(sprintf("Discordant pairs: %d(%.1f%%)\n", discordant, 100*discordant/n_pairs))
cat(sprintf("Tied pairs: %d(%.1f%%)\n", 
            n_pairs - concordant - discordant, 
            100*(n_pairs - concordant - discordant)/n_pairs))

# === STEP 5: Compare with Spearman's rho ===

spearman_result <- cor.test(data$sleep_numeric, data$gpa_numeric, method="spearman")
cat(sprintf("\nSpearman's ρ = %.3f (p = %.4f)\n", 
            spearman_result$estimate, spearman_result$p.value))
cat("Note: Spearman's rho typically 1.5× larger than Kendall's tau for same data\n")
cat(sprintf("Actual ratio: %.2f\n", spearman_result$estimate / tau_val))

# === STEP 6: Sensitivity Analysis (remove outliers) ===

# Identify outliers using boxplot rule
outliers_sleep <- boxplot.stats(data$sleep_numeric)$out
outliers_gpa <- boxplot.stats(data$gpa_numeric)$out

if (length(outliers_sleep) > 0 | length(outliers_gpa) > 0) {
  cat("\n=== Sensitivity Analysis(excluding outliers) ===\n")
  data_no_outliers <- data[!(data$sleep_numeric %in% outliers_sleep | 
                             data$gpa_numeric %in% outliers_gpa), ]
  tau_no_outliers <- cor.test(data_no_outliers$sleep_numeric, 
                              data_no_outliers$gpa_numeric, 
                              method="kendall")$estimate
  cat(sprintf("τb without outliers = %.3f (diff = %.3f)\n", 
              tau_no_outliers, tau_no_outliers - tau_val))
} else {
  cat("\nNo outliers detected.\n")
}

# === APA-Style Reporting ===

cat("\n=== APA-Style Report ===\n")
cat(sprintf(
"A Kendall's tau-b correlation was computed to assess the association between 
sleep quality and academic performance in 150 college students. Both variables 
were measured on ordinal scales(sleep quality: 1-5 Likert; GPA: 6 ordered 
categories). There was a significant positive monotonic association, τb = %.2f, 
95%% CI [%.2f, %.2f], z = %.2f, p %s, indicating that students with better 
sleep quality tended to have higher academic performance. The effect size was %s 
(Cohen, 1988). The analysis revealed that %.0f%% of pairs were concordant 
(both variables increased together), %.0f%% were discordant, and %.0f%% were 
tied. These findings support the link between sleep quality and academic success, 
consistent with prior research(Hershner & Chervin, 2014).\n",
tau_val, tau_ci[2], tau_ci[3], result$statistic, 
ifelse(result$p.value < 0.001, "< .001", sprintf("= %.3f", result$p.value)),
strength,
100*concordant/n_pairs,
100*discordant/n_pairs,
100*(n_pairs - concordant - discordant)/n_pairs
))