卡方分布与卡方检验笔记

发表于 | 分类于
| 字数统计: 9,764 | 阅读时长 ≈ 38

什么是卡方分析

第1个案例

先看一个案例，这个案例会告诉我们什么是卡方分析，这个案例源于《深入浅出统计学》。

某个赌场遇到了一个问题，也就是赌场的老虎机总是出头奖，轮盘总是停在12位，这样赌客就很容易赢钱，如果再这么下赌下去，赌场撑不住，赌场的老板怀疑有人在老虎机上动了手脚，此时就来研究，老虎机是否被动了手脚。下面是某台老虎机的期望概率分布，其中X代表每一局游戏的净收益，如下所示：

x -2 23 48 73 98
P(X=x) 0.977 0.008 0.008 0.006 0.001

其中-2表示，每局的赌注是2美元，如果你没赢的话，你就损失2美元。现在赌场收集了一些统计数据，给出了人们获得某种收益的次数，如下所示：

x -2 23 48 73 98
频数 965 10 9 9 7

现在我们要将每个x值的实际频数与根据概率分布得出的期望频数进行比较，此时时我们定义观察的频数为1000，现在用1000分别乘以概率分布，得到期望频数，如下所示：

x 观察频数 期望频数
-2 965 977
23 10 8
48 9 8
73 9 6
98 7 1

从这个表中我们可以知道，观察频数与期望频数有差别，但是这种差异有多，我们并不清楚。因此我们需要进行某种假设检验，以检验观察频数和期望频数之间的差别，这样我们就能够判定：老虎机被人动过手脚。现在的问题就是，我们使用哪种分布进行这项假设检验？

答案就是使用卡方检验来评估差异，卡方分布通过一个检验统计量来比较期望结果和实际结果之间的差别，然后得出观察频数极值的发生概率。

现在我们先来求一下检验统计量，为此，首先画出一张表，填入相应问题的观察频数和期望频数，然后，用观察频数和期望频数计算下列统计量，其中O代表观察频数，E代表期望频数，如下所示：

mark

也就是说，对于概率分布中的每一个概率，取期望频数和实际频数的差，求差的平方数，再除以期望频数，最后将结果相加，现在我们求一下老虎机问题的检验统计量是多少。

计算过程如下所示：

mark

如果$\chi^2$的值很小，说明观察频数和期望频数之间的差别不显著，卡方越大，差别越显著。

再回到卡方分析的统计量上，检验统计量$\chi^2$提供了一种对观察频数和期望频数之间的差异进行量度的办法，$\chi^2$的数值越小，观察频数和期望频数之间的总差值越小。除数E为期望频数，于是所得的结果与期望频数成反比，如下所示：

mark

$\chi^2$大到什么程度才算得上显著呢？此处先不做解答，现在先了解一下卡方分布。

卡方分布简介

卡方分布是由正态分布推导出来的分布，它的定义为，n个独立标准正态变量的平方和称为有n个自由度的$\chi^2$分布，记为$\chi^2(n)$，$\chi^2(n)$的总体均值为n，总体方差为2n（不知道这是如何推导出来的）。

此外，若干个独立的卡方分布变量的和也能构成卡方方分布，其自由度等于那些卡方分布自由度之和，卡方分布也是一族分布，由该族成员的不同自由度来区分，由于上方分布是正态变量的平方和，因此，它不会取负值（《从数据到结论》吴喜之 ）。

现在再看一下卡方分析的简介。

卡方分析简介

两个样本率的比较通常是使用$\chi^2$检验（Chi-square test）。$\chi^2$检验是以$\chi^2$分布为理论依据，用途广泛的假设检验方法，卡方分析用于比较不同组之间的构成比，它的零假设是假定各组之间的构成是相同的，计算出理论每组的理论构成比，再计算理论值与实际值的差别，如果差别大的话，就拒绝零假设。它用于推断两个总体率或构成比之间有无差别、多个总体率或构成比之间有无差别、多个样品率间的多重比较、两个分类变量之间有无关联性、多维列联表的分析和频数分页拟合优度的$\chi^2$检验，它的扩展分析方法有Fisher精确分析，Ridit分析，CMH分析。

卡方分布用途

卡方分布的两个主要用途。卡方概率分布主要用于检查实际结果与期望结果之间何时存在显著差别，该概率分布使用前前面讲到的检验统计量$\chi^2$进行检验。卡方分布主要有两个用途。

第一用于检验拟合优度，也就是可以检验一组给定的数据与指定分布的吻合程度。例如，可以用它检验老虎机收益的观察频率与我们所期望的分布的吻合程度。$\chi^2$分布的另一个用途是检验两个变量的独立性，通过这个方法可以检查变量之间是否存在某种关联。$\chi^2$分布用到一个参数，即希腊字母$\nu$，读作“纽”，下面是可以看一下$\nu$如何影响概率分布的形状。

卡方分布的概率密度曲线

