Spearman Correlation (ρ)

# Spearman Correlation: Ordinal Likert Scale Data
# Job Satisfaction and Work-Life Balance

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

# Simulate realistic survey data
set.seed(2025)
n <- 120

# Work-life balance: 7-point Likert (skewed toward middle-upper)
wlb_score <- sample(1:7, n, replace=TRUE, 
                    prob=c(0.05, 0.10, 0.15, 0.25, 0.25, 0.15, 0.05))

# Job satisfaction: monotonically associated with WLB (with noise)
latent_satisfaction <- 0.6 * wlb_score + rnorm(n, 0, 1.2)
job_sat <- round(pmax(1, pmin(latent_satisfaction + 2, 7)))

data <- data.frame(
  employee_id = 1:n,
  wlb_score = wlb_score,
  job_sat = job_sat
)

head(data, 10)

cat("\n=== Descriptive Statistics ===\n")
cat("Work-Life Balance Distribution:\n")
print(table(data$wlb_score))
cat("\nJob Satisfaction Distribution:\n")
print(table(data$job_sat))

psych::describe(data[, c("wlb_score", "job_sat")])

# === STEP 1: Visual Exploration ===

cat("\n=== Assumption Checks ===\n")

# Scatterplot with jitter (ordinal data)
ggplot(data, aes(x=wlb_score, y=job_sat)) +
  geom_jitter(width=0.2, height=0.2, alpha=0.4, size=2.5) +
  geom_smooth(method="loess", color="red", se=TRUE) +
  scale_x_continuous(breaks=1:7) +
  scale_y_continuous(breaks=1:7) +
  labs(title="Work-Life Balance vs Job Satisfaction",
       subtitle="Red curve shows monotonic trend(Spearman appropriate)",
       x="Work-Life Balance(1=Very Low to 7=Very High)",
       y="Job Satisfaction(1=Very Low to 7=Very High)") +
  theme_classic(base_size=12)

# Contingency table heatmap
tab <- table(data$wlb_score, data$job_sat)
pheatmap::pheatmap(tab, 
                   main="Frequency Heatmap: WLB × Job Satisfaction",
                   display_numbers=TRUE,
                   cluster_rows=FALSE,
                   cluster_cols=FALSE,
                   color=colorRampPalette(c("white", "orange", "red"))(50))

# Check for ties
cat(sprintf("\nProportion of tied values:\n"))
cat(sprintf("  WLB: %.1f%% (%d unique values out of %d)\n",
            (1 - length(unique(data$wlb_score))/n)*100,
            length(unique(data$wlb_score)), n))
cat(sprintf("  Job Sat: %.1f%% (%d unique values out of %d)\n",
            (1 - length(unique(data$job_sat))/n)*100,
            length(unique(data$job_sat)), n))
cat("Many ties expected with 7-point Likert - Spearman handles via average ranks.\n")

# === STEP 2: Spearman Correlation ===

cat("\n=== Spearman Rank Correlation ===\n")

# Method 1: Base R cor.test
result <- cor.test(data$wlb_score, data$job_sat, 
                   method="spearman", 
                   alternative="two.sided",
                   exact=FALSE)  # Use asymptotic for n>30
print(result)

rho <- result$estimate
cat(sprintf("\nρs = %.3f\n", rho))
cat(sprintf("S = %.0f, p %s\n", 
            result$statistic,
            ifelse(result$p.value < 0.001, "< .001", 
                   sprintf("= %.4f", result$p.value))))

# Method 2: DescTools for CI
rho_ci <- SpearmanRho(data$wlb_score, data$job_sat, conf.level=0.95)
cat(sprintf("\nρs = %.3f, 95%% CI [%.3f, %.3f]\n",
            rho_ci[1], rho_ci[2], rho_ci[3]))

# === STEP 3: Effect Size Interpretation ===

cat("\n=== Effect Size Guidelines ===\n")
cat("Spearman ρ: 0.10=small, 0.30=medium, 0.50=large(Cohen, 1988)\n\n")

if (abs(rho) < 0.10) {
  strength <- "negligible"
} else if (abs(rho) < 0.30) {
  strength <- "small"
} else if (abs(rho) < 0.50) {
  strength <- "medium"
} else {
  strength <- "large"
}

direction <- ifelse(rho > 0, "positive", "negative")
cat(sprintf("Observed: %s(%s monotonic association)\n", strength, direction))

# === STEP 4: Compare with Pearson (Sensitivity Check) ===

cat("\n=== Comparison: Spearman vs Pearson ===\n")
pearson_result <- cor.test(data$wlb_score, data$job_sat, method="pearson")
r_pearson <- pearson_result$estimate

