R语言笔记-分类变量列联表和独立性检验

2022年6月4日 1,003次阅读 来源: PUMC芋圆四号

文章目录

示例数据：
MASS包中的
birthwt数据集。

首先将数据集中的分类变量因子化，具体参考
这里。

频数表与列联表的生成

一维频数表

  • 仅展现1个分类变量频数统计
  • 语法：table(数据框名$变量名)
> table(birthwt$low)

 no yes 
130  59
  • 展现百分比，而非频数：prop.table(频数表)
> prop.table(table(birthwt$low))

       no       yes 
0.6878307 0.3121693  # 可使用round()保留小数，此处不再赘述
  • epiDisplay包中的tab1()函数：给出含百分比、累计百分比的一维频数表，并生成一个频数分布的条形图；语法为tab1(数据框名$变量名)
> tab1(birthwt$low)
birthwt$low : 
        Frequency Percent Cum. percent
no            130    68.8         68.8
yes            59    31.2        100.0
  Total       189   100.0        100.0

《R语言笔记-分类变量列联表和独立性检验》

二维列联表

  • 又称为交叉表，展现了两个分类变量下的频数统计
  • 也可以使用table函数：table(数据框名$变量1,数据框名$变量2)
> table(birthwt$low,birthwt$smoke)
     
      no yes
  no  86  44
  yes 29  30
  • 生成边际频数 （行列subtotal汇总）：addmargins(频数表)
> addmargins(table(birthwt$low,birthwt$smoke))
     
       no yes Sum
  no   86  44 130
  yes  29  30  59
  Sum 115  74 189
  • epiDisplay包中的tabpct()函数：给出含边际频数、按行&按列求百分比的列联表，还给出一个马赛克图；语法为tabpct(数据框名$变量1,数据框名$变量2)
> tabpct(birthwt$low,birthwt$smoke)

Original table 
           birthwt$smoke
birthwt$low   no  yes  Total
      no      86   44    130
      yes     29   30     59
      Total  115   74    189

Row percent 
           birthwt$smoke
birthwt$low      no     yes  Total
        no       86      44    130
             (66.2)  (33.8)  (100)
        yes      29      30     59
             (49.2)  (50.8)  (100)

Column percent 
           birthwt$smoke
birthwt$low   no       %  yes       %
      no      86  (74.8)   44  (59.5)
      yes     29  (25.2)   30  (40.5)
      Total  115   (100)   74   (100)

《R语言笔记-分类变量列联表和独立性检验》

多维列联表

  • 涉及≥3个分类变量的列联表
  • table(), prop.table(), addmargins()等语法同样适用，不再赘述

独立性检验

χ 2 \chi^2 χ2检验

期望频数列联表生成语法：chisq.test(列联表)$expected，用于评价使用何种检验方法 （见下）。各检验方法具体使用场合总结如下：

名称场合语法
Pearson χ 2 \chi^2 χ2检验二维列联表，每个单元格期望频数≥5chisq.test(列联表,correct=FALSE)
Pearson χ 2 \chi^2 χ2检验带连续性校正二维列联表，存在单元格期望频数≥1且≤5chisq.test(列联表) # correct参数默认为TRUE，故无需声明
Fisher精确概率检验二维列联表，总记录数n<40 或 存在单元格期望频数<1fisher.test(列联表)
> chisq.test(table(birthwt$low,birthwt$smoke))$expected
     
            no      yes
  no  79.10053 50.89947
  yes 35.89947 23.10053 # 期望频数均>5，无需进行连续性校正
> chisq.test(table(birthwt$low,birthwt$smoke),correct=FALSE)

	Pearson's Chi-squared test

data:  table(birthwt$low, birthwt$smoke)
X-squared = 4.9237, df = 1, p-value = 0.02649

相对危险度 (relative risk, RR)与比值比 (odds ratio, OR)

RR与OR的定义不再复习，请参考统计学书籍。epiDisplay包中的cs()和cc()函数可用于计算RR、OR，具体语法如下：

