R语言新手求助：批量计算多组fisher.test p值的实现方案

Hey there! Let's walk through this step by step—since you're new to R, I'll keep things clear and avoid jargon where I can.

First, let's align on your goal: you want to compare, for every pair of NYC boroughs and each year, whether the frequency of both boroughs having an event (marked "1") differs between type P and Q using Fisher's exact test. You already melted your data into long format, so we'll build from there.

Step 1: Prep Your Data for Pairwise Comparisons #

First, we need to reshape your melted long-format data into a wide format where each row has a year, type (P/Q), and event values for all boroughs. This makes it easy to compare pairs of boroughs side-by-side.

Assuming your melted data has columns year, borough, type, and event (0/1), use either reshape2 (since you already used melt) or the modern tidyverse approach:

# Using reshape2
library(reshape2)
wide_df <- dcast(your_melted_data, year + type ~ borough, value.var = "event")

# OR using tidyverse (more intuitive for many users)
library(tidyverse)
wide_df <- your_melted_data %>%
  pivot_wider(names_from = borough, values_from = event)

Step 2: Generate All Unique Borough Pairs #

We need every unique pair of boroughs (e.g., Manhattan & Brooklyn, Queens & Bronx). Use combn to create these pairs without duplicates:

# Get unique borough names from your data
all_boroughs <- unique(your_melted_data$borough)
# Generate all 2-borough combinations (no reverse pairs like B&A if we already have A&B)
borough_pairs <- t(combn(all_boroughs, 2))

Step 3: Batch Calculate Fisher's Exact Test (Beginner-Friendly Loop) #

A simple loop will let you iterate over every borough pair and year, calculate key stats, and run the test. We'll also add checks to fix the invalid 2x2 table error you encountered!

# Initialize an empty data frame to store results
results_df <- data.frame(
  borough1 = character(),
  borough2 = character(),
  year = integer(),
  p_both_event = integer(),
  p_total = integer(),
  p_proportion = numeric(),
  q_both_event = integer(),
  q_total = integer(),
  q_proportion = numeric(),
  odds_ratio = numeric(),
  p_value = numeric(),
  stringsAsFactors = FALSE
)

# Loop through each borough pair
for (i in 1:nrow(borough_pairs)) {
  b1 <- borough_pairs[i, 1]
  b2 <- borough_pairs[i, 2]
  
  # Loop through each year
  for (yr in unique(wide_df$year)) {
    # Filter data for the current year
    year_subset <- wide_df[wide_df$year == yr, ]
    
    # Extract P and Q type data for the two boroughs
    p_data <- year_subset[year_subset$type == "P", c(b1, b2)]
    q_data <- year_subset[year_subset$type == "Q", c(b1, b2)]
    
    # Calculate counts for P type
    p_both <- sum(p_data[[b1]] == 1 & p_data[[b2]] == 1, na.rm = TRUE)
    p_total <- nrow(p_data)
    p_not_both <- p_total - p_both
    p_prop <- if(p_total > 0) p_both / p_total else NA
    
    # Calculate counts for Q type
    q_both <- sum(q_data[[b1]] == 1 & q_data[[b2]] == 1, na.rm = TRUE)
    q_total <- nrow(q_data)
    q_not_both <- q_total - q_both
    q_prop <- if(q_total > 0) q_both / q_total else NA
    
    # Build the 2x2 contingency table
    contingency_table <- matrix(
      c(p_both, q_both, p_not_both, q_not_both),
      nrow = 2,
      dimnames = list(Type = c("P", "Q"), Both_Event = c("Yes", "No"))
    )
    
    # Check if the table is valid (fixes your error!)
    # Fisher's test needs at least 2 non-zero cells, no empty rows/columns
    is_valid <- sum(contingency_table > 0) >= 2 && 
      all(rowSums(contingency_table) > 0) && 
      all(colSums(contingency_table) > 0)
    
    # Run test only if table is valid
    if (is_valid) {
      ft <- fisher.test(contingency_table)
      p_val <- ft$p.value
      odds_ratio <- ft$estimate[[1]]
    } else {
      p_val <- NA
      odds_ratio <- NA
      warning(paste("Skipping invalid table for", b1, "&", b2, "in", yr))
    }
    
    # Add results to the data frame
    results_df <- rbind(results_df, data.frame(
      borough1 = b1,
      borough2 = b2,
      year = yr,
      p_both_event = p_both,
      p_total = p_total,
      p_proportion = p_prop,
      q_both_event = q_both,
      q_total = q_total,
      q_proportion = q_prop,
      odds_ratio = odds_ratio,
      p_value = p_val
    ))
  }
}

Step 4: Tidyverse Alternative (More Efficient) #

If you want to use modern R syntax, here's a pipe-based version that avoids explicit loops:

library(tidyverse)
library(broom) # For cleaning test results

# Create all combinations of borough pairs and years (no duplicates)
all_combinations <- expand_grid(
  borough1 = all_boroughs,
  borough2 = all_boroughs
) %>%
  filter(borough1 < borough2) %>% # Skip reverse pairs like B&A
  expand_grid(year = unique(your_melted_data$year))

# Calculate stats and run tests
results_tidy <- all_combinations %>%
  rowwise() %>%
  mutate(
    # Get P and Q data for the pair and year
    p_subset = list(your_melted_data %>%
                      filter(year == !!year, type == "P", borough %in% c(borough1, borough2)) %>%
                      pivot_wider(names_from = borough, values_from = event)),
    q_subset = list(your_melted_data %>%
                      filter(year == !!year, type == "Q", borough %in% c(borough1, borough2)) %>%
                      pivot_wider(names_from = borough, values_from = event)),
    # Calculate counts
    p_both = sum(p_subset[[1]][[borough1]] == 1 & p_subset[[1]][[borough2]] == 1, na.rm = TRUE),
    p_total = nrow(p_subset[[1]]),
    p_prop = if(p_total > 0) p_both / p_total else NA,
    q_both = sum(q_subset[[1]][[borough1]] == 1 & q_subset[[1]][[borough2]] == 1, na.rm = TRUE),
    q_total = nrow(q_subset[[1]]),
    q_prop = if(q_total > 0) q_both / q_total else NA,
    # Build contingency table and check validity
    contingency_table = list(matrix(c(p_both, q_both, p_total - p_both, q_total - q_both), nrow = 2)),
    is_valid = sum(contingency_table[[1]] > 0) >= 2 && all(rowSums(contingency_table[[1]]) > 0) && all(colSums(contingency_table[[1]]) > 0),
    # Run test and extract results
    test_result = list(if(is_valid) fisher.test(contingency_table[[1]]) else NULL),
    p_value = if(is_valid) test_result[[1]]$p.value else NA,
    odds_ratio = if(is_valid) test_result[[1]]$estimate[[1]] else NA
  ) %>%
  select(-p_subset, -q_subset, -contingency_table, -is_valid, -test_result)

Fixing the Invalid 2x2 Table Error #

The error you saw occurs when your contingency table has empty rows/columns (e.g., no events in Q type for a pair, or all observations are "both events"). The is_valid check in both versions skips these cases and marks the p-value as NA, so your code won't crash. You'll get a warning telling you which pair/year caused the issue, so you can investigate if needed.

Example Output #

For your sample case (2012, Manhattan & Brooklyn), the results would look like this in results_df:

内容的提问来源于stack exchange，提问作者Sylvia Rodriguez