Somers' D (Dyx)

# Somers' D: Education Level (IV) → Income Category (DV)
# Asymmetric ordinal association with clear predictor-outcome direction

library(tidyverse)
library(DescTools)   # For SomersDelta
library(vcd)         # For mosaic plots
library(psych)       # For polychoric correlation

set.seed(2025)
n <- 200

# Simulate education levels (IV: 1-5)
# Distribution: roughly normal, centered at 3 (Some College)
education <- sample(1:5, n, replace = TRUE, 
                    prob = c(0.15, 0.25, 0.30, 0.20, 0.10))

# Simulate income categories (DV: 1-4) strongly predicted by education
# Higher education → higher income (monotonic relationship)
latent_income <- 0.8 * education + rnorm(n, 0, 0.9)
income <- cut(latent_income, 
              breaks = c(-Inf, 1.5, 2.5, 3.5, Inf),
              labels = 1:4)
income_numeric <- as.numeric(income)

data <- data.frame(
  id = 1:n,
  education = factor(education, levels = 1:5,
                     labels = c("<HS", "HS", "Some College", "Bachelor's", "Graduate")),
  education_num = education,
  income = factor(income_numeric, levels = 1:4,
                  labels = c("<$30k", "$30-50k", "$50-75k", ">$75k")),
  income_num = income_numeric
)

head(data, 10)

# === STEP 1: Descriptive Statistics ===

cat("=== FREQUENCY DISTRIBUTIONS ===\n")
cat("\nEducation(IV):\n")
table(data$education)
cat("\nIncome(DV):\n")
table(data$income)

# === STEP 2: Cross-Tabulation (DV as rows, IV as columns) ===

cat("\n=== CONTINGENCY TABLE(Rows=DV, Columns=IV) ===\n")
contab <- table(data$income, data$education)
print(contab)

cat("\n=== ROW PERCENTAGES(% within each income level) ===\n")
print(round(prop.table(contab, margin = 1) * 100, 1))

cat("\n=== COLUMN PERCENTAGES(% within each education level) ===\n")
print(round(prop.table(contab, margin = 2) * 100, 1))

# Visualize with heatmap
library(pheatmap)
pheatmap(contab, 
         cluster_rows = FALSE, 
         cluster_cols = FALSE,
         display_numbers = TRUE,
         main = "Income × Education Contingency Table(Frequencies)",
         xlab = "Education(IV)",
         ylab = "Income(DV)")

# Mosaic plot (area proportional to frequency)
mosaic(~ income + education, data = data,
       shade = TRUE, legend = TRUE,
       main = "Mosaic Plot: Income(DV) × Education(IV)")

# === STEP 3: Compute Somers' D (Dyx: DV=income, IV=education) ===

cat("\n=== SOMERS' D COMPUTATION ===\n")

# Somers' Dyx: DV=income (Y), IV=education (X)
# Adjusts for ties on DV (income) only
somers_dyx <- SomersDelta(data$education_num, data$income_num, 
                          conf.level = 0.95)

cat("Somers' Dyx(DV=income, IV=education):\n")
cat(sprintf("  Dyx = %.3f\n", somers_dyx[1]))
cat(sprintf("  95%% CI: [%.3f, %.3f]\n", somers_dyx[2], somers_dyx[3]))

# Reverse: Somers' Dxy (DV=education, IV=income)
somers_dxy <- SomersDelta(data$income_num, data$education_num,
                          conf.level = 0.95)

cat("\nSomers' Dxy(DV=education, IV=income - reversed):\n")
cat(sprintf("  Dxy = %.3f\n", somers_dxy[1]))
cat(sprintf("  95%% CI: [%.3f, %.3f]\n", somers_dxy[2], somers_dxy[3]))

cat("\nNote: Dyx ≠ Dxy because tie adjustments are asymmetric.\n")
cat("Use Dyx when income is the dependent variable(outcome).\n")

# === STEP 4: Compare with Related Measures ===

cat("\n=== COMPARISON WITH SYMMETRIC MEASURES ===\n")

# Kendall's tau-b (symmetric, adjusts for ties on both variables)
tau_b <- cor(data$education_num, data$income_num, method = "kendall")
cat(sprintf("Kendall's tau-b(symmetric): %.3f\n", tau_b))

# Goodman-Kruskal gamma (symmetric, ignores all ties)
library(DescTools)
gamma_val <- GoodmanKruskalGamma(data$education_num, data$income_num,
                                 conf.level = 0.95)
cat(sprintf("Goodman-Kruskal gamma: %.3f [%.3f, %.3f]\n",
            gamma_val[1], gamma_val[2], gamma_val[3]))

cat("\nRelationships:\n")
cat("  - Gamma ignores all ties → typically largest magnitude\n")
cat("  - Tau-b adjusts for ties on both variables → smallest magnitude\n")
cat("  - Somers' Dyx adjusts for DV ties only → intermediate\n")
cat(sprintf("  - Average of Dyx and Dxy ≈ tau-b: (%.3f + %.3f) / 2 = %.3f ≈ %.3f\n",
            somers_dyx[1], somers_dxy[1], 
            (somers_dyx[1] + somers_dxy[1]) / 2, tau_b))

# === STEP 5: Statistical Significance Testing ===

cat("\n=== HYPOTHESIS TEST ===\n")

# Use asymptotic z-test for Somers' D
# Standard error approximation (for large samples)
n_pairs <- n * (n - 1) / 2

# Calculate concordant and discordant pairs manually
concordant <- sum(outer(data$education_num, data$education_num, "<") & 
                  outer(data$income_num, data$income_num, "<")) +
              sum(outer(data$education_num, data$education_num, ">") & 
                  outer(data$income_num, data$income_num, ">"))