卡方分布（Chi-square distribution）是一种连续型分布，$\chi^2$分布只有一个参数，即自由度$\nu$。按$\chi^2$分布的密度函数$f(\chi^2)$可以给出自由度$\nu=1,2,3,\dots$的一簇$\chi^2$分布曲线，由$\chi^2$分布曲线可见，$\chi^2$分布的形状依赖于自由度$\nu$的大小：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
x <- c(seq(0,15,length = 1000)) 
#生成从0到10之间的1000个数字
y1 <- dchisq(x,1) 
# 生成对x的卡方密度值，自由度为1
y2 <- dchisq(x,2) 
# 生成对x的卡方密度值，自由度为1
y6 <- dchisq(x,6) 
# 生成对x的卡方密度值，自由度为6
y10 <- dchisq(x,10) 
# 生成对x的卡方密度值，自由度为10
plot(x,y1,type = "l", ylim = c(0,0.5),xlim=c(0,15),xlab="", ylab="", lty = 1, main = expression(paste(chi^2,"Distribution")))
# type ="l"表示绘制线，如果是n则不绘线；lty=1表示绘制实线
lines(x,y2, type = "l", xlab = "", ylab = "", lty =2) 
lines(x,y6, type = "l", xlab = "", ylab = "", lty =2) 
# lty=2表示绘制虚线lines(x,y10, type = "l", xlab = "", ylab = "", lty= 3) # lty=3表示绘制点线
lines(x,y10, type = "l", xlab = "", ylab = "", lty= 3) # lty=3表示绘制点线
text(0.2,0.3, c(expression(chi^2*(1))))
text(1.9,0.3, c(expression(chi^2*(2))))
text(4,0.15, c(expression(chi^2*(6))))
text(8,0.11, c(expression(chi^2*(10))))

图片如下所示：

mark

  1. 当自由度$\nu \leq 2$时，曲线呈L形，也就是一条先高后低的平滑曲线，检验统计量等于较小的数值的概率远远高于等于较大数值的概率，这就是说，观察频数有可能接近期望频数；
  2. 当$\nu>2$时，卡方分布的曲线就会表现出先低后高，再低的特征，其外形沿着正向扭曲，但是当，随着$\nu$的增加，曲线逐渐趋于对称；
  3. 当自由度$\nu \rightarrow \infty$时，$\chi^2$分布趋于正态分布。

卡方分布的一个基本性质是它的可加性；如果两个独立的随机变量$X_{1}$和$X_{2}$分别服从自由度$\nu_{1}$和$\nu_{2}$的$\chi^2$分布，则它们的和$X_{1}+X_{2}$服从自由度$\nu_{1}$和$nu_{2}$的$\chi^2$分布，即$X_{1}+X_{2}~\chi^2_{\nu1+\nu2}$.

如果使用具有特定参数$\nu$的$\chi^2$分布以及检验统计量${X}^2$，可以简单记为如下形式：

mark

$\chi^2$分布的界值：当自由度$\nu$确定后，$\chi^2$分布曲线下右侧尾部的面积为$\alpha$时，横轴上相应的$\chi^2$值记为$\chi^2_{\alpha,\nu}$，即$\chi^2$分布的界值。$\chi^2$值越大，p值越小，反之，$\chi^2$值越小，界值越大。

卡方分析的显著性

现在回答前面的那个问题，即检验统计量$\chi^2$大到什么程度才算得上显著呢？

这要取决于显著性水平，用卡方分布进行的检验为单尾检验，右尾被称作拒绝域，于是，通过查看检验统计量是否位于右尾的拒绝域内，就可以判定根据期望分布得出目前结果的可能性，如果用显著性水平$\alpha$进行检验，就写为$\chi^2_{\alpha}(\nu)$，在以前没有计算机的时候，通常是使用卡方概率表来计算卡方分布的拒绝域，下面先讲一下如何使用卡方概率表，理解了这个，就方便理解软件的计算结果。

