Phi Coefficient (φ)

# Phi Coefficient: Gender × Voting Preference (2×2 Table)
# Demonstrates phi coefficient for association between two binary variables

library(tidyverse)
library(psych)       # For phi function
library(vcd)         # For assocstats and mosaic plots
library(DescTools)   # For odds ratio with CI

set.seed(2025)
n <- 250

# Simulate realistic data: gender gap in environmental policy support
# Females more likely to support environmental policy
gender <- sample(c("Male", "Female"), n, replace=TRUE, prob=c(0.48, 0.52))

# Policy support depends on gender
support <- character(n)
for (i in 1:n) {
  if (gender[i] == "Male") {
    support[i] <- sample(c("No", "Yes"), 1, prob=c(0.58, 0.42))  # Males: 42% support
  } else {
    support[i] <- sample(c("No", "Yes"), 1, prob=c(0.38, 0.62))  # Females: 62% support
  }
}

# Binary numeric coding (needed for some functions)
gender_num <- ifelse(gender == "Male", 0, 1)
support_num <- ifelse(support == "No", 0, 1)

data <- data.frame(
  id = 1:n,
  gender = gender,
  support = support,
  gender_num = gender_num,
  support_num = support_num
)

cat("=== Political Survey: Gender × Environmental Policy Support ===\n")
head(data, 10)

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

cat("\n=== 2×2 Contingency Table ===\n")
tab <- table(data$gender, data$support)
print(addmargins(tab))

cat("\nTable dimensions:", nrow(tab), "×", ncol(tab), "(2×2 - perfect for phi coefficient)\n")

# Check for zero cells
if (any(tab == 0)) {
  cat("\n⚠ WARNING: Zero cells detected - may affect interpretation\n")
} else {
  cat("\n✓ No zero cells - all cells have observations\n")
}

# === STEP 2: Check Expected Frequencies ===

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

if (all(expected >= 5)) {
  cat("\n✓ All expected frequencies ≥5 - chi-square approximation valid\n")
  cat("   Phi coefficient p-value from chi-square is reliable\n")
} else if (any(expected < 5)) {
  cat("\n⚠ Some expected frequencies <5 - use Fisher's exact test\n")
}

# === STEP 3: Compute Phi Coefficient ===

cat("\n=== Phi Coefficient Analysis ===\n")

# Method 1: Using psych package
phi_val <- phi(tab)
cat(sprintf("φ = %.3f\n", phi_val))

# Method 2: Manual calculation using χ²/n relationship
chi_stat <- chi_result$statistic
phi_manual <- sqrt(chi_stat / n)
if (phi_val < 0) phi_manual <- -phi_manual  # Preserve sign
cat(sprintf("Manual calculation: φ = √(χ²/n) = √(%.2f/%d) = %.3f\n", 
            chi_stat, n, phi_manual))

# Method 3: Via vcd package (comprehensive output)
assoc <- assocstats(tab)
cat("\n=== Association Statistics(vcd package) ===\n")
print(assoc)

# === STEP 4: Significance Test ===

cat("\n=== Chi-square Test of Independence ===\n")
cat(sprintf("χ²(1) = %.2f, p %s\n",
            chi_stat,
            ifelse(chi_result$p.value < 0.001, "< .001",
                   sprintf("= %.4f", chi_result$p.value))))

# Fisher's exact test (for comparison/validation)
fisher_result <- fisher.test(tab)
cat("\n=== Fisher's Exact Test(exact p-value) ===\n")
cat(sprintf("p %s\n",
            ifelse(fisher_result$p.value < 0.001, "< .001",
                   sprintf("= %.4f", fisher_result$p.value))))
cat(sprintf("Odds Ratio = %.2f, 95%% CI [%.2f, %.2f]\n",
            fisher_result$estimate,
            fisher_result$conf.int[1],
            fisher_result$conf.int[2]))

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

cat("\n=== Effect Size Interpretation ===\n")

if (abs(phi_val) < 0.1) {
  strength <- "negligible"
} else if (abs(phi_val) < 0.3) {
  strength <- "small"
} else if (abs(phi_val) < 0.5) {
  strength <- "medium"
} else {
  strength <- "large"
}

cat(sprintf("φ = %.2f is a %s effect(Cohen's benchmarks)\n", phi_val, strength))
cat("Cohen(1988) benchmarks: 0.1=small, 0.3=medium, 0.5=large\n")

# Proportion of variance explained
phi_squared <- phi_val^2
cat(sprintf("\nφ² = %.3f (%.1f%% of variance explained)\n", 
            phi_squared, 100*phi_squared))

# === STEP 6: Complementary Effect Sizes ===

cat("\n=== Complementary Effect Sizes ===\n")

# Cramér's V (should equal |phi| for 2×2 tables)
cramers_v <- sqrt(chi_stat / (n * (min(dim(tab)) - 1)))
cat(sprintf("Cramér's V = %.3f (equals |φ| for 2×2 tables)\n", cramers_v))

# Odds ratio with CI
OR <- OddsRatio(tab, conf.level=0.95)
cat(sprintf("\nOdds Ratio = %.2f, 95%% CI [%.2f, %.2f]\n",
            OR[1], OR[2], OR[3]))
