R中卡方的事后检验答案

【问题标题】：Post-Hoc tests for chi-sq in RR中卡方的事后检验
【发布时间】：2017-02-07 17:31:22
【问题描述】：

我有一张像这样的桌子。

> dput(theft_loc)
structure(c(13704L, 14059L, 14263L, 14450L, 14057L, 15503L, 14230L, 
16758L, 15289L, 15499L, 16066L, 15905L, 18531L, 19217L, 12410L, 
13398L, 13308L, 13455L, 13083L, 14111L, 13068L, 19569L, 18771L, 
19626L, 20290L, 19816L, 20923L, 20466L, 20517L, 19377L, 20035L, 
20504L, 20393L, 22409L, 22289L, 7997L, 8106L, 7971L, 8437L, 8246L, 
9090L, 8363L, 7934L, 7874L, 7909L, 8150L, 8191L, 8746L, 8277L, 
27194L, 25220L, 26034L, 27080L, 27334L, 30819L, 30633L, 10452L, 
10848L, 11301L, 11494L, 11265L, 11985L, 11038L, 12104L, 13368L, 
14594L, 14702L, 13891L, 12891L, 12939L), .Dim = c(7L, 10L), .Dimnames = structure(list(
    c("Sunday", "Monday", "Tuesday", "Wednesday", "Thursday", 
    "Friday", "Saturday"), c("BAYVIEW", "CENTRAL", "INGLESIDE", 
    "MISSION", "NORTHERN", "PARK", "RICHMOND", "SOUTHERN", "TARAVAL", 
    "TENDERLOIN")), .Names = c("", "")), class = "table")

我运行了chisq.test，结果显着。我现在想运行一些成对测试，看看意义在哪里。我尝试使用fifer 包和chisq.post.test 函数，但我收到一个错误，上面写着out of workspace。

我还可以通过哪些其他方式进行多重比较测试？

【问题讨论】：

标签： r statistics

【解决方案1】：

这将起作用（在事后测试中尝试chisq.test 而不是默认的fisher.test（精确））：

(Xsq <- chisq.test(theft_loc))  # Prints test summary, p-value very small,
#       Pearson's Chi-squared test
# data:  theft_loc
# X-squared = 1580.1, df = 54, p-value < 2.2e-16 # reject null hypothesis for independence

library(fifer)
chisq.post.hoc(theft_loc, test='chisq.test')

有输出

 Adjusted p-values used the fdr method.

               comparison  raw.p  adj.p
1       Sunday vs. Monday 0.0000 0.0000
2      Sunday vs. Tuesday 0.0000 0.0000
3    Sunday vs. Wednesday 0.0000 0.0000
4     Sunday vs. Thursday 0.0000 0.0000
5       Sunday vs. Friday 0.0000 0.0000
6     Sunday vs. Saturday 0.0000 0.0000
7      Monday vs. Tuesday 0.0000 0.0000
8    Monday vs. Wednesday 0.0000 0.0000
9     Monday vs. Thursday 0.0000 0.0000
10      Monday vs. Friday 0.0000 0.0000
11    Monday vs. Saturday 0.0000 0.0000
12  Tuesday vs. Wednesday 0.1451 0.1451
13   Tuesday vs. Thursday 0.0000 0.0000
14     Tuesday vs. Friday 0.0000 0.0000
15   Tuesday vs. Saturday 0.0000 0.0000
16 Wednesday vs. Thursday 0.0016 0.0017
17   Wednesday vs. Friday 0.0000 0.0000
18 Wednesday vs. Saturday 0.0000 0.0000
19    Thursday vs. Friday 0.0000 0.0000
20  Thursday vs. Saturday 0.0000 0.0000
21    Friday vs. Saturday 0.0000 0.0000

正如我们所见，除了一对之外的所有成对测试都很重要，我们也可以使用不同的p-value-correction（通过将control从默认的fdr更改为bonferroni）。

【讨论】：

