我有一个大约有20000个观测值的数据框架。由此,我创建了一个频率为两个变量的列联表。
有了这个,我想对独立性进行卡方检验,看看我的两个变量之间是否存在关系。通常情况下,这很容易,但尽管原始数据帧的大小很大,但许多单元格的预期值为0。我想删除任何包含频率小于5的行。
我已经广泛搜索了堆栈交换,但我找不到解决这个特定问题的方法,我要么a)理解(我对R相对陌生),要么b)使用列联表而不是原始数据帧。
非常感谢您的帮助。
编辑:
谢谢你的回复贾斯汀。
根据请求,我上传了数据帧和列联表的摘录。我还上传了我迄今为止尝试过的少量代码,并取得了结果。
数据帧
Department Super
AAP 1
ACS 4
ACE 1
AMA 1
APS 3
APS 2
APS 1
APS 1
ARC 5
ARC 7
ARC 1
BIB 6
BIB 6
BMS 2
所以有两列,第一列是三个字母的部门代码,第二列是一位数的整数(1-7)。
应急表
table(department,super)
1 2 3 4 5 6 7 8
ACS 32 10 7 24 50 7 24 14
AMA 0 4 2 6 10 3 11 1
...
所以一个标准的频率列联表。
到目前为止,我知道我可以创建一个逻辑测试,测试单元格内容是否小于5:
depSupCrosstab <- depSupCrosstab[,2:8]>5
我不知道的是,如果有任何FALSE条目,如何使用这行代码创建的矩阵来删除整行。
希望能有所帮助。恐怕我是新手,但只有一种方法可以学习。。。
恐怕您的问题更复杂。卡方检验的假设是每个细胞的预期频率大于5。在您的示例中,您试图选择列联表中每个单元格的计数,这是观察到的频率。预期频率(在零假设下)是根据行和列总数计算的,如这里的基本示例所示。
举个例子,一个假设的列联表可能看起来像:
ACS <- c(32, 10, 7, 24, 50, 7, 24, 14)
AMA <- c(0, 4, 2, 6, 10, 3, 11, 1)
ARC <- c(6, 10, 12, 3, 12, 23, 10, 2)
tab <- rbind(ACS, AMA, ARC)
如果您筛选观察到的计数等于或小于5,您将删除AMA和ARC:
apply(tab,1, function(x) any(x<=5))
ACS AMA ARC
FALSE TRUE TRUE
这在概念上是错误的,因为如上所述,预期频率取决于整个数据。获取exp.计数:
chisq.test(tab, correct=F)$expected
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
ACS 22.558304 14.247350 12.466431 19.590106 42.742049 19.590106 26.713781
AMA 4.968198 3.137809 2.745583 4.314488 9.413428 4.314488 5.883392
ARC 10.473498 6.614841 5.787986 9.095406 19.844523 9.095406 12.402827
[,8]
ACS 10.091873
AMA 2.222615
ARC 4.685512
Warning message:
In chisq.test(tab, correct = F): Chi-squared approximation may be incorrect
卡方测试会发出警告消息,因为确实有一些单元格的exp.计数小于5。但是,如果只删除AMA,则表的动态(行和列总数)会发生变化,并且所有exp.计数都在5:以上
chisq.test(tab[-2,], correct=F)$expected
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
ACS 25.95122 13.658537 12.97561 18.439024 42.34146 20.487805 23.21951
ARC 12.04878 6.341463 6.02439 8.560976 19.65854 9.512195 10.78049
[,8]
ACS 10.926829
ARC 5.073171
因此,如果你同时删除AMA和ARC,你就会丢失一个重要信息。
您可以尝试运行Fisher精确测试(请参阅下面的解释):
fisher.test(tab,simulate.p.value=TRUE,B=10000)
总结:
- 单个观察到的频率是预期频率的不良指标。观察到的频率可能低于5,但该小区的预期频率将高于5
- 在大型列联表中,可以接受高达20%的经验频率低于5,但结果是统计能力的损失,因此测试可能无法检测到真正的效果。即使在这种情况下,exp.频率也不应该低于1
- 低经验频率的另一种测试方法是Fisher精确测试。卡方检验统计量近似卡方分布。如果样本量大,则该近似变得更准确,因此需要exp.frequency>5.即使样本量很小,Fisher精确检验也会计算卡方统计量的精确概率,但它可能更需要计算。不幸的是,对于大于2x2的列联表,您可能需要模拟p值,这有其自身的局限性(此处没有讨论空间,但这是一个很好的研究主题)。选择大量复制进行模拟(B),并对其进行调整,以查看解决方案的健壮性
我想我在一个相关的问题中找到了答案。在这种情况下,apply
是您的朋友,因为它可以迭代列或行。
要创建与您类似的数据帧,然后只选择所有列>5的行,可以使用以下方法:
set.seed(1985)
tosub <- data.frame(matrix(round(runif(n = 80, min = 0, max = 100)), ncol = 8))
head(tosub,2)
x <- apply(tosub[,1:8] > 5, MARGIN = 1, all)
summary(x)
tosub[which(x),]
X1 X2 X3 X4 X5 X6 X7 X8
1 66 30 72 59 26 69 76 47
2 27 42 26 95 66 14 67 18
4 42 28 93 7 35 35 95 23
5 38 89 69 91 98 91 60 69
9 89 31 91 72 28 31 58 58
10 53 87 27 89 95 37 98 20