cat("Interpretation: Odds of supporting policy are %.1fx higher for females vs males\n",
    OR[1])

# Risk ratio (relative risk)
risk_female <- tab["Female", "Yes"] / sum(tab["Female", ])
risk_male <- tab["Male", "Yes"] / sum(tab["Male", ])
RR <- risk_female / risk_male
cat(sprintf("\nRelative Risk = %.2f\n", RR))
cat(sprintf("Females are %.1fx more likely to support policy than males\n", RR))

# === STEP 7: Visualizations ===

# Mosaic plot
mosaic(~ gender + support, data=data,
       shade=TRUE,  # Color by Pearson residuals
       main=sprintf("Gender × Policy Support\nφ = %.2f (%s effect)", phi_val, strength),
       labeling=labeling_border(rot_labels=c(0, 0, 0, 0)))

# Grouped bar plot
ggplot(data, aes(x=gender, fill=support)) +
  geom_bar(position="dodge") +
  scale_fill_manual(values=c("No"="coral", "Yes"="steelblue")) +
  labs(title="Policy Support by Gender",
       subtitle=sprintf("φ = %.2f, p < .001 (%s association)", phi_val, strength),
       x="Gender",
       y="Count",
       fill="Support Policy") +
  theme_classic() +
  theme(text=element_text(size=12))

# Stacked proportions
ggplot(data, aes(x=gender, fill=support)) +
  geom_bar(position="fill") +
  scale_fill_manual(values=c("No"="coral", "Yes"="steelblue")) +
  scale_y_continuous(labels=scales::percent) +
  labs(title="Policy Support Distribution by Gender",
       subtitle=sprintf("Females: %.0f%% support | Males: %.0f%% support",
                        100*risk_female, 100*risk_male),
       x="Gender",
       y="Proportion",
       fill="Support Policy") +
  theme_classic()

# === STEP 8: Detailed Cross-tabulation ===

cat("\n=== Detailed Cross-Tabulation ===\n")
prop_table <- prop.table(tab, margin=1) * 100  # Row percentages
cat("\nRow percentages(% within gender):\n")
print(round(prop_table, 1))

cat("\nInterpretation:\n")
cat(sprintf("  • %.1f%% of females support the policy\n", prop_table["Female", "Yes"]))
cat(sprintf("  • %.1f%% of males support the policy\n", prop_table["Male", "Yes"]))
cat(sprintf("  • Gender gap: %.1f percentage points\n", 
            prop_table["Female", "Yes"] - prop_table["Male", "Yes"]))

# === STEP 9: Confidence Interval for Phi ===

cat("\n=== Bootstrap 95% CI for Phi Coefficient ===\n")

# Bootstrap CI
set.seed(2025)
B <- 1000
phi_boot <- numeric(B)

for (b in 1:B) {
  # Resample with replacement
  idx <- sample(1:n, n, replace=TRUE)
  boot_tab <- table(data$gender[idx], data$support[idx])
  if (all(dim(boot_tab) == c(2, 2))) {  # Ensure 2×2 table
    phi_boot[b] <- phi(boot_tab)
  } else {
    phi_boot[b] <- NA
  }
}

phi_boot <- na.omit(phi_boot)
phi_ci <- quantile(phi_boot, c(0.025, 0.975))

cat(sprintf("φ = %.3f, 95%% Bootstrap CI [%.3f, %.3f]\n",
            phi_val, phi_ci[1], phi_ci[2]))

if (sign(phi_ci[1]) == sign(phi_ci[2])) {
  cat("✓ CI does not cross zero - association is significant\n")
} else {
  cat("⚠ CI crosses zero - association not significant\n")
}

# === APA-Style Reporting ===

cat("\n=== APA-Style Report ===\n")
cat(sprintf(
"A phi coefficient was computed to assess the association between gender and 
support for a new environmental policy in a sample of 250 registered voters. 
The 2×2 contingency table showed that %.1f%% of females(n=%d) supported the 
policy compared to %.1f%% of males(n=%d). There was a significant positive 
association, φ = %.2f, 95%% CI [%.2f, %.2f], χ²(1) = %.2f, p < .001, indicating 
that gender and policy support were related. The effect size was %s according to 
Cohen(1988) guidelines. The phi coefficient accounted for %.1f%% of the variance 
in policy support(φ² = %.3f). The odds of supporting the policy were %.1f times 
higher for females compared to males(OR = %.2f, 95%% CI [%.2f, %.2f]). These 
findings suggest a gender gap in environmental policy preferences, with females 
showing greater support. Assumptions were met: both variables were genuinely 
binary, observations were independent(simple random sample), and all expected 
cell frequencies exceeded 5 (minimum = %.1f).\n",
100*risk_female, sum(tab["Female", ]),
100*risk_male, sum(tab["Male", ]),
phi_val, phi_ci[1], phi_ci[2],
chi_stat, strength,
100*phi_squared, phi_squared,
OR[1], OR[1], OR[2], OR[3],
min(expected)
))