Cramer's V

# Cramér's V: Gender × Political Party Association
# 2×3 contingency table with small-medium effect size

library(tidyverse)
library(vcd)           # For assocstats (Cramér's V)
library(DescTools)     # For CramerV with CI
library(effectsize)    # For cramers_v with CI

set.seed(2025)
n <- 300

# Simulate gender-party association (women lean Democrat, men lean Republican)
# Realistic probabilities based on polling data
gender <- sample(c("Male", "Female"), n, replace=TRUE, prob=c(0.48, 0.52))

party <- character(n)
for (i in 1:n) {
  if (gender[i] == "Female") {
    party[i] <- sample(c("Democrat", "Republican", "Independent"), 1, 
                       prob=c(0.48, 0.28, 0.24))  # Women lean Democrat
  } else {
    party[i] <- sample(c("Democrat", "Republican", "Independent"), 1, 
                       prob=c(0.32, 0.42, 0.26))  # Men lean Republican
  }
}

data <- data.frame(
  respondent_id = 1:n,
  gender = factor(gender, levels=c("Male", "Female")),
  party = factor(party, levels=c("Democrat", "Republican", "Independent"))
)

head(data, 10)

# === STEP 1: Create Contingency Table ===

cat("=== Contingency Table: Gender × Political Party ===\n")
contingency <- table(data$gender, data$party)
print(contingency)

cat("\n=== Row Percentages(% within each gender) ===\n")
row_pct <- prop.table(contingency, margin=1) * 100
print(round(row_pct, 1))

cat("\n=== Column Percentages(% within each party) ===\n")
col_pct <- prop.table(contingency, margin=2) * 100
print(round(col_pct, 1))

# === STEP 2: Check Chi-Square Assumptions ===

chi_result <- chisq.test(contingency)
cat("\n=== Expected Frequencies ===\n")
print(round(chi_result$expected, 2))

min_expected <- min(chi_result$expected)
cat(sprintf("\nMinimum expected frequency: %.2f\n", min_expected))
cat(sprintf("All cells ≥5? %s\n", ifelse(min_expected >= 5, "YES ✓", "NO ✗")))

pct_cells_below_5 <- 100 * sum(chi_result$expected < 5) / length(chi_result$expected)
cat(sprintf("Cells with expected freq <5: %.0f%% (should be <20%%)\n", pct_cells_below_5))

if (pct_cells_below_5 < 20) {
  cat("Chi-square assumptions met ✓\n")
} else {
  cat("WARNING: Chi-square assumptions violated - consider Fisher's exact test\n")
}

# === STEP 3: Compute Chi-Square Test ===

cat("\n=== Chi-Square Test of Independence ===\n")
print(chi_result)

cat(sprintf("\nχ²(%d) = %.2f, p %s\n", 
            chi_result$parameter,
            chi_result$statistic,
            ifelse(chi_result$p.value < 0.001, "< .001", 
                   sprintf("= %.3f", chi_result$p.value))))

if (chi_result$p.value < 0.05) {
  cat("Conclusion: Significant association between gender and party affiliation\n")
} else {
  cat("Conclusion: No significant association detected\n")
}

# === STEP 4: Compute Cramér's V ===

# Method 1: vcd package (classic)
assoc_stats <- assocstats(contingency)
cat("\n=== Association Measures(vcd::assocstats) ===\n")
print(assoc_stats)

# Method 2: DescTools (with CI)
v_ci <- CramerV(contingency, conf.level=0.95)
cat(sprintf("\nCramér's V = %.3f, 95%% CI [%.3f, %.3f]\n", 
            v_ci[1], attr(v_ci, "lwr.ci"), attr(v_ci, "upr.ci")))

# Method 3: effectsize package (recommended for modern reporting)
v_effect <- cramers_v(contingency, ci=0.95)
cat("\n=== Cramér's V(effectsize package) ===\n")
print(v_effect)

v_value <- as.numeric(v_effect$Cramers_v)

# === STEP 5: Interpret Effect Size ===

cat("\n=== Effect Size Interpretation ===\n")
cat("For df = min(r-1, c-1) = min(2-1, 3-1) = 1:\n")
cat("Cohen(1988) benchmarks: small=0.10, medium=0.30, large=0.50\n\n")

df_cramers <- min(nrow(contingency)-1, ncol(contingency)-1)

if (df_cramers == 1) {
  # Use Cohen's standard benchmarks for df=1
  if (v_value < 0.10) {
    effect_interp <- "negligible"
  } else if (v_value < 0.30) {
    effect_interp <- "small"
  } else if (v_value < 0.50) {
    effect_interp <- "medium"
  } else {
    effect_interp <- "large"
  }
} else {
  # Adjust benchmarks upward for df>1
  if (v_value < 0.07) {
    effect_interp <- "negligible"
  } else if (v_value < 0.21) {
    effect_interp <- "small"
  } else if (v_value < 0.35) {
    effect_interp <- "medium"
  } else {
    effect_interp <- "large"
  }
}

cat(sprintf("Cramér's V = %.3f → %s effect size\n", v_value, effect_interp))
cat(sprintf("Variance explained(approximate): %.1f%%\n", v_value^2 * 100))

# === STEP 6: Standardized Residuals (which cells drive association?) ===