cat(sprintf("Spearman ρ = %.3f\n", rho))
cat(sprintf("Pearson r  = %.3f\n", r_pearson))
cat(sprintf("Difference = %.3f\n", abs(rho - r_pearson)))
cat("\nInterpretation:\n")
if (abs(rho - r_pearson) < 0.05) {
  cat("- Small difference: Linear and monotonic associations similar.\n")
  cat("- Either method acceptable; Spearman preferred for ordinal data.\n")
} else {
  cat("- Moderate difference suggests non-linearity or non-normality.\n")
  cat("- Spearman preferred as more robust for ordinal Likert scales.\n")
}

# === STEP 5: Visualize Ranks (What Spearman Analyzes) ===

data$wlb_rank <- rank(data$wlb_score)
data$jobsat_rank <- rank(data$job_sat)

ggplot(data, aes(x=wlb_rank, y=jobsat_rank)) +
  geom_point(alpha=0.5, size=2) +
  geom_smooth(method="lm", color="blue", se=TRUE) +
  labs(title="Rank Space: What Spearman Correlation Analyzes",
       subtitle=sprintf("ρ = %.2f (Pearson correlation on ranks)", rho),
       x="Work-Life Balance Rank",
       y="Job Satisfaction Rank") +
  theme_classic()

# === STEP 6: Bootstrap CI (Robustness Check) ===

cat("\n=== Bootstrap 95% CI(n=1000 resamples) ===\n")
set.seed(2025)
boot_rho <- replicate(1000, {
  indices <- sample(1:n, n, replace=TRUE)
  cor(data$wlb_score[indices], data$job_sat[indices], method="spearman")
})

boot_ci <- quantile(boot_rho, c(0.025, 0.975))
cat(sprintf("Bootstrap CI: [%.3f, %.3f]\n", boot_ci[1], boot_ci[2]))
cat(sprintf("Parametric CI: [%.3f, %.3f]\n", rho_ci[2], rho_ci[3]))
cat("Bootstrap and parametric CIs should be similar for adequate sample size.\n")

# === STEP 7: Compare with Kendall's Tau (Alternative Rank Method) ===

kendall_result <- cor.test(data$wlb_score, data$job_sat, method="kendall")
tau <- kendall_result$estimate

cat(sprintf("\nKendall's τb = %.3f (alternative rank correlation)\n", tau))
cat(sprintf("Ratio ρ/τ = %.2f (typically ≈1.5 for same data)\n", rho/tau))
cat("Spearman and Kendall should have same sign and consistent magnitude.\n")

# === STEP 8: Detailed Visualizations ===

# Grouped bar plot by WLB category
data_summary <- data %>%
  group_by(wlb_score) %>%
  summarize(mean_jobsat = mean(job_sat),
            se_jobsat = sd(job_sat)/sqrt(n()),
            .groups='drop')

ggplot(data_summary, aes(x=factor(wlb_score), y=mean_jobsat)) +
  geom_bar(stat="identity", fill="steelblue", alpha=0.7) +
  geom_errorbar(aes(ymin=mean_jobsat - 1.96*se_jobsat,
                    ymax=mean_jobsat + 1.96*se_jobsat),
                width=0.3) +
  labs(title="Monotonic Trend: Work-Life Balance → Job Satisfaction",
       subtitle=sprintf("Spearman ρ = %.2f***", rho),
       x="Work-Life Balance Score",
       y="Mean Job Satisfaction(±95% CI)") +
  theme_classic()

# === APA-Style Reporting ===

cat("\n=== APA-Style Report ===\n")
cat(sprintf(
"A Spearman rank-order correlation was computed to assess the relationship 
between work-life balance and job satisfaction in %d employees. Both variables 
were measured on 7-point Likert scales(ordinal data), making Spearman 
correlation appropriate over Pearson correlation which assumes interval-level 
measurement. Scatterplot inspection confirmed a monotonic positive relationship. 
There was a significant positive correlation between work-life balance and job 
satisfaction, ρs(%d) = %.2f, 95%% CI [%.2f, %.2f], p < .001. The effect size 
was %s(Cohen, 1988), indicating that employees with better work-life balance 
tended to report higher job satisfaction. This monotonic association was robust, 
as confirmed by Kendall's tau-b(τb = %.2f). Spearman's rank correlation was 
preferred over Pearson(r = %.2f) due to the ordinal nature of Likert scales. 
These findings support organizational interventions targeting work-life balance 
to improve employee job satisfaction, consistent with meta-analytic research 
showing moderate-to-strong associations(ρ ≈ 0.45-0.55).\n",
n, n-2, rho, rho_ci[2], rho_ci[3], strength, tau, r_pearson
))