卡方检验在计数资料中的应用，包括推断两个总体率或构成比之间有无差别、多个总体率或构成比之间有无差别、多个样本率间的多重比较、两个分类变量之间有无关联性、多维列联表的分析和频数分布拟合优度的卡方检验。

下面分别介绍计数资料怎么进行卡方检验。

---

目  录

• 1. 四格表资料的卡方检验

  • 1.1 数据集数据

  • 1.2 向量或矩阵数据

  • 1.3 chisq.test()函数

• 2. Fisher确切概率法

• 3. 配对四格表资料的卡方检验

  • 3.1 b+c < 40

  • 3.2 b+c ≥ 40

  • 3.3 mcnemar.test()函数

• 4. Cochran-Mantel-Haenszel检验

  • 4.1 mantelhaen.test()函数

  • 4.2 数据集形式的数据

  • 4.3 向量或数组数据

  • 4.4 计算各层OR值

  • 4.5 Breslow-Day检验

  • 4.6 BreslowDayTest()函数

  • 4.7 WoolfTest()函数

• 5. 计算列联系数

• 6. 频数分布拟合优度的卡方检验

---

1. 四格表资料的卡方检验

1.1 数据集数据

选用survival包的colon数据集。

library(survival) #加载内置数据集的包data(colon) # 加载数据集

mytable chisq.test(mytable)  # 进行连续性校正chisq.test(mytable, correct = FALSE) # 不进行连续性校正

无论是否进行连续性校正，结果都显示p值＞0.05，接受原假设，说明sex和status无关。

还可以输出卡方检验摘要：

chisq.test(mytable)$observed   # 实际频数(和mytable一样)chisq.test(mytable)$expected   # 期望频数chisq.test(mytable)$residuals  # Pearson 残差chisq.test(mytable)$stdres     # 标准化残差

1.2 向量或矩阵数据

compare706,792,184,176), nr = 2,                dimnames = list(c("male", "female"),                                c("Yes", "No")));comparechisq.test(compare)

注意：在使用chisq.test()函数计算时，要注意单元格的期望频数。如果所有单元格频数都不为零，并且所有单元的期望频数≥5，那么Pearson卡方检验是合理的，否则会显示警告信息。

如果在计算时出现警告信息，表明表中单元格期望频数有＜5的值，这可能会使卡方近似无效。

如果数据不满足卡方检验的条件时，应使用Fisher精确检验。

1.3 chisq.test()函数

chisq.test()函数可以进行卡方列联表检验和拟合优度检验。

chisq.test(x, y = NULL, # x是由数据构成的向量或矩阵，y是数据向量(当x为矩阵时，y忽略)           correct = TRUE, # 逻辑词，默认为TRUE，在计算2x2列联表的检验统计量时是否使用连续性校正           p = rep(1/length(x), length(x)), # p是原假设落在小区间的理论概率，默认值表示均匀分布.            # P是和x长度相同的概率向量，P值如果输入有负值，会返回错误信息。           rescale.p = FALSE, # 逻辑词，为TRUE，则p将缩放成和为1的向量；为FALSE，如果p向量和≠1，则会返回错误信息           simulate.p.value = FALSE, # 为TRUE则不进行连续性校正，根据蒙特卡洛检验计算p值           B = 2000) # B为整数，在蒙特卡洛检验中使用的重复次数。

对于2x2的列联表，参数correct的默认值为TRUE，即使用Yate连续修正，目的是提高P值，避免"有显著差异"不可靠的情况发生；

2. Fisher确切概率法

在样本较小时(单元的期望频数＜5)，需要用Fisher精确检验来完成独立性检验。Fisher检验开始是针对2x2四格表提出的，当卡方检验的条件不满足时，可以使用精确检验。

fisher.test(x, y = NULL, # 参数x为二维列联表形式的矩阵，或者由因子构成的对象。y为因子构成的向量，当x为矩阵时，此值无效。            workspace = 200000, # 正整数，表示用于网络算法工作空间的大小            hybrid = FALSE,  # 逻辑词，仅用于2x2列联表，为FALSE时(默认值)表示精确计算概率，取TRUE表示用混合算法计算概率。            hybridPars = c(expect = 5, percent = 80, Emin = 1),            control = list(), or = 1, # 为列表，指定低水平算法的组成，or为优势比的原假设，默认值为1，仅用于2x2列联表。            alternative = "two.sided", # 备择假设，默认值"two.sided"，表示双侧检验(不独立)，“less”表示单侧小于检验(负相关)；“greater”表示单侧大于检验(正相关)。            conf.int = TRUE, # 逻辑词，为TRUE则给出优势比的置信区间。            conf.level = 0.95,  # 置信水平，默认0.95            simulate.p.value = FALSE, B = 2000) # 逻辑词，为TRUE，表示用Monto Carlo方法计算P值，B为正整数，表示Monto Carlo重复的次数。

例：研究COPD与吸烟之间的关系，调查了52名COPD患者和33名非COPD患者吸烟的人数。

compare49,28,3,5),nr = 2, byrow = TRUE,                dimnames = list(c("smoke", "Notsmoke"),                                c("COPD", "normal")))comparechisq.test(compare)$expected  # 查看期望频数fisher.test(compare) # 进行fisher精确检验

