Contingency Coefficient (C)

# Contingency Coefficient: Smoking × Lung Disease
# 3×3 table demonstrating C_max interpretation

library(tidyverse)
library(vcd)           # For assocstats (includes C)
library(DescTools)     # For ContCoef with CI

set.seed(2025)
n <- 250

# Simulate smoking-lung disease association
smoking_categories <- c("Never", "Former", "Current")
disease_categories <- c("None/Mild", "Moderate", "Severe")

smoking <- sample(smoking_categories, n, replace=TRUE, prob=c(0.40, 0.35, 0.25))

# Lung disease severity depends strongly on smoking
disease <- character(n)
for (i in 1:n) {
  if (smoking[i] == "Never") {
    disease[i] <- sample(disease_categories, 1, prob=c(0.70, 0.25, 0.05))
  } else if (smoking[i] == "Former") {
    disease[i] <- sample(disease_categories, 1, prob=c(0.45, 0.40, 0.15))
  } else {  # Current
    disease[i] <- sample(disease_categories, 1, prob=c(0.20, 0.45, 0.35))
  }
}

data <- data.frame(
  participant_id = 1:n,
  smoking = factor(smoking, levels=smoking_categories),
  disease = factor(disease, levels=disease_categories)
)

head(data, 10)

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

cat("=== Contingency Table: Smoking Status × Lung Disease ===\n")
contingency <- table(data$smoking, data$disease)
print(contingency)

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

# === STEP 2: Chi-Square Test ===

chi_result <- chisq.test(contingency)
cat("\n=== Chi-Square Test of Independence ===\n")
print(chi_result)

cat("\nExpected Frequencies:\n")
print(round(chi_result$expected, 2))

min_expected <- min(chi_result$expected)
cat(sprintf("\nMinimum expected frequency: %.2f (should be ≥5)\n", min_expected))

if (min_expected >= 5 & sum(chi_result$expected < 5) / length(chi_result$expected) < 0.20) {
  cat("Chi-square assumptions satisfied ✓\n")
} else {
  cat("WARNING: Chi-square assumptions may be violated\n")
}

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

# === STEP 3: Calculate Contingency Coefficient C ===

# Manual calculation
chi2 <- chi_result$statistic
c_raw <- sqrt(chi2 / (chi2 + n))

cat("\n=== Contingency Coefficient Calculation ===\n")
cat(sprintf("C = √(χ²/(χ²+n)) = √(%.2f/(%.2f+%d)) = %.3f\n", 
            chi2, chi2, n, c_raw))

# Calculate C_max
r <- nrow(contingency)
c_cols <- ncol(contingency)
min_dim <- min(r, c_cols)
c_max <- sqrt((min_dim - 1) / min_dim)

cat(sprintf("\nC_max = √((min(r,c)-1)/min(r,c)) = √((%d-1)/%d) = %.3f\n",
            min_dim, min_dim, c_max))
cat(sprintf("\nFor %d×%d table: C can range from 0 to %.3f (not 0 to 1)\n",
            r, c_cols, c_max))

# Adjusted Contingency Coefficient C*
c_adjusted <- c_raw / c_max
cat(sprintf("\nAdjusted C* = C/C_max = %.3f/%.3f = %.3f\n", 
            c_raw, c_max, c_adjusted))
cat("C* rescales C to 0-1 range for interpretability\n")

# Using DescTools (if available)
if (requireNamespace("DescTools", quietly=TRUE)) {
  c_ci <- DescTools::ContCoef(contingency, conf.level=0.95)
  cat(sprintf("\nContingency Coefficient with CI: C = %.3f, 95%% CI [%.3f, %.3f]\n",
              c_ci[1], attr(c_ci, "lwr.ci"), attr(c_ci, "upr.ci")))
}

# === STEP 4: Compare with Cramér's V ===

cat("\n=== Comparison: C vs Cramér's V ===\n")
assoc_stats <- assocstats(contingency)
print(assoc_stats)

v_value <- assoc_stats$cramer
cat(sprintf("\nContingency Coefficient C = %.3f (max = %.3f)\n", c_raw, c_max))
cat(sprintf("Adjusted C* = %.3f\n", c_adjusted))
cat(sprintf("Cramér's V = %.3f (always ranges 0-1)\n", v_value))
cat("\nNote: Cramér's V preferred over C for standardized reporting\n")
cat("V does not have the C_max limitation\n")

# === STEP 5: Effect Size Interpretation ===

cat("\n=== Effect Size Interpretation ===\n")
cat("Interpreting C is complex due to C_max variation by table size\n")
cat("Better to interpret adjusted C* on 0-1 scale:\n\n")

if (c_adjusted < 0.30) {
  c_interp <- "small"
} else if (c_adjusted < 0.50) {
  c_interp <- "medium"
} else {
  c_interp <- "large"
}

cat(sprintf("C* = %.3f → %s effect(using 0.30/0.50 thresholds)\n", 
            c_adjusted, c_interp))
cat(sprintf("Relative to maximum: C is %.0f%% of C_max\n", 
            100 * c_raw / c_max))

