R语言系列第四期：③R语言表格数据率的比较
本文首发于 ”百味科研芝士“ 微信公众号，转载请注明：百味科研芝士，Focus科研人的百味需求
连续型数据的组间比较往往可以采用t检验/wilcoxon检验或者ANOVA方差分析/KW检验来完成。但是对于分类资料来说，这些方法就是行不通的了。详情点击：
R语言系列第四期：①R语言单样本双样本差异性检验
R语言系列第四期：②R语言多组样本方差分析与KW检验
在这个部分我们会介绍一系列用于分析表格数据的函数，我们会着重看prop.test()，binom.test() , chisq.test() 以及 fisher.test()函数。
A▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼与前面连续型资料分析的顺序一样，我们先来介绍一下单样本率的检验
▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲我们考虑这样一个例子（Altman，1991），这里215名病人中39名被观测到患有哮喘，然后有人对“随机病人”患有哮喘的概率是0.15这个假设做检验。我们可以用函数prop.test()来做检验：
> prop.test(39,215,0.15)
1-sample proportions test with continuity correction
data: 39 out of 215, null probability 0.15
X-squared = 1.425, df = 1, p-value = 0.2326
alternative hypothesis: true p is not equal to 0.15
95 percent confidence interval:
0.1335937 0.2408799
sample estimates:
p
0.1813953
#Tips：函数prop.test()中的三个参数分别是阳性观测数，总数，以及参考概率值。最后一个参数如果没有设定的话，那么默认是0.5。我们可以从当前的结果中得到的结论是无差别（P>0.05），而且在95%置信区间里我们可以看出范围包含了0.15，所以也可认为无差别。
这里的0.15是人为构造出来的。但是如果我们有一组这样的数据，往往更希望得到这个概率参数的置信区间，这里输出结果的结尾已经给我们算好了。
除此之外，还可以利用函数binom.test()在二项分布下做检验。这是你能得到精确的检验概率，所以一般比上一个prop.test()更被人喜欢。不过prop.test()除了单比例检验之外还能做其他的事情。为了得到这个p值，我们先计算出x取每一个可能的点的概率，然后再将观测到的小于等于x的概率都加在一起。
> binom.test(39,215,0.15)
Exact binomial test
data: 39 and 215
number of successes = 39, number of trials = 215, p-value = 0.2135
alternative hypothesis: true probability of success is not equal to 0.15
95 percent confidence interval:
0.1322842 0.2395223
sample estimates:
probability of success
0.1813953
#Tips：在检验水准为0.05水平下的置信区间其实是从0.025水平下的双边检验中得到的。就是说两变分别是0.025和0.975的水平。这两个计算方式的结果是一致的。
B▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼两个独立的比例比较
▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲前面提到的函数prop.test()也能够用于比较两个或多个比例。这种情况下它的参数就应该是两个向量，前一个向量里存放各组阳性观测数，后一个向量里是每组总数。
我们还是拿Altman的例子来演示（可以理解成两位医生手术成功率比较）：
> success<-c(9,4)
> total<-c(12,13)
> prop.test(success,total)
2-sample test for equality of proportions with continuity correction
data: success out of total
X-squared = 3.2793, df = 1, p-value = 0.07016
alternative hypothesis: two.sided
95 percent confidence interval:
0.01151032 0.87310506
sample estimates:
prop 1 prop 2
0.7500000 0.3076923
#Tips：可以看出来R会将数据自动识别成两样本率的差异检验，求的是阳性结果的概率差异。然后我们得到的P值是0.07016，这个值与0.05很近，得到的是无差异的阴性结果。但是95%的置信区间为[0.011531032,0.87310506]这个范围没有包含0，这个置信区间是比例之差的置信区间，它的结论是不可以认为两个医生的手术成功率是一样的阳性结果，二者的差异是由置信区间和假设检验使用的是不同的近似方法导致的。
我们还可以用Fisher精确概率法检验。这个检验在给定行和列的边际值的情况下计算2*2表格的条件分布。简单来讲Fisher确切概率法的结果是在这个四格表里出现当前数据情况或者更有利于阳性结果的概率之和，它会计算每一种情况的概率，然后再累加起来。
相关的函数就是fisher.test()，他要求输入的数据是矩阵形式的，如下：
> a<-matrix(c(9,4,3,9),nrow=2)
> a
[,1] [,2]
[1,] 9 3
[2,] 4 9
> fisher.test(a)
Fisher's Exact Test for Count Data
data: a
p-value = 0.04718
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.9006803 57.2549701
sample estimates:
odds ratio
6.180528
#Tips：这里的置信区间是比值比的结果，也就是计算(p1/(1-p1))/(p2/(1-p2))的区间，是一个衡量Fisher检验中相关程度的指标，得到的结果可以跟1比较。不过这里的结果同样和假设检验的结果相矛盾，原因同上。
和fisher.test()一样，在chisq.test()中的标准χ2检验需要矩阵类型的数据源。而作为一个2*2表格来说，这个检验与prop.test()的结果是完全一致的。
> chisq.test(a)
Pearson's Chi-squared test with Yates' continuity correction
data: a
X-squared = 4.3672, df = 1, p-value = 0.03664
C▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼趋势检验，多组率
▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲有的时候你想要比较多于两个部分，而这些部分又是有序的分类，所以你想找到一个随着分组序号递增或者递减的趋势。
这个部分我们使用Altman的数据，这个例子记录了一组女性是否使用剖腹产生育孩子，以及对应产妇鞋子码数的数据，数据在R语言ISwR数据包里。
> library(ISwR)
> caesar.shoe
<4 4 4.5 5 5.5 6+
Yes 5 7 6 7 8 10
No 17 28 36 41 46 140
#Tips：这里呢，有一个基于正态近似的检验可以使用。这个检验计算每组的观测比例和所有组的比例之间的加权平方和的偏差。检验统计量近似服从自由度为k-1的χ2分布。
为了使用prop.test()，我们需要将数据转化成两个分别放有阳性数据和总数的变量里：
> caesar.shoe.yes<-caesar.shoe['Yes',]
> caesar.shoe.total<-margin.table(caesar.shoe,2)
> caesar.shoe.yes
<4 4 4.5 5 5.5 6+
5 7 6 7 8 10
> caesar.shoe.total
<4 4 4.5 5 5.5 6+
22 35 42 48 54 150
之后我们就可以检验了：
> prop.test(caesar.shoe.yes,caesar.shoe.total)
6-sample test for equality of proportions without continuity correction
data: caesar.shoe.yes out of caesar.shoe.total
X-squared = 9.2874, df = 5, p-value = 0.09814
alternative hypothesis: two.sided
sample estimates:
prop 1 prop 2 prop 3 prop 4 prop 5 prop 6
0.22727273 0.20000000 0.14285714 0.14583333 0.14814815 0.06666667
Warning message:
In prop.test(caesar.shoe.yes, caesar.shoe.total) :
Chi-squared近似算法有可能不准
#Tips：得到结果是大于0.05的，结果并不显著，但是从细节上来说，剖腹产这一组的数值都相对较小，而在R语言的警告里，给我们提示，有一些格子的理论频数小于5，算法结果可能不准确。
同时可以使用函数prop.tend.test()函数来检测不同部分的趋势。它有3个参数：x,n和score。前两个与prop.test()中一致，而最后一个是赋予每组的分数，默认是简单的1,2，...，k，这些数据使得我们前后每组之间是有顺序的。这个检验的本质是一个用分数对不同部分进行的加权线性回归，我们对当前的数据进行检验，就成为了一个自由度为1的χ2检验。
> prop.trend.test(caesar.shoe.yes,caesar.shoe.total)
Chi-squared Test for Trend in Proportions
data: caesar.shoe.yes out of caesar.shoe.total ,
using scores: 1 2 3 4 5 6
X-squared = 8.0237, df = 1, p-value = 0.004617
所以，如果我们假设鞋子码数是线性的，那么我们可以看到一个显著的趋势。
D▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼r × c表格
▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲为了分析行列都多于两个分类的表格数据，可以使用函数chisq.test()和 fisher.test()，但是后者在每一格数字比较大而且超出两行或两列时的计算量非常大。我们使用我们之前在分类数据描述的章节中采用的例子，那个婚姻状况与咖啡因消费情况的数据：
> caff.marital<-matrix(c(652,1537,598,242,36,46,38,21,218,327,106,67),
+ nrow=3,byrow=T)
> colnames(caff.marital)<-c("0","1-150","151-300",">300")
> rownames(caff.marital)<-c("Married","Divorced","Single")
> caff.marital
0 1-150 151-300 >300
Married 652 1537 598 242
Divorced 36 46 38 21
Single 218 327 106 67
> chisq.test(caff.marital)
Pearson's Chi-squared test
data: caff.marital
X-squared = 51.656, df = 6, p-value = 2.187e-09
检验结果显示高度显著。当然，我们也可以查看chisq.test()函数的一些额外的返回值。比如每个格子的期望值：
> chisq.test(caff.marital)$expected
0 1-150 151-300 >300
Married 705.83179 1488.01183 578.06533 257.09105
Divorced 32.85648 69.26698 26.90895 11.96759
Single 167.31173 352.72119 137.02572 60.94136
然后如果你有兴趣可以自己手动计算卡方值统计量，按照χ2检验的计算公式，我们可以用（A-T）2/T来表示（A实际频数，T理论频数）：
> chisq.test(caff.marital)$expected->T
> chisq.test(caff.marital)$observed->A
> (A-T)^2/T
0 1-150 151-300 >300
Married 4.1055981 1.612783 0.6874502 0.8858331
Divorced 0.3007537 7.815444 4.5713926 6.8171090
Single 15.3563704 1.875645 7.0249243 0.6023355
#Tips：我们把每个格子里的数值加起来就是总的统计量大小，而每个格子里的数值代表的是各自的贡献量。
也可以对原始数据使用chisq.test()，这里我们使用之前的juul数据作为例子：
> attach(juul)
> chisq.test(tanner,sex)
Pearson's Chi-squared test
data: tanner and sex
X-squared = 28.867, df = 4, p-value = 8.318e-06
> chisq.test(tanner,sex)$observed
sex
tanner 1 2
1 291 224
2 55 48
3 34 38
4 41 40
5 124 204
#Tips：数据会被汇总成为一个表格里。然后通过χ2检验来计算差异。
关于表格数据的统计分析就介绍到这里了，我们下期再见。
参考资料：
1.《R语言统计入门（第二版）》 人民邮电出版社 Peter Dalgaard著
2.《R语言初学者指南》 人民邮电出版社 Brian Dennis著
3.Vicky的小笔记本《blooming for you》 by Vicky