为了求出临界值，首选找出自由度$\nu$以及显著性水平$\alpha$，在第一列查找$\nu$，第一行查找$\alpha$，交点即x值，从$P(\chi^2\alpha(\nu) \geq(x)=\alpha$中得出临界值。

以5%为显著性水平，8为自由度进行检验，若要得到临界值，那么在第一列查找8，第一行查找0.05，查出数值为15.51，因此只要统计量$\chi^2$大于15.51，则在显著性水平为5%，自由度为8的情况下，检验统计量就位于拒绝域以内，如下所示：

mark

我们在前面通过卡方分析计算过老虎机的理论期望频数与实际频数的卡方值是38.272，它的值大于15.51，也就是位于拒绝域内，因为我们可以判定，老虎机做了手脚，事实上，我们还在前面使用了R来计算这个卡方值，结论也是，p值远小于0.05，如下所示：

1
2
3
p <- c(0.977,0.008,0.008,0.006,0.001)  # 期望概率
x <- c(965,10,9,9,7) # 实际的频数
chisq.test(x,p=p)

计算结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
> p <- c(0.977,0.008,0.008,0.006,0.001) 
> x <- c(965,10,9,9,7) 
> chisq.test(x,p=p)
        Chi-squared test for given probabilities
data:  x
X-squared = 38.272, df = 4, p-value = 9.845e-08
Warning message:
In chisq.test(x, p = p) : Chi-squared近似算法有可能不准

结果显示，卡方值为38.272，自由度为4，p值远小于0.01。

卡方分析的假设检验步骤

利用卡方分布进行假设检验的步骤与其他方法进行假设检验的步骤基本上一样，包括哪下步好：

  1. 确定进行检验的假设及其备选假设；
  2. 求出期望频数和自由度；
  3. 确定用于决策的拒绝域；
  4. 计算检验统计量$\chi^2$；
  5. 查看检验统计量是否位于拒绝域以内；
  6. 作出决策。

第2个案例

再看一个案例，这个案例还是来源于《深入浅出统计学》，案例是这个样子的：

赌场中的某个赌桌上的骰子有问题，赌场老板列出了这个股子的观察频数，如下所示：

数值 1 2 3 4 5 6
频数 107 198 192 125 132 248

现在我们要计算一下，这个股子到底有没有问题，我们按照卡方检验的步骤来进行。

第一步：决定要进行检验的假设和备选假设。

H0（假设）：骰子没问题，每一面数值出现的概率都相同，也就是1/6。

H1（备选假设）：骰子有问题。

第二步：求出期望频数和自由度。

根据上面的数据，我们知道，所有的观察频数相加为1002，现在我们计算一下理论上每个面出现的频数，如下所示：

x 观察频数 期望频数
1 107 167
2 198 167
3 192 167
4 125 167
5 132 167
6 248 167

自由度为5，这个很好理解，6个期望频数，总和为1002，即我们必须求出6个信息，同时受到1个限制，因此$\nu$=6-1=5.

第三：确定用于做决策的拒绝域。

从概率表中可以查出$\chi^2_{1\%}(5)=15.09$，那么拒绝域为$\chi^2>15.09$的范围。

第四：计算检验统计量$\chi^2$，如下所示：

mark

第五：查提检验统计量是否位于拒绝域以内

拒绝域由$\chi^2>15.09$决定，由于我们计算出来的$\chi^2=88.24$，因此检验统计量位于拉拒绝域内。

第六：做出决策。

由于我们计算出来的检验统计量位于拒绝域内，因此说明在显著性水平为1%的情况下，有足够的证据拒绝原假设，于是我们接受备选假设，骰子有问题。

现在用R来过计算一下，如下所示：

1
2
3
4
5
6
7
8
> p <- c(1/6,1/6,1/6,1/6,1/6,1/6) 
> x <- c(107,198,192,125,132,248) 
> chisq.test(x,p=p)
        Chi-squared test for given probabilities
data:  x
X-squared = 88.24, df = 5, p-value < 2.2e-16

在这个案例中，我们可以发现，卡方分析不仅能够比较两个率的大小，还能计算出理论频数与实际频数的差异，实际上，这就是拟合优度（Goodness of Fit）的思想，即看实际值与理论值（预测值）的差异有多大，在不同的场合下，拟合优度有不同的用处，例如，在组间比较时，通过比较实际值与理论值的差异大小，做出是否有统计学差异的结论；在模型评价中，根据实际值与模型值的差异大小，判断模型是否很好地拟合了数据；在判断数据是否服从某一分布时，根据其理论分布与实际数据的差异做出判断，等等。

以下是其他的一些案例练习。

四格表资料的$\chi^2$检验

案例1

某医院欲比较异梨醇口服液（试验组）和氢氯噻嗪+地塞米松（对照组）降低颅内压的疗效。将200例颅内压增高症患者随机分为两组，结果见表7-1。问两组降低颅内压的总体有效率有无差别？（《医学统计学》第四版，孙振球，第七章 卡方检验 例7-1，P96）


1
2
data71<-matrix(c(99,75,5,21),nr = 2,dimnames = list(c("试验组", "对照组"),c("有效", "无效")))
data71

数据的四格表如下所示：

这张图里的四个数是该表的基本数据，其余的数据都是由这4个基本数据推算出来的，称为四格表（fourfold table）资料。

1
2
chi <- chisq.test(data71)
chi$expected #　显示理论数

理论数如下所示：

1
2
chisq.test(data71,correct=F) 
# correct = F表示不进行Yates连续性校正

结果：卡方值为12.8571，p值= 0.0003362小于0.01，按sig = 0.01的水准，拒绝原假设，接受备用假设，可以认为异山梨醇口服液降低颅内压的有效率高于氢氯噻嗪+地塞米松的有效率。这个结果与SAS计算的结果相同。R语言中默认的卡方检验是进行Yates连续性校正，因此上题中进行检验时，要加入correct=F。

四格表资料的Fisher确切概率法

四格表中的计数资料的实际频数都是分类资料，是不连续的，因此计算出的卡方值都是离散型分布，而卡方界值表的依据是卡方分布，卡方分布是连续分布，因此在实际工作中，对于四格表资料，需要进行Yates校正，其规则如下：

  1. 当$n\geq 40$，且所有的$T\geq 5$时，用卡方检验的基本公式。
  2. 当$n\geq 40$，但有$1\leq T \geq 5$时，用校正公式，或改用四格表资料的Fisher确切概率法。
  3. 当$n\leq 40$，或$T \leq 1$时，用四格表资料的Fisher确切概率法。

案例2

某医师欲比较胞磷胆碱与神经节苷酯治疗脑血管疾病的疗效，将78例脑血管疾病患者随机分为两组，结果见表7-2。问两种药物治疗脑血管疾病的有效率是否有差别？（《医学统计学》第四版，孙振球，第七章 卡方检验，例7-2）

1
2
3
4
5
6
7
8
9
10
data72<-matrix(c(46,18,6,8),
               nr = 2,
               dimnames = list(c("胞磷胆碱","神经节苷脂"),c("有效","无效")))
chisq.test(data72,correct=F) # 不进行Yates校正的结果
chisq.test(data72,correct=T) # 进行Yates校正的结果
# 这个结果与下面的结果一样
# 因为R中进行卡方检验时是默认进行Yates连续性校正的
chisq.test(data72)

结果如下所示：

不校正的情况下，卡方值为4.35，p值为0.03695小于0.05，进行卡方连续性校正的情况下，卡方值为3.1448,p=0.07617。校正与不校正得出的结果在0.05附近，因此下结论要谨慎。

配对四格表卡方检验

案例3

某实验室分别用乳胶凝集法和免疫荧光法对58名可疑系统红斑狼疮患者血清中抗核抗体进行测定，结果见表7-3。问两种方法的检测结果有无差别？（《医学统计学》第三版，孙振球，第七章 卡方检验，例7-3）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
data73<-matrix(c(11,2,12,33),
               nr = 2,dimnames = list(c("+", "-"),c("+", "-")))
library(epiR) # 载入epiR包，使用其中的epi.kappa函数
mcnemar.test(data73) # 配对检验
chisq.test(data73)$expect
epi.kappa(data73)$kappa 
# 一致性检验，即考察这两个方法在检验方面是否吻合
# 或者使用vcd包中的Kappa
library(vcd)
Kappa(data73)

结果如下所示：

结果分析：McNemar配对卡方检验的p值为0.01616,在0.05水平上，有显著性差异，免疫荧光法的总体阳性检出率要高于乳胶凝集法。

Kappa检验结果的p值小于0.01，估计值为0.455，说明两种方法的吻合度有统计学意义，但吻合度一般。

简单来说，McNemar检验是比较两种方法A与B对同一个对象的检验是否有差异，判断哪种方法比较好，比较容易检验出来结果。而Kappa检验就是比较A与B这两种方法检验出来的结果是不是一致的。在有些配对案例中可能会出现A的方法优于B，但A与B检验的结果不一致的情况。

四格表精确Fisher检验

案例4:

某医师为研究乙肝免疫球蛋白预防胎儿宫内感染HBV的效果，将33例HBsAg阳性孕妇随机分为预防注射组和非预防组，结果见表7-4。问两组新生儿的HBV总体感染率有无差别？（《医学统计学》第四版，孙振球，第七章 卡方检验，例7-4)