discordant <- sum(outer(data$education_num, data$education_num, "<") & 
                  outer(data$income_num, data$income_num, ">")) +
              sum(outer(data$education_num, data$education_num, ">") & 
                  outer(data$income_num, data$income_num, "<"))

cat(sprintf("Total 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))

# Asymptotic test (large sample approximation)
# Note: Exact SE calculation complex; using bootstrap or DescTools CI
ci_width <- somers_dyx[3] - somers_dyx[2]
se_approx <- ci_width / (2 * 1.96)
z_stat <- somers_dyx[1] / se_approx
p_value <- 2 * pnorm(-abs(z_stat))

cat(sprintf("\nAsymptotic z-test:\n"))
cat(sprintf("  z = %.2f\n", z_stat))
cat(sprintf("  p-value = %.4f\n", p_value))

if (p_value < 0.001) {
  cat("  Result: Highly significant(p < .001)\n")
} else if (p_value < 0.05) {
  cat("  Result: Significant(p < .05)\n")
} else {
  cat("  Result: Not significant(p ≥ .05)\n")
}

# === STEP 6: Effect Size Interpretation ===

cat("\n=== EFFECT SIZE INTERPRETATION ===\n")
cat("Somers' D magnitude guidelines(similar to tau):\n")
cat("  |D| < 0.1: negligible\n")
cat("  0.1 ≤ |D| < 0.3: small\n")
cat("  0.3 ≤ |D| < 0.5: moderate\n")
cat("  |D| ≥ 0.5: large\n\n")

D_val <- abs(somers_dyx[1])
if (D_val < 0.1) {
  strength <- "negligible"
} else if (D_val < 0.3) {
  strength <- "small"
} else if (D_val < 0.5) {
  strength <- "moderate"
} else {
  strength <- "large"
}

cat(sprintf("Observed Dyx = %.3f: %s effect size\n", somers_dyx[1], strength))

# PRE (Proportional Reduction in Error) interpretation
cat("\nPRE Interpretation:\n")
cat(sprintf("Knowing education(IV) reduces error in predicting income(DV) by %.1f%%\n",
            abs(somers_dyx[1]) * 100))

# === STEP 7: Ordinal Regression Comparison ===

cat("\n=== ORDINAL LOGISTIC REGRESSION(Proportional Odds Model) ===\n")

library(MASS)
model <- polr(income ~ education, data = data, Hess = TRUE)
summary(model)

cat("\nNote: Ordinal regression provides coefficient estimates and tests\n")
cat("whether education significantly predicts income, controlling for\n")
cat("proportional odds assumption. Consistent with Somers' D results.\n")

# === STEP 8: Visualization of Association ===

cat("\n=== VISUALIZATIONS ===\n")

# Stacked bar chart showing income distribution by education
ggplot(data, aes(x = education, fill = income)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent) +
  labs(title = "Income Distribution by Education Level",
       subtitle = "Clear monotonic trend: Higher education → Higher income",
       x = "Education Level(IV)",
       y = "Proportion",
       fill = "Income(DV)") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Cumulative proportion plot to visualize monotonicity
cumulative_props <- data %>%
  group_by(education_num, income_num) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(education_num) %>%
  mutate(prop = n / sum(n),
         cum_prop = cumsum(prop))

ggplot(cumulative_props, aes(x = education_num, y = cum_prop, 
                             color = factor(income_num))) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 3) +
  labs(title = "Cumulative Income Proportions by Education",
       subtitle = "Monotonic pattern confirms positive association",
       x = "Education Level(1=<HS to 5=Graduate)",
       y = "Cumulative Proportion",
       color = "Income") +
  theme_classic()

# === APA-STYLE REPORTING ===

cat("\n=== APA-STYLE REPORT ===\n")
cat(sprintf(
"Somers' D was computed to assess the asymmetric ordinal association between 
education level(independent variable: 5 ordered categories from <HS to Graduate 
degree) and income category(dependent variable: 4 ordered categories from <$30k 
to >$75k) in a sample of 200 working adults. Somers' D designates income as the 
dependent variable and adjusts for ties on the dependent variable only, making it 
appropriate when a clear predictor-outcome relationship is hypothesized.

A strong positive monotonic association was found, Dyx = %.2f, 95%% CI [%.2f, %.2f], 
z = %.2f, p < .001, indicating that higher education levels were strongly associated 
with higher income categories. The effect size was %s according to standard 
interpretation guidelines. Using the proportional reduction in error(PRE) 
interpretation, knowledge of education level reduced prediction error for income 
category by approximately %.0f%%.

The analysis revealed that %.0f%% of observation pairs were concordant(both 
education and income ranked in same direction), %.0f%% were discordant, and %.0f%% 
involved ties. Somers' Dyx(%.2f) was slightly larger than Kendall's tau-b(%.2f) 
but smaller than Goodman-Kruskal gamma(%.2f), reflecting the asymmetric 
adjustment for ties on the dependent variable only. Results were consistent with 
ordinal logistic regression analysis, supporting the conclusion that educational 
attainment is a strong predictor of income category, consistent with U.S. Census 
data patterns(U.S. Census Bureau, 2020).\n",
somers_dyx[1], somers_dyx[2], somers_dyx[3], z_stat, strength,
abs(somers_dyx[1]) * 100,
100 * concordant / n_pairs,
100 * discordant / n_pairs,
100 * (n_pairs - concordant - discordant) / n_pairs,
somers_dyx[1], tau_b, gamma_val[1]
))