因为p-value=0.2518＞0.05,并且优势比的置信区间包含1，故接受原假设, 认为吸烟与COPD无关。

3. 配对四格表资料的卡方检验

计数资料的配对设计常用于两种检验方法、培养方法、诊断方法的比较。特点是对样本中各观察单位分别用两种方法处理，然后观察两种处理方法的某两分类变量的计数结果。

假设检验为McNemar卡方检验，检验统计量有两种(分b+c ≥ 40和b+c < 40的情况)。该法一般用于样本含量不太大的资料。

• 当b+c < 40时，使用连续性校正，即correct=TRUE。
• 当b+c ≥ 40时，不使用连续性校正，即correct=FALSE。

3.1 b+c < 40

例：某实验室分别用乳胶凝集法和免疫荧光法对58名可疑系统性红斑狼疮患者血清中抗核抗体进行测定，问这两种方法的检测结果有无差别？

X 11,2,12,33),             nrow = 2,            byrow = FALSE,            dimnames = list(免疫荧光法 = c("+","-"),乳胶凝集法 = c("+","-")));X#           乳胶凝集法# 免疫荧光法   +  -#          + 11 12#          -  2 33

mcnemar.test(X,correct=TRUE)# McNemar's Chi-squared test with continuity correction# data:  X# McNemar's chi-squared = 5.7857, df = 1, p-value = 0.01616结论：p-value = 0.01616 ＜ 0.05，可以认为两种方法的检测结果不同，免疫荧光法的阳性检测率较高

3.2 b+c ≥ 40

例：某医院同时用A、B两种方法测定160份痰标本中的抗酸杆菌，问A、B两种方法的检出率有无显著性差异？

X 52,20,35,53),             nrow = 2,            byrow = FALSE,            dimnames = list(A方法 = c("+","-"),B方法 = c("+","-")));X#       B方法# A方法  +  -#     + 52 35#     - 20 53mcnemar.test(X,correct=FALSE) # McNemar's Chi-squared test# data:  X# McNemar's chi-squared = 4.0909, df = 1, p-value = 0.04311结论：p-value = 0.04311 ＜ 0.05，可以认为两种方法的检测结果不同。

3.3 mcnemar.test()函数

使用mcnemar.test()函数进行McNemar检验。

