Goodman-Kruskal Gamma (γ)

# Goodman-Kruskal Gamma: Political Ideology × Policy Support
# Symmetric ordinal association (no clear DV/IV)

library(tidyverse)
library(DescTools)   # For GoodmanKruskalGamma
library(vcd)         # For mosaic plots
library(MASS)        # For ordinal regression

set.seed(2025)
n <- 180

# Political ideology (1-5: very liberal to very conservative)
ideology <- sample(1:5, n, replace = TRUE,
                   prob = c(0.20, 0.25, 0.30, 0.15, 0.10))

# Policy support (1-5) strongly associated with ideology
# More liberal (low ideology) → more support (high policy)
latent_support <- 6 - 0.9 * ideology + rnorm(n, 0, 0.8)
policy_support <- cut(latent_support,
                      breaks = c(-Inf, 1.5, 2.5, 3.5, 4.5, Inf),
                      labels = 1:5)
policy_numeric <- as.numeric(policy_support)

data <- data.frame(
  respondent_id = 1:n,
  ideology = factor(ideology, levels = 1:5,
                    labels = c("Very Liberal", "Liberal", "Moderate",
                              "Conservative", "Very Conservative")),
  ideology_num = ideology,
  policy_support = factor(policy_numeric, levels = 1:5,
                         labels = c("Strongly Oppose", "Oppose", "Neutral",
                                   "Support", "Strongly Support")),
  policy_num = policy_numeric
)

head(data, 10)

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

cat("=== FREQUENCY DISTRIBUTIONS ===\n")
cat("\nPolitical Ideology:\n")
table(data$ideology)
cat("\nEnvironmental Policy Support:\n")
table(data$policy_support)

# === STEP 2: Cross-Tabulation ===

cat("\n=== CONTINGENCY TABLE ===\n")
contab <- table(data$ideology, data$policy_support)
print(contab)

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

# Heatmap
library(pheatmap)
pheatmap(contab,
         cluster_rows = FALSE,
         cluster_cols = FALSE,
         display_numbers = TRUE,
         main = "Ideology × Policy Support(Frequencies)",
         xlab = "Policy Support",
         ylab = "Political Ideology")

# Mosaic plot
mosaic(~ ideology + policy_support, data = data,
       shade = TRUE, legend = TRUE,
       main = "Mosaic: Political Ideology × Policy Support")

# === STEP 3: Compute Goodman-Kruskal Gamma ===

cat("\n=== GOODMAN-KRUSKAL GAMMA ===\n")

# Gamma (ignores all tied pairs)
gamma_result <- GoodmanKruskalGamma(data$ideology_num, data$policy_num,
                                    conf.level = 0.95)

cat(sprintf("Gamma(γ) = %.3f\n", gamma_result[1]))
cat(sprintf("95%% CI: [%.3f, %.3f]\n", gamma_result[2], gamma_result[3]))

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

n_pairs <- n * (n - 1) / 2

# Calculate manually
concordant <- sum(outer(data$ideology_num, data$ideology_num, "<") &
                  outer(data$policy_num, data$policy_num, "<")) +
              sum(outer(data$ideology_num, data$ideology_num, ">") &
                  outer(data$policy_num, data$policy_num, ">"))

discordant <- sum(outer(data$ideology_num, data$ideology_num, "<") &
                  outer(data$policy_num, data$policy_num, ">")) +
              sum(outer(data$ideology_num, data$ideology_num, ">") &
                  outer(data$policy_num, data$policy_num, "<"))

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))

cat("\nNote: Gamma = (C - D) / (C + D), ignoring tied pairs\n")
cat(sprintf("Manual calculation: (%.0f - %.0f) / (%.0f + %.0f) = %.3f\n",
            concordant, discordant, concordant, discordant,
            (concordant - discordant) / (concordant + discordant)))

# === STEP 5: Compare with Kendall's Tau-b ===

cat("\n=== COMPARISON WITH KENDALL'S TAU-B ===\n")

tau_b <- cor(data$ideology_num, data$policy_num, method = "kendall")
cat(sprintf("Kendall's tau-b = %.3f\n", tau_b))
cat(sprintf("Gamma = %.3f\n", gamma_result[1]))
cat(sprintf("Ratio γ/τb = %.2f\n", gamma_result[1] / tau_b))