函数语法1语法2
cs()：用于计算RRcs(数据框名$结局变量,数据框名$因素变量)cs(cctable=列联表)
cc()：用于计算ORcc(数据框名$结局变量,数据框名$因素变量)cc(cctable=列联表)

需要注意的是：

  1. 若将列联表作为输入，需要将其作为参数cctable声明
  2. 若将列联表作为输入，列联表必须为结局变量在前、因素变量在后的形式
> cc(birthwt$low,birthwt$smoke) # 语法1

           birthwt$smoke
birthwt$low  no yes Total
      no     86  44   130
      yes    29  30    59
      Total 115  74   189

OR =  2.02 
95% CI =  1.08, 3.78  
Chi-squared = 4.92, 1 d.f., P value = 0.026
Fisher's exact test (2-sided) P value = 0.036 

> mytable<-table(birthwt$low,birthwt$smoke) # 注意该列联表生成时将关心的结局变量放在前面
> cc(cctable=mytable) # 语法2

       
         no yes Total
  no     86  44   130
  yes    29  30    59
  Total 115  74   189

OR =  2.02 
95% CI =  1.08, 3.78  
Chi-squared = 4.92, 1 d.f., P value = 0.026
Fisher's exact test (2-sided) P value = 0.036

cc()除了给出检验结果外，还能生成一个不同组比值 (Odds)的变化图：
《R语言笔记-分类变量列联表和独立性检验》

分层情形下的独立性检验：Mantel-Haenszel检验

M-H检验：检验两分类变量在第三个变量 （分层变量）的调整下是否仍然独立。
列联表写法：table(数据框名$待检验变量1,数据框名$待检验变量2,数据框名$用于分层的变量3)

函数语法
M-H检验，不进行连续性校正mantelhaen.test(列联表,correct=FALSE)
M-H检验，进行连续性校正mantelhaen.test(列联表) # 默认校正，无需声明
M-H检验，连续性校正，并生成OR图 （需要epiDisplay包）mhor(mhtable=列联表)
> mytable<-table(birthwt$low,birthwt$smoke,birthwt$race)
> mantelhaen.test(mytable,correct=FALSE) # 声明不进行连续性校正

	Mantel-Haenszel chi-squared test without continuity
	correction

data:  mytable
Mantel-Haenszel X-squared = 9.4134, df = 1, p-value =
0.002154
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
 1.490740 6.389949
sample estimates:
common odds ratio 
         3.086381 

> mhor(mhtable=mytable)

Stratified analysis by  Var3 
               OR lower lim. upper lim. P value
Var3 white   5.66      1.657      25.14 0.00179 # 可见，白人亚组中吸烟对低体重的OR最大
Var3 black   3.14      0.487      23.45 0.22797
Var3 other   1.25      0.273       5.28 0.75103
M-H combined 3.09      1.491       6.39 0.00215

M-H Chi2(1) = 9.41 , P value = 0.002 
Homogeneity test, chi-squared 2 d.f. = 2.98 , P value = 0.225

mhor()函数还额外生成了一个不同层中OR差异的图：
《R语言笔记-分类变量列联表和独立性检验》

配对列联表的一致性检验：McNemar检验

对同一对象进行两种处理，然后观察这两种处理的分类变量结果之间的一致性。最典型的例子是比较两种诊断方法的一致性。

检验适用场合语法
McNemar检验，无连续性校正–mcnemar.test(列联表,correct=FALSE)
McNemar检验，有连续性校正不一致结果的个案总数<40mcnemar.test(列联表)

以下面的诊断一致性检验作为示例：

诊断试验A(+)诊断试验A(-)
诊断试验B(+)112
诊断试验B(-)1233
> my.matrix<-matrix(c(11,2,12,33),nrow=2)
> mcnemar.test(my.matrix)

	McNemar's Chi-squared test with continuity correction

data:  my.matrix
McNemar's chi-squared = 5.7857, df = 1, p-value = 0.01616
    原文作者：PUMC芋圆四号
    原文地址: https://blog.csdn.net/ethmery/article/details/109210213
    本文转自网络文章，转载此文章仅为分享知识，如有侵权，请联系博主进行删除。
点赞