mcnemar.test(x, y=NULL, #参数x为二维列联表形式的矩阵，或者由因子构成的对象。y为因子构成的向量，当x为矩阵时，此值无效。             correct=TRUE)  #逻辑词，默认为TRUE，仅在2x2列联表下进行连续校正。

4. Cochran-Mantel-Haenszel检验

CMH检验可以理解为分层卡方检验，CMH检验用于高维列联表的分析，即在控制了某一个或几个混杂因素(分层变量)之后，检验二维RxC表中行变量与列变量之间是否存在统计学关联。

假设检验：H0：为任一层的行变量X与列变量Y均不相关；  H1：为至少有一层X与Y存在统计学关联。

当H0成立时，CMH统计量渐近卡方分布。需要注意的是，当各层间行变量与列变量相关的方向不一致时，CMH统计量的检验效能较低。

根据行变量X和列变量Y的类型不同，CMH卡方统计量包括：

1.相关统计量：适用于双向有序分类变量；  2.方差分析统计量：也称行平均得分统计量，适用于列变量Y为有序分类变量；  3.一般关联统计量：适用双向无序分类变量，目的是检验X和Y是否存在关联性。

4.1 mantelhaen.test()函数

mantelhaen.test()函数用来进行CMH卡方检验。其原假设是，两个分类变量在第三个变量的每一层中都是条件独立的。

mantelhaen.test(x, y = NULL, # x为数组形式的三维列联表，行与列的维度至少为2，最后一个维度为分层变量;                # 或者x为至少有2个水平的因子，y也为至少有2个水平的因子，如果x为三维数组时，y忽略;                z = NULL, # 至少有2个水平的因子，表示哪一层对应x中的元素，哪一层对应y的元素，如果x为三维数组，则z忽略;                alternative = c("two.sided", "less", "greater"), # 备择假设，默认为"two.sided"双侧检验;                correct = TRUE, # 逻辑词，计算检验统计量时使用连续性校正;                exact = FALSE, # 逻辑词，是否计算精确检验;                conf.level = 0.95) # 置信水平，默认0.95

注意：数组或向量不允许存在缺失值，x，y，z必须为长度相同的数字向量。

4.2 数据集形式的数据

选用survival包的colon数据集。xtabs()函数可以基于三个或多个分类变量生成多维列联表。

mytable mantelhaen.test(mytable)

结果中p-value = 2.921e-08，表明患者接受的治疗与结局状态在性别的每一水平下并不独立。

4.3 向量或数组数据

为研究心肌梗死与近期使用避孕药之间的关系，采用病例对照研究方法调查了234 名心肌梗死病人与1742名对照使用口服避孕药的状况。考虑到年龄是混杂因素，按照年龄分层后结果见下表。试分析排除了年龄因素的影响后，心肌梗死是否与近期使用口服避孕药有关。

mydata 17, 121, 47, 944,                  12, 14, 158, 663),                dim = c(2, 2, 2),                dimnames = list(心肌梗死 = c("病例", "对照"),                                    近期使用口服避孕药 = c("是", "否"),                                     年龄 = c(", "≥ 40岁")))mydatamantelhaen.test(mydata, correct = F) # 经典的2x2xk水平 CMH 检验

结果显示P<0.001，按a=0.05的检验水准拒绝H0，接受H1，可认为控制了年龄的影响后，心肌梗死与近期服用口服避孕药有关。

4.4 计算各层OR值

apply(mydata, 3, function(x) (x[1,1]*x[2,2])/(x[1,2]*x[2,1]))#  # 2.821874 3.596745 

4.5 Breslow-Day检验

对于分层病例对照研究或队列研究资料，通常应用Breslow-Day检验对各层的效应值(OR或RR)进行齐性检验。

若不拒绝齐性假设(p＞0.05)，才可依据CMH检验的结果推断出暴露因素是否与疾病相关。如果相关，可进一步用Mantel-Haenszel法估计OR或RR值及其可信区间。

若拒绝了齐性假设(p<0.05)，则提示分层变量与暴露因素间存在交互作用，此时CMH检验的结果不能说明问题，可进行多元logistic回归分析。

4.6 BreslowDayTest()函数

BreslowDayTest(x, # 2x2xk列联表               OR = NA, # 要检验的OR值,默认使用Mantel-Haenszel估算值。               correct = FALSE) # 为TRUE，则使用Tarone调整的Breslow-Day检验

library(DescTools) # 加载包BreslowDayTest(mydata)  # 进行Breslow-Day检验# Breslow-Day test on Homogeneity of Odds Ratios# data:  mydata# X-squared = 0.23409, df = 1, p-value = 0.6285BreslowDayTest(mydata, correct = TRUE)  # Tarone校正的Breslow-Day检验# Breslow-Day Test on Homogeneity of Odds Ratios (with Tarone correction)# data:  mydata# X-squared = 0.23369, df = 1, p-value = 0.6288

结果为P=0.6285，可认为两年龄组口服避孕药对心肌梗死的总体OR值同质：用Mantel-Haenszel法估计OR值及其95%可信区间为3.09(1.93-4.93)。

4.7 WoolfTest()函数

除了BreslowDayTest()函数外，WoolfTest()函数也可以对2x2xk列联表的同质性进行检验。

WoolfTest(x)  # x为2x2xk列联表，最后一个维度为分层变量WoolfTest(mydata)# Woolf Test on Homogeneity of Odds Ratios (no 3-Way assoc.)# data:  mydata# X-squared = 0.23358, df = 1, p-value = 0.6289

结论和BreslowDayTest检验一样。

5. 计算列联系数

若行x列表资料的两个分类变量都为无序分类变量，进行多个样本率的比较、样本构成比的比较都可以用卡方检验，也就是使用chisq.test()函数来计算。

若得知两个分类变量之间有关联性，需进一步分析关系的密切程度时，可计算Pearson列联系数。

列联系数C (Contingency Coeff) 取值为0-1之间，0表示完全独立, 1表示完全相关；愈接近于1，关系愈密切。

vcd包中的assocstats()函数可以用来计算二维列联表的phi系数、列联系数和Cramer’s V系数。

例：测得某地 5801 人的 ABO 血型和 MN 血型结果如下，问这两种血型之间是否有关联。

mytable 431, 388, 495, 137,                    490, 410, 587, 179,                    902, 800, 950, 32),                  nrow = 4, byrow = FALSE,                  dimnames = list(ABO血型 = c("O","A","B","AB"),                                  MN血型 = c("M","N","MN")))mytablechisq.test(mytable) # 判断两分类变量有无关联性

p-value < 2.2e-16，可以认为两种血型系统间有关联，可进一步计算Pearson列联系数，以分析其密切程度。

library(vcd)assocstats(mytable) # 查看关联性的强弱

从上面可以看出，列联系数为0.188，虽然有关联性，但列联系数比较小，虽然有统计学意义，但关系不太密切。

6. 频数分布拟合优度的卡方检验

频数分布拟合优度(goodness-of-fit)检验还是使用chisq.test()函数。

chisq.test(x, # x为向量或单行/单列矩阵，参数 x 必须全部是非负整数。           p, # P是和x长度相同的概率向量，P值如果输入有负值，会返回错误信息。              # p中向量的和应为1，p不给出则向量值全部相等。           rescale.p = TRUE) # 逻辑词，为TRUE，则p将缩放成和为1的向量；为FALSE,如果和≠1，则会返回错误信息。

x 10,25,15) chisq.test(x)  # chisq.test(as.table(x))等效

x 82,47,20,18,22) # x为向量p 45,25,25,18,15) # p为概率chisq.test(x, p = p)  # 返回错误信息,概率总和必须为1chisq.test(x, p = p, rescale.p = TRUE) # 运行正常

---

相关文章

R语言统计与绘图：ROC曲线上32种截断值的计算 2020-06-18 R语言统计与绘图：ROC曲线的统计计算 2020-06-17 R语言统计与绘图：绘制不同坐标轴置信区间的森林图 2020-06-05 R语言统计与绘图：方差分析 2020-06-01 R语言统计与绘图：pROC包绘制ROC曲线 2020-05-25