1
2
3
4
5
6
7
8
9
10
11
data74<-matrix(c(4,5,18,6),nr = 2,
               dimnames = list(c("预防注射组", 
                                 "非预防组"),c("阳性", "阴性")))
chi<-chisq.test(data74) 
# 理论频数小于5的情况，这里用Yates连续校正，
# 默认的就是，不用添加correct = T
chi$expected 
# 查看理论频数，有小于5的情况
fisher.test(data74)

结果如下所示：

案例5:

某单位研究胆囊腺癌、腺瘤的P53基因表达，对同期手术切除的胆囊腺癌、腺瘤标本各10份，用免疫组化法检测P53基因，资料见表7-6。问胆囊腺癌和胆囊腺瘤的P53基因表达阳性率有无差别？（《医学统计学》第四版，孙振球，第七章 卡方检验，例7-5）

1
2
3
4
5
data75<-matrix(c(6,1,4,9),nr = 2,
               dimnames = list(c("胆囊腺癌", "胆囊腺瘤"),
                               c("阳性", "阴性")))
fisher.test(data75)

结果如下所示：

无序R×C列表

案例6:列与行的数目都大于2

无序R×C列表这种情况是行、列都无序，并且行与列的数目都大于2，例如： 某医师研究物理疗法、药物治疗和外用膏药三种疗法治疗周围性面神经麻痹的疗效，资料见表7-8。问三种疗法的有效率有无差别（《医学统计学》第四版，孙振球，第七章 卡方检验，例7-6）？
如下所示：