cat("\n=== Standardized Residuals(|z| > 2 indicates significant cell) ===\n")
std_resid <- chi_result$stdres
print(round(std_resid, 2))

cat("\nInterpretation: Positive residuals = more obs than expected;")
cat(" Negative = fewer than expected\n")

# Identify significant cells
significant_cells <- which(abs(std_resid) > 2, arr.ind=TRUE)
if (nrow(significant_cells) > 0) {
  cat("\nCells with |standardized residual| > 2:\n")
  for (i in 1:nrow(significant_cells)) {
    row_idx <- significant_cells[i, 1]
    col_idx <- significant_cells[i, 2]
    cat(sprintf("  %s × %s: z = %.2f\n", 
                rownames(contingency)[row_idx],
                colnames(contingency)[col_idx],
                std_resid[row_idx, col_idx]))
  }
}

# === STEP 7: Visualizations ===

# Mosaic plot (area proportional to frequency)
par(mfrow=c(1,1))
mosaic(contingency, 
       shade=TRUE, 
       legend=TRUE,
       main="Mosaic Plot: Gender × Party\n(Blue=more than expected, Red=fewer)")

# Grouped bar plot
library(ggplot2)
data_summary <- data %>%
  count(gender, party) %>%
  group_by(gender) %>%
  mutate(pct = n / sum(n) * 100)

ggplot(data_summary, aes(x=gender, y=pct, fill=party)) +
  geom_bar(stat="identity", position="dodge", color="black") +
  geom_text(aes(label=sprintf("%.1f%%", pct)), 
            position=position_dodge(width=0.9), vjust=-0.5, size=3) +
  scale_fill_manual(values=c("Democrat"="#0015BC", "Republican"="#E81B23", 
                             "Independent"="#808080")) +
  labs(title="Political Party Affiliation by Gender",
       subtitle=sprintf("Cramér's V = %.3f (%s effect)", v_value, effect_interp),
       x="Gender", y="Percentage", fill="Party") +
  theme_classic() +
  theme(legend.position="bottom")

# Heatmap of observed frequencies
library(pheatmap)
pheatmap(contingency, 
         display_numbers=TRUE, 
         cluster_rows=FALSE, 
         cluster_cols=FALSE,
         main="Frequency Heatmap: Gender × Party",
         color=colorRampPalette(c("white", "steelblue"))(50))

# === STEP 8: Bias-Corrected Cramér's V (for small samples) ===

cat("\n=== Bias-Corrected Cramér's V ===\n")
# Bias correction: V_corrected = sqrt(max(0, φ² - ((r-1)(c-1))/(n-1)))
phi_squared <- chi_result$statistic / n
r <- nrow(contingency)
c <- ncol(contingency)
bias_term <- ((r - 1) * (c - 1)) / (n - 1)
v_corrected <- sqrt(max(0, phi_squared - bias_term))

cat(sprintf("Uncorrected V = %.3f\n", v_value))
cat(sprintf("Bias-corrected V = %.3f\n", v_corrected))
cat(sprintf("Difference = %.4f (negligible for n=%d)\n", v_value - v_corrected, n))
cat("Note: Bias correction more important for small n or large tables\n")

# === STEP 9: Power Analysis ===

cat("\n=== Post-hoc Power Analysis ===\n")
library(pwr)
# For chi-square, effect size w = V when df=1; adjust otherwise
effect_size_w <- v_value
power_result <- pwr.chisq.test(w=effect_size_w, N=n, df=(r-1)*(c-1), sig.level=0.05)
cat(sprintf("Achieved power = %.2f (for detecting V=%.2f at α=.05)\n", 
            power_result$power, v_value))

if (power_result$power < 0.80) {
  cat("WARNING: Power below 0.80 - results may be underpowered\n")
}

# === APA-Style Reporting ===

cat("\n=== APA-Style Report ===\n")
cat(sprintf(
"A chi-square test of independence was conducted to examine the association 
between gender and political party affiliation in a sample of %d US voters. 
All expected cell frequencies exceeded 5, meeting chi-square assumptions. 
There was a significant association between gender and party affiliation, 
χ²(%d, N = %d) = %.2f, p %s. Cramér's V was computed as an effect size 
measure, V = %.2f, 95%% CI [%.2f, %.2f], indicating a %s association. 
Examination of row percentages revealed that women were more likely to 
identify as Democrat(%.0f%%) compared to men(%.0f%%), while men were 
more likely to identify as Republican(%.0f%%) compared to women(%.0f%%). 
Standardized residuals identified Female×Democrat(z = %.2f) and 
Male×Republican(z = %.2f) as cells contributing most strongly to the 
association. These findings align with documented gender gaps in US political 
party affiliation(Pew Research Center, 2023).\n",
n,
chi_result$parameter,
n,
chi_result$statistic,
ifelse(chi_result$p.value < 0.001, "< .001", sprintf("= %.3f", chi_result$p.value)),
v_value,
attr(v_ci, "lwr.ci"),
attr(v_ci, "upr.ci"),
effect_interp,
row_pct["Female", "Democrat"],
row_pct["Male", "Democrat"],
row_pct["Male", "Republican"],
row_pct["Female", "Republican"],
std_resid["Female", "Democrat"],
std_resid["Male", "Republican"]
))