cat("\nWhy gamma > tau-b:\n")
cat("  Gamma ignores tied pairs(focuses on concordant vs discordant)\n")
cat("  Tau-b adjusts for ties, diluting the association\n")
cat("  Gamma typically 1.2-1.5× larger than tau-b for same data\n")

# === STEP 6: Statistical Significance Testing ===

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

# Asymptotic test using CI
ci_width <- gamma_result[3] - gamma_result[2]
se_approx <- ci_width / (2 * 1.96)
z_stat <- gamma_result[1] / se_approx
p_value <- 2 * pnorm(-abs(z_stat))

cat(sprintf("z-statistic = %.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 7: Effect Size Interpretation ===

cat("\n=== EFFECT SIZE INTERPRETATION ===\n")
cat("Gamma magnitude guidelines:\n")
cat("  |γ| < 0.1: negligible\n")
cat("  0.1 ≤ |γ| < 0.3: small\n")
cat("  0.3 ≤ |γ| < 0.6: moderate\n")
cat("  |γ| ≥ 0.6: large\n\n")

gamma_val <- abs(gamma_result[1])
if (gamma_val < 0.1) {
  strength <- "negligible"
} else if (gamma_val < 0.3) {
  strength <- "small"
} else if (gamma_val < 0.6) {
  strength <- "moderate"
} else {
  strength <- "large"
}

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

# PRE interpretation for gamma
cat("\nPRE interpretation(for concordant vs discordant pairs only):\n")
cat(sprintf("Of pairs that are not tied, %.0f%% more are concordant than discordant\n",
            abs(gamma_result[1]) * 100))

# === STEP 8: Visualization ===

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

# Stacked bar chart
ggplot(data, aes(x = ideology, fill = policy_support)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent) +
  scale_fill_brewer(palette = "RdYlGn", direction = -1) +
  labs(title = "Policy Support by Political Ideology",
       subtitle = sprintf("Gamma = %.2f (%s negative association)",
                         gamma_result[1], strength),
       x = "Political Ideology",
       y = "Proportion",
       fill = "Policy Support") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Mean policy support by ideology
mean_support <- data %>%
  group_by(ideology_num) %>%
  summarise(mean_policy = mean(policy_num),
            se = sd(policy_num) / sqrt(n()))

ggplot(mean_support, aes(x = ideology_num, y = mean_policy)) +
  geom_line(linewidth = 1.2, color = "blue") +
  geom_point(size = 3, color = "blue") +
  geom_errorbar(aes(ymin = mean_policy - 1.96*se,
                    ymax = mean_policy + 1.96*se),
                width = 0.2) +
  labs(title = "Mean Policy Support by Ideology",
       subtitle = "Clear negative monotonic trend",
       x = "Political Ideology(1=Very Liberal to 5=Very Conservative)",
       y = "Mean Policy Support(±95% CI)") +
  theme_classic()

# === APA-STYLE REPORTING ===

cat("\n=== APA-STYLE REPORT ===\n")
cat(sprintf(
"Goodman-Kruskal gamma was computed to assess the symmetric ordinal association
between political ideology(5 ordered categories from very liberal to very
conservative) and environmental policy support(5 ordered categories from strongly
oppose to strongly support) in 180 voters. Gamma was chosen over Kendall's tau-b
because it focuses on the pure concordant vs discordant relationship, ignoring
tied pairs, making it useful for comparing associations across studies with varying
tie structures.

A large negative association was found, γ = %.2f, 95%% CI [%.2f, %.2f], z = %.2f,
p < .001, indicating that more liberal ideology was strongly associated with greater
environmental policy support(negative because liberal=1 and high support=5). The
effect size was %s according to standard interpretation guidelines. Of observation
pairs that were not tied, %.0f%% more were concordant(both variables moved in same
direction) than discordant.

The analysis revealed %.0f%% concordant pairs, %.0f%% discordant pairs, and %.0f%%
tied pairs(on at least one variable). Gamma(%.2f) was larger than Kendall's tau-b
(%.2f) by a factor of %.2f, reflecting gamma's exclusion of ties. These findings
align with political science research showing strong associations between liberal
ideology and environmental support(Dunlap et al., 2016).\n",
gamma_result[1], gamma_result[2], gamma_result[3], z_stat, strength,
abs(gamma_result[1]) * 100,
100 * concordant / n_pairs,
100 * discordant / n_pairs,
100 * (n_pairs - concordant - discordant) / n_pairs,
gamma_result[1], tau_b, abs(gamma_result[1] / tau_b)
))