1
2
3
4
5
data76<-matrix(c(199,164,118,7,18,26),
               nr = 3,dimnames = list(c("物理疗法组", "药物治疗组","外用膏药组"),
                                      c("有效", "无效")))
chi.result <- chisq.test(data76)
chi.result # 查看卡方分析结果

结果如下所示：

双向无序R×C列表的卡方检验

案例7:

某医师在研究血管紧张素I转化酶(ACE)基因I/D多态与2型糖尿病肾病(DN)的关系时，将249例2型糖尿病患者按有无糖尿病肾病分为两组，资料见表7-9，问两组2型糖尿病患者的ACE基因型总体分布有无差别（《医学统计学》第四版，孙振球，第七章 卡方检验，例7-7）？

1
2
3
4
5
6
data77 <- matrix(c(42,30,48,72,21,36),
                 nrow=2,dimnames=list(c("DN组","无DN组"),
                                      c("DD","ID","II")))
chi.result77 <- chisq.test(data77)
chi.result77$expected
chi.result77

结果如下所示：

双向无序R×C列表的关联性检验

案例8

测得某地5801人的ABO血型和MN血型结果如表7-10，问两种血型系统之间是否有关联？（《医学统计学》第四版，孙振球，第七章 卡方检验，例7-8）

没有找到关于Pearson列联系数的函数，手工计算如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
data78 <- matrix(c(431,388,495,137,490,410,587,179,902,800,950,32),nrow=4,dimnames=list(c("O","A","B","AB"),
                                                                                        c("M","N","MN")))
chi.result <- chisq.test(data78)
chi.result
sqrt(chi.result$statistic/(sum(data78)+chi.result$statistic)) # 列联系数
# 根据手工计算结果，可以构建一个计算Pearson列联系数的自定义函数，如下所示：
Pearson_Continuity <- function(x){
  chisq_result <- chisq.test(x)$statistic 
  # 提取计算的卡方值
  Pcoc <- sqrt(chisq_result/(sum(x)+chisq_result))
  print.noquote("Persaon列联系数") 
# print.noquote的作用是显示字符，不带引号，
# 而print则是带引号
  return(as.numeric(Pcoc))}
Pearson_Continuity(data78)

结果如下所示：

两两比较——所有组进行两两比较

当对所有组进行两两比较时，通常采用的方法是Bonferroni法，此方法对原有的检验水平进行调整，即原有的检验水平α除以比较次数得到调整后的检验水平。书本（《医学统计学》第三版，孙振球，P122）中的公式在分母中加了1，不清楚为什么，没有提出采用的是Bonferroni法。而另外两本参考书，即《医学统计学》第6版，马斌荣，P91与《医学统计学》第6版，颜虹，2010年，P155这两本书中所用到的两两比较的方法是Bonferroni法，在下面的题目中使用Bonferroni法。

案例9

对例7-6中的资料进行两两比较，以推断是否任两种疗法治疗周围性面神经麻痹的有效率均有差别？（《医学统计学》第四版，孙振球，第七章 卡方检验，例7-9）

例7-6中共有三组疗法，需要两两比较的次数是3次，因此检验水平是0.05/3=0.0167。

1
2
3
4
5
6
7
8
data76<-matrix(c(199,164,118,7,18,26),
               nr = 3,
               dimnames = list(c("物理疗法组", "药物治疗组","外用膏药组"),
                               c("有效", "无效")))
chisq.test(data76[1:2,],correct=F) 
# 比较物理疗法组与药物治疗组，
# R中默认四格表进行了Yates校正，
# 此处设为F，不进行校正，因为实际值与观测值都大于5

结果如下所示：

从结果可以看出物理疗法组与药物治疗组有统计学差异。

1
2
chisq.test(data76[2:3,], correct=F)
#比较药物治疗组与外用膏药组


从结果来看，药物治疗组与外用膏药组无统计学上的差异。

1
2
chisq.test(data76[-2,],correct=F)
# 比较物理疗法组与外用膏药组

从结果来看，物理疗法组与外用膏药组有统计学上的差异。

两两比较——各个实验组与一个对照组进行比较

两两比较——各个实验组与一个对照组进行比较
当用各个实验组与一个对照组分别进行比较时，需要对检验水平α进行校正，即α/(k-1),其中k是组别。在此案例中，组别k=3,因此，调整后的检验水平=0.05/(3-2)=0.025。

案例10:

以表7-8资料中的药物治疗组为对照组，物理疗法组与外用膏药组为试验组，试分析两试验组与对照组的总体有效率有无差别？（《医学统计学》第四版，孙振球，第七章 卡方检验，例7-10）

1
2
3
4
5
chisq.test(data76[1:2,],correct=F) 
#比较物理疗法组与药物治疗组
chisq.test(data76[2:3,], correct=F) 
#比较药物治疗组与外用膏药组

