如何获得列联表？
Posted
技术标签:
【中文标题】如何获得列联表？【英文标题】：How do I get a contingency table? 【发布时间】：2011-11-18 12:34:24 【问题描述】：我正在尝试根据特定类型的数据创建列联表。这对于循环等是可行的......但是因为我的最终表格将包含超过 10E5 个单元格，所以我正在寻找一个预先存在的函数。
我的初始数据如下：
PLANT ANIMAL INTERACTIONS
---------------------- ------------------------------- ------------
Tragopogon_pratensis Propylea_quatuordecimpunctata 1
Anthriscus_sylvestris Rhagonycha_nigriventris 3
Anthriscus_sylvestris Sarcophaga_carnaria 2
Heracleum_sphondylium Sarcophaga_carnaria 1
Anthriscus_sylvestris Sarcophaga_variegata 4
Anthriscus_sylvestris Sphaerophoria_interrupta_Gruppe 3
Cerastium_holosteoides Sphaerophoria_interrupta_Gruppe 1
我想创建一个这样的表：
Propylea_quatuordecimpunctata Rhagonycha_nigriventris Sarcophaga_carnaria Sarcophaga_variegata Sphaerophoria_interrupta_Gruppe
---------------------- ----------------------------- ----------------------- ------------------- -------------------- -------------------------------
Tragopogon_pratensis 1 0 0 0 0
Anthriscus_sylvestris 0 3 2 4 3
Heracleum_sphondylium 0 0 1 0 0
Cerastium_holosteoides 0 0 0 0 1
也就是说，所有植物物种在行中，所有动物物种在列中，有时没有相互作用（而我的初始数据只列出了发生的相互作用）。
【问题讨论】：
列联表中的 10E5 个单元格！！！你在做什么分析？如果您使用卡方检查交互作用，则每个单元格中至少需要有 5 个观察值。 【参考方案1】：与dplyr / tidyr:
df <- read.table(text='PLANT ANIMAL INTERACTIONS
Tragopogon_pratensis Propylea_quatuordecimpunctata 1
Anthriscus_sylvestris Rhagonycha_nigriventris 3
Anthriscus_sylvestris Sarcophaga_carnaria 2
Heracleum_sphondylium Sarcophaga_carnaria 1
Anthriscus_sylvestris Sarcophaga_variegata 4
Anthriscus_sylvestris Sphaerophoria_interrupta_Gruppe 3
Cerastium_holosteoides Sphaerophoria_interrupta_Gruppe 1', header=TRUE)
library(dplyr)
library(tidyr)
df %>% spread(ANIMAL, INTERACTIONS, fill=0)
# PLANT Propylea_quatuordecimpunctata Rhagonycha_nigriventris Sarcophaga_carnaria Sarcophaga_variegata Sphaerophoria_interrupta_Gruppe
# 1 Anthriscus_sylvestris 0 3 2 4 3
# 2 Cerastium_holosteoides 0 0 0 0 1
# 3 Heracleum_sphondylium 0 0 1 0 0
# 4 Tragopogon_pratensis 1 0 0 0 0
【讨论】：
【参考方案2】：只需使用“reshape2”包的dcast()函数：
ans = dcast( df, PLANT~ ANIMAL,value.var = "INTERACTIONS", fill = 0 )
这里“PLANT”将在左列，“ANIMALS”在顶行，表格将使用“INTERACTIONS”填充，“NULL”值将使用 0 填充。
【讨论】：
【参考方案3】：我想指出的是，我们可以在不使用 with 函数的情况下获得与 Andrie 发布的相同结果：
R 基础包
# 3 options
table(warpbreaks[, 2:3])
table(warpbreaks[, c("wool", "tension")])
table(warpbreaks$wool, warpbreaks$tension, dnn = c("wool", "tension"))
tension
wool L M H
A 9 9 9
B 9 9 9
封装 gmodel：
library(gmodels)
# 2 options
CrossTable(warpbreaks$wool, warpbreaks$tension)
CrossTable(warpbreaks$wool, warpbreaks$tension, dnn = c("Wool", "Tension"))
Cell Contents
|-------------------------|
| N |
| Chi-square contribution |
| N / Row Total |
| N / Col Total |
| N / Table Total |
|-------------------------|
Total Observations in Table: 54
| warpbreaks$tension
warpbreaks$wool | L | M | H | Row Total |
----------------|-----------|-----------|-----------|-----------|
A | 9 | 9 | 9 | 27 |
| 0.000 | 0.000 | 0.000 | |
| 0.333 | 0.333 | 0.333 | 0.500 |
| 0.500 | 0.500 | 0.500 | |
| 0.167 | 0.167 | 0.167 | |
----------------|-----------|-----------|-----------|-----------|
B | 9 | 9 | 9 | 27 |
| 0.000 | 0.000 | 0.000 | |
| 0.333 | 0.333 | 0.333 | 0.500 |
| 0.500 | 0.500 | 0.500 | |
| 0.167 | 0.167 | 0.167 | |
----------------|-----------|-----------|-----------|-----------|
Column Total | 18 | 18 | 18 | 54 |
| 0.333 | 0.333 | 0.333 | |
----------------|-----------|-----------|-----------|-----------|
【讨论】：
【参考方案4】：基本 R 中的 xtabs 应该可以工作，例如：
dat <- data.frame(PLANT = c("p1", "p2", "p2", "p4", "p5", "p5", "p6"),
ANIMAL = c("a1", "a2", "a3", "a3", "a4", "a5", "a5"),
INTERACTIONS = c(1,3,2,1,4,3,1),
stringsAsFactors = FALSE)
(x2.table <- xtabs(dat$INTERACTIONS ~ dat$PLANT + dat$ANIMAL))
dat$ANIMAL
dat$PLANT a1 a2 a3 a4 a5
p1 1 0 0 0 0
p2 0 3 2 0 0
p4 0 0 1 0 0
p5 0 0 0 4 3
p6 0 0 0 0 1
chisq.test(x2.table, simulate.p.value = TRUE)
我认为这应该很容易满足您的需求。我不确定它如何在效率方面扩展到 10E5 列联表，但这可能是一个单独的统计问题。
【讨论】：
【参考方案5】：reshape 包应该可以解决问题。
> library(reshape)
> df <- data.frame(PLANT = c("Tragopogon_pratensis","Anthriscus_sylvestris","Anthriscus_sylvestris","Heracleum_sphondylium","Anthriscus_sylvestris","Anthriscus_sylvestris","Cerastium_holosteoides"),
ANIMAL= c("Propylea_quatuordecimpunctata","Rhagonycha_nigriventris","Sarcophaga_carnaria","Sarcophaga_carnaria","Sarcophaga_variegata","Sphaerophoria_interrupta_Gruppe","Sphaerophoria_interrupta_Gruppe"),
INTERACTIONS = c(1,3,2,1,4,3,1),
stringsAsFactors=FALSE)
> df <- melt(df,id.vars=c("PLANT","ANIMAL"))
> df <- cast(df,formula=PLANT~ANIMAL)
> df <- replace(df,is.na(df),0)
> df
PLANT Propylea_quatuordecimpunctata Rhagonycha_nigriventris
1 Anthriscus_sylvestris 0 3
2 Cerastium_holosteoides 0 0
3 Heracleum_sphondylium 0 0
4 Tragopogon_pratensis 1 0
Sarcophaga_carnaria Sarcophaga_variegata Sphaerophoria_interrupta_Gruppe
1 2 4 3
2 0 0 1
3 1 0 0
4 0 0 0
我仍在研究如何解决 order 问题，有什么建议吗？
【讨论】：
你可以用一个命令替换最后三行：cast(PLANT~ANIMAL, data = df, value = "INTERACTIONS", fill = 0) 如果您想根据输入数据帧的排序顺序对结果进行排序，您可以在行上使用order(unique(df$PLANT))，在列上使用明显的模拟。您的示例不需要unique，但每个配对有多个条目且其值相加的版本可能需要它。【参考方案6】：
在基础 R 中，使用 table 或 xtabs：
with(warpbreaks, table(wool, tension))
tension
wool L M H
A 9 9 9
B 9 9 9
xtabs(~wool+tension, data=warpbreaks)
tension
wool L M H
A 9 9 9
B 9 9 9
gmodels 包有一个函数 CrossTable，它提供类似于 SPSS 或 SAS 用户所期望的输出：
library(gmodels)
with(warpbreaks, CrossTable(wool, tension))
Cell Contents
|-------------------------|
| N |
| Chi-square contribution |
| N / Row Total |
| N / Col Total |
| N / Table Total |
|-------------------------|
Total Observations in Table: 54
| tension
wool | L | M | H | Row Total |
-------------|-----------|-----------|-----------|-----------|
A | 9 | 9 | 9 | 27 |
| 0.000 | 0.000 | 0.000 | |
| 0.333 | 0.333 | 0.333 | 0.500 |
| 0.500 | 0.500 | 0.500 | |
| 0.167 | 0.167 | 0.167 | |
-------------|-----------|-----------|-----------|-----------|
B | 9 | 9 | 9 | 27 |
| 0.000 | 0.000 | 0.000 | |
| 0.333 | 0.333 | 0.333 | 0.500 |
| 0.500 | 0.500 | 0.500 | |
| 0.167 | 0.167 | 0.167 | |
-------------|-----------|-----------|-----------|-----------|
Column Total | 18 | 18 | 18 | 54 |
| 0.333 | 0.333 | 0.333 | |
-------------|-----------|-----------|-----------|-----------|
【讨论】：
你能解释一下这些带小数位的数字是什么意思吗？我使用 gmodels 创建了一个列联表，我设置为 TRUE 的唯一参数是 prop.c（即，其他所有参数都设置为 FALSE）。我仍然得到一个额外的数字，它与列百分比和单元格的实际 n 值一起显示......我一生都无法弄清楚它是什么（是的，我已经搜索了很多如何解释输出！）。谢谢。 您的答案在上面的输出中。输出顶部是一个名为Cell Contents 的框。它解释了每个数字的含义，即卡方和各种行和列分数。以上是关于如何获得列联表？的主要内容，如果未能解决你的问题，请参考以下文章