完美！！难怪它为什么不起作用。有没有办法比较一周中的几天和 pdDistrict？我的主要问题是试图了解数据中的差异所在。有没有犯罪率较高的一天和地点。
嗯，通过事后测试，我认为我们可以得出哪些对是显着的，为了衡量哪一天/哪个位置的犯罪率更高，我们需要一个不同的测试。
我明白了......所以也许 1) aov(n ~ place+day, data = df)->fit , 2) TukeyHSD(fit) 使用与上述评论相同的 df。跨度>
刚刚实现了一个未被充分利用的 ggplot facet 技巧。您可以像这样查看一天/地点组合的罪行：ggplot(df,aes(x=place,y=n))+geom_bar(stat="identity")+facet_grid(.~day+place,scale="免费”)+主题(axis.text.x=element_text(angle=90))+主题(strip.text.x=element_text(angle=90))。南方的犯罪率一直很高。
使用泊松分布的 GLM 来确定较高频率出现的位置是否更有意义？ glm(freq~PdDistrict+DayOfWeek,df,family ='poisson')

【解决方案2】：

由于此处不再维护fifer，因此使用RVAideMemoire 的解决方案（此处更详细地描述https://rdrr.io/cran/RVAideMemoire/src/R/chisq.multcomp.R）：

install.packages("RVAideMemoire")
library(RVAideMemoire)
chisq.multcomp(theft_loc, p.method = "none")
>      7874    7909    7934    7971    7997    8106    8150    8191    8246    8277    8363    8437    8746    9090    10452   10848   11038   11265   11301  
7909  0.78056 -       -       -       -       -       -       -       -       -       -       -       -       -       -       -       -       -       -      
7934  0.63321 0.84256 -       -       -       -       -       -       -       -       -       -       -       -       -       -       -       -       -      
7971  0.44095 0.62272 0.76923 -       -       -       -       -       -       -       -       -       -       -       -       -       -       -       -      
7997  0.32889 0.48533 0.61768 0.83698 -       -       -       -       -       -       -       -       -       -       -       -       -       -       -      
8106  0.06647 0.11954 0.17444 0.28701 0.39036 -       -       -       -       -       -       -       -       -       -       -       -       -       -      
8150  0.02923 0.05720 0.08854 0.15860 0.22857 0.73002 -       -       -       -       -       -       -       -       -       -       -       -       -      
8191  0.01238 0.02625 0.04298 0.08354 0.12732 0.50552 0.74841 -       -       -       -       -       -       -       -       -       -       -       -      
8246  0.00339 0.00802 0.01417 0.03081 0.05073 0.27360 0.45342 0.66793 -       -       -       -       -       -       -       -       -       -       -      
8277  0.00152 0.00382 0.00706 0.01637 0.02817 0.18156 0.32174 0.50276 0.80943 -       -       -       -       -       -       -       -       -       -      
8363  0.00012 0.00037 0.00078 0.00216 0.00422 0.04522 0.09741 0.18128 0.36396 0.50497 -       -       -       -       -       -       -       -       -      
8437  1.0e-05 3.6e-05 8.5e-05 0.00027 0.00060 0.01007 0.02585 0.05643 0.13921 0.21586 0.56805 -       -       -       -       -       -       -       -      
8746  1.3e-11 8.8e-11 3.2e-10 2.0e-09 7.1e-09 8.2e-07 4.5e-06 2.0e-05 0.00013 0.00032 0.00341 0.01841 -       -       -       -       -       -       -      
9090  < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 6.2e-14 8.1e-13 8.0e-12 1.5e-10 6.9e-10 3.7e-08 8.1e-07 0.01000 -       -       -       -       -       -      
10452 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 -       -       -       -       -

类别替换为每个类别的计数。我不喜欢对多重比较进行更正（请参阅下面的参考文献进行讨论），但fdr 是可用的。

莫兰，医学博士（2003 年）。在生态学研究中拒绝顺序 Bonferroni 的论据。奥伊克斯。

【讨论】：

再次 RVAideMemoire 依赖于 mixOmics，这在我使用的 R 版本 3.6.2 中不可用。你必须使用安装包("package_url",repo=NULL,type="source")
谢谢，不知道为什么这么基本的东西在别处没有实现，可能写起来很鸡肋。