使用 R，将多个卡方列联表测试应用于分组数据框，并添加包含测试 p 值的新列

我有一个类似于下面的示例的数据框（这是我的实际数据框的一小部分摘录）。

frequencies <- data.frame(sex=c("female", "female", "male", "male", "female", "female", "male", "male", "female", "female", "male", "male", "female", "female", "male", "male"),
                      ecotype=c("Crab", "Wave", "Crab", "Wave", "Crab", "Wave", "Crab", "Wave", "Crab", "Wave", "Crab", "Wave", "Crab", "Wave", "Crab", "Wave"),
                      contig_ID=c("Contig100169_2367", "Contig100169_2367", "Contig100169_2367", "Contig100169_2367", "Contig100169_2367", "Contig100169_2367", "Contig100169_2367", "Contig100169_2367", 
                                  "Contig100169_2481", "Contig100169_2481", "Contig100169_2481", "Contig100169_2481", "Contig100169_2481", "Contig100169_2481", "Contig100169_2481", "Contig100169_2481"),
                      allele=c("p", "p", "p", "p", "q", "q", "q", "q", "p", "p", "p", "p", "q", "q", "q", "q"),
                      frequency=c(157, 98, 140, 65, 29, 8, 26, 9, 182, 108, 147, 80, 46, 4, 49, 4))

我想对“contig_ID”和“ecotype”的每个组合进行单独的卡方应急测试，测试“性别”和“等位基因”之间的关联。然后，我想将这些结果总结在一个表中，其中包括“contig_ID”和“ecotype”每种组合的 p 值。例如，从给出的示例表中，我期望有 4 个 p 值的结果表，如下例所示。

results <- data.frame(ecotype=c("Crab", "Wave", "Crab", "Wave"),
                  contig_ID=c("Contig100169_2367", "Contig100169_2367", "Contig100169_2481", "Contig100169_2481"),
                  pvalue=c("pval", "pval", "pval", "pval"))

或者，仅向原始表添加 p 值列也可以，每个组合的 p 值仅在所有相关行中重复。

我一直在尝试使用诸如lapply() and summarise()结合chisq.test()实现这一目标，但到目前为止还没有运气。我也尝试过使用类似的方法：表中每一行的 R 卡方检验（3x2 列联表） https://stackoverflow.com/questions/34232869/r-chi-squared-test-3x2-contingency-table-for-each-row-in-a-table，但也无法完成这项工作。

我们可以将contig_ID and ecotype列并创建一个嵌套数据框，并将数据转换为矩阵，如下所示。

library(tidyverse)

frequencies2 <- frequencies %>%
  group_by(contig_ID, ecotype) %>%
  nest() %>%
  mutate(M = map(data, function(dat){
    dat2 <- dat %>% spread(sex, frequency)
    M <- as.matrix(dat2[, -1])
    row.names(M) <- dat2$allele
    return(M)
  }))

如果我们看一下第一个元素M列中，我们会发现每组的数据都转换为矩阵。

frequencies2$M[[1]]
#   female male
# p    157  140
# q     29   26

从这里，我们可以应用chisq.test到每个矩阵并提取 p 值。frequencies3是最终的输出。

frequencies3 <- frequencies2 %>%
  mutate(pvalue = map_dbl(M, ~chisq.test(.x)$p.value)) %>%
  select(-data, -M) %>%
  ungroup()
frequencies3
# # A tibble: 4 x 3
#   contig_ID         ecotype pvalue
#   <fct>             <fct>    <dbl>
# 1 Contig100169_2367 Crab     1.00 
# 2 Contig100169_2367 Wave     0.434
# 3 Contig100169_2481 Crab     0.284
# 4 Contig100169_2481 Wave     0.958