# For comparison, interpret Cramér's V (df=2 benchmarks: 0.07/0.21/0.35)
df <- min(r-1, c_cols-1)
cat(sprintf("\nFor comparison, Cramér's V = %.3f\n", v_value))
if (df == 2) {
  if (v_value < 0.21) {
    v_interp <- "small"
  } else if (v_value < 0.35) {
    v_interp <- "medium"
  } else {
    v_interp <- "large"
  }
  cat(sprintf("With df=2, V = %.3f → %s effect\n", v_value, v_interp))
}

# === STEP 6: Standardized Residuals ===

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

cat("\nCells with |standardized residual| > 2:\n")
for (i in 1:nrow(std_resid)) {
  for (j in 1:ncol(std_resid)) {
    if (abs(std_resid[i,j]) > 2) {
      cat(sprintf("  %s × %s: z = %.2f\n",
                  rownames(contingency)[i],
                  colnames(contingency)[j],
                  std_resid[i,j]))
    }
  }
}

# === STEP 7: Visualizations ===

# Mosaic plot
par(mfrow=c(1,1))
mosaic(contingency, shade=TRUE, legend=TRUE,
       main=sprintf("Smoking × Lung Disease\nC = %.2f (max = %.2f), C* = %.2f",
                    c_raw, c_max, c_adjusted))

# Heatmap
library(pheatmap)
pheatmap(contingency, 
         display_numbers=TRUE, 
         cluster_rows=FALSE, 
         cluster_cols=FALSE,
         main=sprintf("Smoking × Disease(C* = %.2f, %s effect)", 
                      c_adjusted, c_interp),
         color=colorRampPalette(c("white", "orange", "red"))(50))

# Stacked bar plot
data_long <- data %>%
  count(smoking, disease) %>%
  group_by(smoking) %>%
  mutate(pct = n / sum(n) * 100)

ggplot(data_long, aes(x=smoking, y=pct, fill=disease)) +
  geom_bar(stat="identity", color="black") +
  geom_text(aes(label=sprintf("%.0f%%", pct)),
            position=position_stack(vjust=0.5), color="white", size=4) +
  scale_fill_manual(values=c("None/Mild"="#2ECC71", 
                             "Moderate"="#F39C12", 
                             "Severe"="#E74C3C")) +
  labs(title="Lung Disease Severity by Smoking Status",
       subtitle=sprintf("C = %.2f (C* = %.2f after adjustment), V = %.2f",
                        c_raw, c_adjusted, v_value),
       x="Smoking Status", y="Percentage", fill="Disease Severity") +
  theme_classic()

# === STEP 8: Table of C_max Values ===

cat("\n=== Reference: C_max for Different Table Dimensions ===\n")
table_sizes <- data.frame(
  Table = c("2×2", "3×3", "4×4", "5×5", "3×4", "4×5", "2×5"),
  min_dim = c(2, 3, 4, 5, 3, 4, 2),
  C_max = c(sqrt(1/2), sqrt(2/3), sqrt(3/4), sqrt(4/5), 
            sqrt(2/3), sqrt(3/4), sqrt(1/2))
)
table_sizes$C_max <- round(table_sizes$C_max, 3)
print(table_sizes)

cat("\nNote: C_max = √((min(r,c)-1)/min(r,c)) varies by table dimensions\n")
cat("This is why Cramér's V(always 0-1) is preferred over C\n")

# === APA-Style Reporting ===

cat("\n=== APA-Style Report ===\n")
cat(sprintf(
"A chi-square test of independence examined the association between smoking 
status and lung disease severity in %d adults aged 40-70. All expected cell 
frequencies exceeded 5, satisfying chi-square assumptions. There was a 
significant association, χ²(%d, N = %d) = %.2f, p < .001. The Contingency 
Coefficient was C = %.2f, which must be interpreted relative to its maximum 
possible value of C_max = %.2f for a 3×3 table. The adjusted coefficient 
C* = C/C_max = %.2f indicates a %s association on the 0-1 scale. For 
comparison, Cramér's V = %.2f (which always ranges 0-1), also indicating a 
%s effect. Examination of row percentages revealed that current smokers had 
higher rates of severe lung disease(%.0f%%) compared to former smokers 
(%.0f%%) and never smokers(%.0f%%). Standardized residuals identified 
Current Smokers×Severe Disease(z = %.2f) and Never Smokers×None/Mild 
(z = %.2f) as cells contributing most strongly to the association. These 
findings align with epidemiological evidence linking smoking to lung disease 
severity. Note: Cramér's V is generally preferred over the Contingency 
Coefficient for effect size reporting due to its consistent 0-1 range across 
all table dimensions.\n",
n,
chi_result$parameter,
n,
chi_result$statistic,
c_raw,
c_max,
c_adjusted,
c_interp,
v_value,
v_interp,
row_pct["Current", "Severe"],
row_pct["Former", "Severe"],
row_pct["Never", "Severe"],
std_resid["Current", "Severe"],
std_resid["Never", "None/Mild"]
))