使用R对分组数据框执行多组卡方列联表检验并添加p值列

First, let's align on your data structure: your dataframe has repeated entries for each combination of contig_ID, sex, and ecotype (either as pre-calculated frequencies or individual-level observations). We’ll group by each contig, build a 2x2 contingency table of sex vs ecotype, run a chi-squared test, and append the resulting p-value to every row of your original data.

Step 1: Prepare Your Data #

If your data is at the individual level (each row represents one organism), start by aggregating to get frequency counts. Skip this step if you already have a count column.

# Load required libraries
library(dplyr)
library(tidyr)

# Example individual-level data (replace with your actual data)
set.seed(123)
individual_data <- expand.grid(
  sex = c("female", "male"),
  ecotype = c("Crab", "Wave"),
  contig_ID = paste0("Contig", 1:5)
) %>% 
  slice(rep(1:n(), sample(10:50, n(), replace = TRUE))) %>% 
  select(sex, ecotype, contig_ID)

# Aggregate to get counts per contig, sex, ecotype
frequencies <- individual_data %>% 
  group_by(contig_ID, sex, ecotype) %>% 
  summarise(count = n(), .groups = "drop")

Step 2: Define a Function to Calculate Chi-Squared P-Values #

This function takes a subset of data for a single contig, reshapes it into a 2x2 table, runs the chi-squared test (with continuity correction, recommended for 2x2 tables), and returns the p-value.

get_chi_p_value <- function(group_data) {
  # Reshape to wide format to create the contingency table
  cont_table <- group_data %>%
    pivot_wider(names_from = ecotype, values_from = count) %>%
    select(-sex) %>%
    as.matrix()
  
  # Run chi-squared test with continuity correction
  chi_result <- chisq.test(cont_table, correct = TRUE)
  
  # Return the p-value
  return(chi_result$p.value)
}

Step 3: Apply the Function and Add P-Values to Your Dataframe #

Use dplyr grouping to apply the function to each contig, then add a new column with the p-value to every row of your original data.

# Add p-value column to the frequencies dataframe
frequencies_with_p <- frequencies %>%
  group_by(contig_ID) %>%
  mutate(chi_p_value = get_chi_p_value(cur_data())) %>%
  ungroup()

# View the result
head(frequencies_with_p)

Key Notes: #

• Small Sample Sizes: If you get warnings about low expected frequencies, replace chisq.test with fisher.test(cont_table)$p.value in the function for Fisher's exact test, which is more robust for small counts.
• Interpretation: A p-value < 0.05 typically indicates a significant association between sex and ecotype for that contig.

Example Output: #

# A tibble: 6 × 5
  contig_ID sex    ecotype count chi_p_value
  <chr>     <fct>  <fct>   <int>       <dbl>
1 Contig1   female Crab       38       0.248
2 Contig1   female Wave       21       0.248
3 Contig1   male   Crab       20       0.248
4 Contig1   male   Wave       14       0.248
5 Contig2   female Crab       10       0.689
6 Contig2   female Wave       17       0.689

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