第一个结果显示，P=0.009343小于0.025，表明物理疗法组与药物治疗组有统计学意义
第二个结果显示，P=0.03214大于0.025，表明药物治疗组与外用膏药组无统计学意义。

如果不使用校正p值的这种方法，可以使用其它的R包来进行两两比较，其中要用到的R包是rcompanion。现在使用这个包来进行两两比较，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
install.packages("rcompanion")
library(rcompanion)
data76<-matrix(c(199,164,118,7,18,26),
               nr = 3,
               dimnames = list(c("物理疗法组", "药物治疗组","外用膏药组"),
                               c("有效", "无效"))
pairwiseNominalIndependence(data76,
                            fisher = FALSE,
                            gtest  = FALSE,
                            chisq  = TRUE,
                            method = "fdr",)

计算结果如下所示：

1
2
3
4
               Comparison  p.Chisq p.adj.Chisq
1 物理疗法组 : 药物治疗组 1.68e-02    0.025200
2 物理疗法组 : 外用膏药组 9.34e-06    0.000028
3 药物治疗组 : 外用膏药组 4.78e-02    0.047800

单向R×C列联表分析——列有序

案例11

某公司研制一种治疗风湿性关节炎的中药，选择试验组和对照组各60例，试验组服用新药，对照组服用常规药物。以中医证侯疗效为主要疗效指标，分为无效、有效、显效、痊愈四个等级。治疗两个疗程后，两组的疗效情况如表5.4所示（略，见代码），欲分析两种药物的疗效是否有差异（《医学案例统计分析与SAS应用》第五章例5.4，冯国双，北京大学医学出版社，2010:93）。

第一种方法：Ridit分析

1
2
3
4
5
6
7
8
9
10
11
12
13
data54 <- c(33,18,24,31,1,9,2,2)
data54m <- matrix(data54,nrow=2,
                  dimnames=list(c("对照组","实验组"),
                                c("无效","有效","显效","痊愈")))
data54m
#从数据来看，列变量，即无效、有效、显效、痊愈这四个等级是有序的，
# 即越来越好，而行变量，对照组与实验组，这是无效的，适合采用Ridit分析：
library(Ridit)
ridit(data54m,1) 
# 函数ridit的参数中1表示分组信息是行，即行是无序的，列是有序的，
# 如果是2，则表示分组信息在列，即列是无序的，行是有序的

第二种方法：CMH分析

用CMH也能进行RC列联表的单向有序分析，CMHtest函数在vcdExtra包中，需要载入

1
2
library(vcdExtra)
CMHtest(data54m)

CMH的单向有序分析参见第3行，即cmeans，这一行表示列有序，行无序的分析结果。

双向有序R×C列联表分析-CMH检验

对于双向有序的RC列联表进行分析，常用于CMH检验，卡方检验常用于分类变量资料的一种假设检验，该法主要用于两个或多个样本率的比较，也可以用于两变量间的关联分析、频属分析的拟合优度检验等。CMH检验是在1959年提出的mh统计方法的基础上提出来的，1988年才发展完善，现在习惯上称之为扩展的MH卡方检验。他们的主要区别在于，如果用同样的方法做同样的实验，由于实验的地点不同，若将这些不同地点的实验数据简单合并后有卡方检验，会参杂很多混杂因素，这会使检验结果产生很大的偏差—主要是由不同实验的点的软、硬条件不同造成的。采用CMH检验可以解决这类问题，采用的是多中心或分层分析方法，多用R×C列联表的分层统计处理1。此外，CMH还应用于双向有序的RC列联表的分析。如下所示：

案例12

某研究者欲研究年龄与冠状动脉粥样硬化等级之间的关系，将278例尸解资料整理成表7-13，问年龄与冠状动脉粥样硬化等级之间是否存在线性变化趋势？（《医学统计学》第四版，孙振球，第七章 卡方检验，例7-9）

第一种方法：

1
2
3
4
5
6
library(vcdExtra) 
# 载入vcdExtra包
data711 <- matrix(c(70,27,16,9,22,24,23,20,4,9,13,15,2,3,7,14),nrow=4)
dimnames(data711) <- list(age=c("20~","30~","40~","≥50"),degree=c("-","+","++","+++"))
CMHtest(data711) # 分析数据
chisq.test(data711)

第二种方法：

通过coin包中的lbl_test函数也能进行CMH双向有序表的计算：

1
2
3
library(coin)
data711t <- as.table(data711)
lbl_test(data711t,scores=list(age=c(1,2,3,4),degree=c(1,2,3,4)))

其中，scores是各个有序变量的得分，与教材中相同。

第三种方法：

1
2
3
4
5
library(reshape2)
datax <- melt(data711)
datax$age <- factor(datax$age,order=T)
datax$degree <- factor(datax$degree,order=T)
CMHtest(xtabs(value~age+degree,data=datax))

结果解读：(1)第1行的cor指的是行与列都为有序变量时的结果，即教材中的线性回归分量，为63.389，这个与教材的结果不一致，但R计算的结果与教材附带的SAS程序计算出来的结果一致，教材可能是印刷错误；(2)第2行的rmeans表示行是有序变量时的结果；(3)第3行的cmeans表示列是有序变量时的结果；(4)general表示列与行都是无序变量时的结果。

常规卡方计算的结果是71.4325。线性回归分量的结果是63.389，二者之差是教材中的偏离线性回归分量，即71.4325-63.389=8.0435，相应的P值为:

1
1-pchisq(8.0435,8)

p值结果为0.4292325，可以认为年龄与冠状动物粥样硬化之间不仅存在相关关系，并且为线性关系，从数据中可以看出来，冠状动脉粥样硬化的等级随着年龄的增加而增高。
注：pchisq(8.0435,8)表示，计算自由度为8，卡方值是8.0435的左侧卡方分布曲线下面积。

第4种方法：

如下所示：

1
2
3
install.packages("multiCA")
library(multiCA)
multiCA.test(data711)

计算结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
> multiCA.test(data711)
$overall
	Multinomial Cochran-Armitage trend test
data:  data711
W = 63.68, df = 3, p-value = 9.611e-14
alternative hypothesis: true slope for outcomes 1:nrow(x) is not equal to 0
$individual
[1] 1.020527e-10 4.248539e-01 1.963178e-02 2.613094e-09
attr(,"method")
[1] "Holm-Shaffer"

CMH检验对多层分类资料的分析

案例13

某公司推出一种治疗颈椎病的治疗仪，为评价其疗效，进行多中心的随机对照试验。共选择两家医院，每家医院60例研究对象，试验组与对照组各30例。以颈痛疗效作为主要疗效指标，根据缓解程度分为有效和无效。治疗20天后，各中心每组的疗效情况如下所示（参见代码），欲分析两种治疗仪的疗效是否有差异（《医学案例统计分析与SAS应用》第五章例5.4，冯国双，北京大学医学出版社，2010:97）。

1
2
3
4
5
6
7
8
9
data56 <- c(8,4,22,26,11,9,19,21)
data56a <- array(data56,dim=c(2,2,2))
data56a
dimnames(data56a) <- list(
  c("对照组","试验组"),
  c("无效","有效"),
  c("中心1","中心2"))
CMHtest(data56a,overall=TRUE)

结果由三部分构成，第一部分是中心1的结果，第二部分是中心2 的结果，第三部分是控制了中心的影响后，两组疗效差异无统计学意义。

拟合优势的卡方检验

案例14：判断数据是否服从泊松分布

观察某克山病区克山病患者的空间分布情况，调查者将该地区划分为279个取样单位，统计各取样单位历年累计病例数，资料见表7-15的第（1）、（2）栏，问此资料是否服从Poisson分布？（《医学统计学》第三版，孙振球，第七章 卡方检验，例7-13）

1
2
3
4
5
6
7
8
9
10
11
12
x<-0:6
y<-c(26,51,75,63,38,17,9)
mean<-mean(rep(x,y))
q<-ppois(x,mean)
n<-length(y)
p<-c()
p[1]<-q[1]
p[n]<-1-q[n-1]
for(i in 2:(n-1))
  p[i]<-q[i]-q[i-1]
chisq.test(y, p=p,correct=F)

案例15 ：判断数据是否服从某一比例

假设H0：总体具有某分布F 备择假设H1：总体不具有该分布。
我们将数轴分成若干个区间，所抽取的样本会分布在这些区间中。在原假设成立的条件下，我们便知道每个区间包含样本的个数的期望值。用实际值Ni 与期望值Npi可以构造统计量K 。皮尔森证明，n趋向于无穷时，k收敛于m-1的塔防分布。m为我们分组的个数。有了这个分布，我们就可以做假设检验。

某消费者协会为了确定市场上消费者对5咎品牌啤酒的喜好情况，随机抽取了1000名啤酒爱好者作为样本进行如下试验：每个人得到5种品牌的啤酒各一瓶，但未标明牌子。这5种啤酒按分别写着A、B、C、D、E字母的5张纸片随机的顺序送给每一个人，下表是根据样本资料整理得的各种品牌啤酒爱好者的频数分页，试根据这些数据判断消费者对这5种品牌啤酒的爱好有无明显差异？

1
2
3
4
5
x <- c(210,312,170,85,223) 
n <- sum(x); m <- length(x) 
p <- rep(1/m,m) 
K <- sum((x-n*p)^2/(n*p)); K #计算出K值 
p <- 1-pchisq(K,m-1); p #计算出p值

在R语言中 chisq.test()，可以完成拟合优度检验。默认就是检验是否为均匀分布，如果是其他分布，需要自己分组，并在参数p中指出。上面题目的解法如下所示：

1
chisq(x)

同样拒绝原假设。

案例16 ：判断数据是否服从正态分布

抽取31名学生的成绩，检验是否为正态分布。

1
2
3
4
5
x <- c(25,45,50,54,55,61,64,68,72,75,75,78,79,81,83,84,84,84,85,86,86,86,87,89,89,89,90,91,91,92,100) 
A <- table(cut(x,breaks=c(0,69,79,89,100))) ;A#对样本数据进行分组。 
p <- pnorm(c(70,80,90,100),mean(x),sd(x)) #获得理论分布的概率值 
p <- c(p[1],p[2]-p[1],p[3]-p[2],1-p[3]) 
chisq.test(A,p=p)


结果显示，其p值小于0.05，成绩不属于正态分布。

再用shirpo.test()和qq图检验一下：

1
2
shapiro.test(x)
qqnorm(x);qqline(x, col = 2)


结果显示，p值小于0.05，不是正态分页，并且QQ图明显表示出成绩不符合正态分布。

案例17 ：判断数据是否符合原比例

大麦杂交后关于芒性的比例应该是 无芒：长芒：短芒=9:3:4 。我们的实际观测值是335：125：160 。请问观测值是否符合预期？

1
2
3
p <- c(9/16,3/16,4/16) 
x <- c(335,125,160) 
chisq.test(x,p=p)

在分组的时候要注意，每组的频数要大于等于5.如果理论分布依赖于多个未知参数，只能先由样本得到参数的估计量。然后构造统计量K，此时K的自由度减少位置参数的数量个。

ks检验

R语言中提供了ks.test()函数，理论上可以检验任何分布。他既能够做单样本检验，也能做双样本检验。

案例18：判断数据是否服从指数分布

记录一台设备无故障工作时常，并从小到大排序420 500 920 1380 1510 1650 1760 2100 2300 2350。问这些时间是否服从拉姆达=1/1500的指数分布？

1
2
x <- c(420,500,920,1380,1510,1650,1760,2100,2300,2350) 
ks.test(x,"pexp",1/1500)

案例19：判断两批数据是否服从相同分布

有两个分布，分别抽样了一些数据，问他们是否服从相同的分布。

1
2
3
4
5
6
7
X<-scan() 
0.61 0.29 0.06 0.59 -1.73 -0.74 0.51 -0.56 0.39 1.64 0.05 -0.06 0.64 -0.82 0.37 1.77 1.09 -1.28 2.36 1.31 1.05 -0.32 -0.40 1.06 -2.47
Y<-scan() 
2.20 1.66 1.38 0.20 0.36 0.00 0.96 1.56 0.44 1.50 -0.30 0.66 2.31 3.29 -0.27 -0.37 0.38 0.70 0.52 -0.71
ks.test(X,Y)

将频数表转换为原始数据

在做卡方分析时，我们经常使用频数表，如下所示：

1
2
3
4
5
6
7
8
> rc <- matrix(c(3,3,4,2),nrow=4,byrow=T)
> rc <- matrix(c(3,3,4,2),nrow=2,byrow=T)
> colnames(rc) <- c("somker","nonsmoker")
> rownames(rc) <- c("g1","g2")
> rc
   somker nonsmoker
g1      3         3
g2      4         2

但有时候我们需要把频数表转换为原始数据，此时可以使用NCStats这个包，这个包无法使用install.packages命令安装，但可以去它的官网进行安装，如下所示：

http://www.rforge.net/NCStats/files

mark

下载最新版的即可，安装，不过这个包还依赖于其它包的安装，根据R的提示安装即可，我在安装过程中，系统提示还要安装FSA和plotrix这两个包，安装后，加载此包，将rc这个数据还原于原始数据，如下所示：

1
2
3
library(NCStats)
rc
frc <- expandTable(rc,c("group","somker status"))

结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
> rc
   somker nonsmoker
g1      3         3
g2      4         2
> frc <- expandTable(rc,c("group","somker status"))
> frc
   group somker status
1     g1        somker
2     g1        somker
3     g1        somker
4     g2        somker
5     g2        somker
6     g2        somker
7     g2        somker
8     g1     nonsmoker
9     g1     nonsmoker
10    g1     nonsmoker
11    g2     nonsmoker
12    g2     nonsmoker

参考资料

  1. http://doctor.51daifu.com/2007/0305/1ADAD919C959T148648.shtml
  2. http://blog.sina.com.cn/s/blog_686dcf990100z9q3.html
  3. http://medstats.blog.163.com/blog/static/17974801620089108192736/
  4. http://blog.sina.com.cn/s/blog_4b361fe00100y222.html
  5. 孙振球. 医学统计学[M]. 人民卫生出版社, 2010.
  6. 徐俊晓. 统计学与R读书笔记(第六版)
  7. 葛毅, 胡良平. 高维列联表资料的统计分析与 SAS 软件实现 (一)[J]. 中西医结合学报, 2009, 7(11): 1086-1089.
  8. http://bbs.pinggu.org/thread-521307-1-1.html↩
  9. 白话统计.冯国双
  10. 医学统计学与R语言（公众号）