R语言统计3：独立性检验（卡方检验）

作者: cc的生信随记 | 来源:发表于2023-03-28 23:53 被阅读0次

R语言学习（六）基本统计分析--下
python做数理统计
R语言与统计-3：卡方检验
卡方检验
pyspark卡方检验用于特征选择
2020-07-08 R基础绘图+统计
R语言卡方检验
卡方独立性检验
卡方优度检测 (Python 实现) --基于jupyter
R语言卡方检验大全

基本原理：卡方检验就是统计样本的实际观测值与理论推断值之间的 $\color{green}{偏离程度}$ ，实际观测值与理论推断值之间的偏离程度就决定卡方值的大小，如果卡方值越大，二者偏差程度越大；反之，二者偏差越小；若两个值完全相等时，卡方值就为0，表明理论值完全符合。（★ 卡方检验针对 $\color{green}{分类变量}$ ）
通用公式： $x^{2} =∑ \frac {|observed - expected|^{2}} {expected}$
四格卡方值快速计算公式（又称：拟合度公式）：

观察频数	字段2	字段3	合计
row1	a	b	（a+b）
row2	c	d	（c+d）
合计	（a+c）	（b+d）	（a+b+c+d）

$\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $x^{2} = \frac {n(a*d - b*c)^{2}} {(a+b)(c+d)(a+c)(b+d)}$

$\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $\$ $n = (a+b+c+d)$

不同类型卡方检验的选择：
1） $\color{green}{Pearson Chi-Square}$ 条件:
2 * 2四格表：N≥40，且所有理论频数T≥5
R * C表：低于1/5单元格理论频数T＜5

$\$ $\$ 2） $\color{green}{Continuity}$ $\color{green}{Correction}$ 条件
$\$ $\$ 仅适用于2 * 2四格表：当N≥40但有1≤理论频数T＜5时
$\$ $\$ R * C表: 没有校正的卡方检验

$\$ $\$ 3） $\color{green}{Fisher's Exact Test}$ 条件
$\$ $\$ 当N＜40时或当有理论频数T＜1时使用。
$\$ $\$ SPSS默认在四格表中输出，R * C表需要手动命令后才提供。

$\$ $\$ 4） $\color{green}{Likelihood Ratio}$ 条件
$\$ $\$ 当N＞40，最小期望频数＞5时，结论与皮尔逊卡方基本一致

四格表资料的卡方检验:

# 创建数据集
ID <- seq(1,200)
treat <- c(rep("treated",104),rep("placebo",96))
treat <- factor(treat)
impro <- c(rep("marked",99),rep("none",5),rep("marked",75),rep("none",21))
impro <- as.factor(impro)
data1 <- data.frame(ID,treat,impro)

table(data1$treat,data1$impro)
##          
##           marked none
##   placebo     75   21
##   treated     99    5

示例1 - $\color{green}{CrossTable()}$

# CrossTable()函数可以直接对原始数据记录进行交叉表创建，以及卡方检验
library(gmodels)
CrossTable(data1$treat, data1$impro, expected = T, chisq = T, fisher = T, mcnemar = T, format = "SPSS")

Chi-squared-1

★ 本例符合pearson卡方，卡方值为12.85707，p<0.01（拒绝原假设，说明treat和impro相关）。
★ 四格表资料卡方检验的专用公式/四格表资料卡方检验的校正公式/配对四格表资料的卡方检验/四格表资料的Fisher精确概率法，均可用此方法

示例2 - 把数据变成2x2列联表，然后用 $\color{green}{chisq.test()}$ 函数做卡方检验

mytable <- table(data1$treat,data1$impro)
mytable
##          
##           marked none
##   placebo     75   21
##   treated     99    5

chisq.test(mytable,correct = F)
## 
##  Pearson's Chi-squared test
## 
## data:  mytable
## X-squared = 12.857, df = 1, p-value = 0.0003362

配对四格表资料的卡方检验（ $\color{green}{McNemar}$ 检验）:

ana <- matrix(c(11,12,2,33), nrow = 2, byrow = T,
              dimnames = list(免疫荧光 = c("阳性","阴性"),
                              乳胶凝集 = c("阳性","阴性"))
              )

ana
##         乳胶凝集
## 免疫荧光 阳性 阴性
##     阳性   11   12
##     阴性    2   33

mcnemar.test(ana, correct = T)
## 
##  McNemar's Chi-squared test with continuity correction
## 
## data:  ana
## McNemar's chi-squared = 5.7857, df = 1, p-value = 0.01616

★ 当b(免疫荧光) + c(免疫荧光) < 40时，使用连续性校正，即correct=TRUE。
★ 当b(免疫荧光) + c(免疫荧光) ≥ 40时，不使用连续性校正，即correct=FALSE。

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

hbv <- matrix(c(4,18,5,6), nrow = 2, byrow = T,
              dimnames = list(组别 = c("预防注射组","非预防组"),
                              效果 = c("阳性","阴性"))
              )
hbv
##             效果
## 组别       阳性 阴性
##   预防注射组    4   18
##   非预防组      5    6

fisher.test(hbv)
## 
##  Fisher's Exact Test for Count Data
## 
## data:  hbv
## p-value = 0.121
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.03974151 1.76726409
## sample estimates:
## odds ratio 
##  0.2791061

行 x 列表资料的卡方检验：

示例1 - $\color{green}{多个样本率}$ 的比较，首先要对数据格式转换一下，变成table或者matrix

M <- as.table(rbind(c(199, 7), 
                    c(164, 18),
                    c(118, 26)))
M
##       effect
## trt    有效 无效
##   物理  199    7
##   药物  164   18
##   外用  118   26

dimnames(M) <- list(trt = c("物理", "药物", "外用"),
                    effect = c("有效","无效"))

chisq.test(M, correct = F)
## 
##  Pearson's Chi-squared test
## 
## data:  M
## X-squared = 21.038, df = 2, p-value = 2.702e-05

# 物理治疗组和药物治疗组的卡方检验
chisq.test(M[1:2,], correct = F)
## 
##  Pearson's Chi-squared test
## 
## data:  M[1:2, ]
## X-squared = 6.756, df = 1, p-value = 0.009343

# 物理治疗组和外用膏药组的卡方检验
chisq.test(M[c(1,3),], correct = F)
## 
##  Pearson's Chi-squared test
## 
## data:  M[c(1, 3), ]
## X-squared = 21.323, df = 1, p-value = 3.881e-06

# 药物治疗组和外用膏药组的卡方检验
chisq.test(M[2:3,], correct = F)
## 
##  Pearson's Chi-squared test
## 
## data:  M[2:3, ]
## X-squared = 4.591, df = 1, p-value = 0.03214

示例2 - $\color{green}{样本构成比}$ 的比较

ace <- matrix(c(42,48,21,30,72,36),nrow = 2,byrow = T,
              dimnames = list(dn = c("dn组","非dn组"),
                              idi = c("dd","id","ii"))
              )
ace
##         idi
## dn       dd id ii
##   dn组   42 48 21
##   非dn组 30 72 36

chisq.test(ace, correct = F)
## 
##  Pearson's Chi-squared test
## 
## data:  ace
## X-squared = 7.9127, df = 2, p-value = 0.01913

Cochran-Mantel-Haenszel 卡方统计量检验(分层卡方检验）

简称CMH检验，主要用于高维列联表的分析。通过对分层因素进行控制，从而考察调整混杂因素之后暴露/处理因素与结局事件之间的关联性

假设检验：H0：为任一层的行变量X与列变量Y均不相关； H1：为至少有一层X与Y存在统计学关联。
当H0成立时，CMH统计量渐近卡方分布。需要注意的是，当各层间行变量与列变量相关的方向不一致时，CMH统计量的检验效能较低。

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

myo <- array(c(17,47,
               121,944,
               12,158,
               14,663),
             dim = c(2,2,2),
             dimnames = list(心肌梗死 = c("病例","对照"),
                             口服避孕药 = c("是","否"),
                             年龄分层 = c("<40岁","≥40岁")
                             )
             )
myo
## , , 年龄分层 = <40岁
## 
##         口服避孕药
## 心肌梗死 是  否
##     病例 17 121
##     对照 47 944
## 
## , , 年龄分层 = ≥40岁
## 
##         口服避孕药
## 心肌梗死  是  否
##     病例  12  14
##     对照 158 663

mantelhaen.test(myo,correct = F)
## 
##  Mantel-Haenszel chi-squared test without continuity correction
## 
## data:  myo
## Mantel-Haenszel X-squared = 24.184, df = 1, p-value = 8.755e-07
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.930775 4.933840
## sample estimates:
## common odds ratio 
##          3.086444

★ p-value = 8.755e-07 (p<0.001)，按a=0.05的检验水准拒绝H0，接受H1，可认为控制了年龄的影响后，心肌梗死与近期服用口服避孕药有关

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

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

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

# Breslow-Day对各层的效应值进行齐性检验
library(DescTools)
## Registered S3 method overwritten by 'DescTools':
##   method         from 
##   reorder.factor gdata

BreslowDayTest(myo)
## 
##  Breslow-Day test on Homogeneity of Odds Ratios
## 
## data:  myo
## X-squared = 0.23409, df = 1, p-value = 0.6285

★ p-value = 0.6285，可认为两年龄组口服避孕药对心肌梗死的总体OR值同质

# woolf法检验不同分层之间的效应值的同质性
woolf <- function(x) {
  x <- x + 1 / 2
  k <- dim(x)[3]
  or <- apply(x, 3, function(x) (x[1,1]*x[2,2])/(x[1,2]*x[2,1]))
  w <-  apply(x, 3, function(x) 1 / sum(1 / x))
  1 - pchisq(sum(w * (log(or) - weighted.mean(log(or), w)) ^ 2), k - 1)
}

woolf(myo)
## [1] 0.6400154

★ p-value = 0.6400154，和BreslowDayTest差不多

参考：
1.R语言和医学统计学(3): 卡方检验 (https://blog.csdn.net/Ayue0616/article/details/127587730)

R语言卡方检验大全 (https://www.jianshu.com/p/f1abe50b819b)

R语言学习（六）基本统计分析--下
独立性检验卡方独立性检验使用chisq.test()函数对二维表的行变量和列变量进行卡方独立性检验，代码如下 ...
python做数理统计
卡方分析--独立性检验
R语言与统计-3：卡方检验
R语言与统计-1：t检验与秩和检验[https://www.jianshu.com/p/ba629f6ae85d]...
卡方检验
白话统计学—卡方检验基本原理R语言实现卡方检验的替换组内两两比较等级资料的比较单向R×C列联表分析——列有序双向有...
pyspark卡方检验用于特征选择
卡方检验特征选择原理计算特征变量与因变量的卡方独立性检验统计量，如果特征变量与因变量独立，说明其对预测因变量效果...
2020-07-08 R基础绘图+统计
R统计安装加载必要的R包 QUALITATIVE DATA QUANTITATIVE DATA 卡方检验，Fis...
R语言卡方检验
卡方检验是一种确定两个分类变量之间是否存在显着相关性的统计方法。这两个变量应该来自相同的人口，他们应该是类似 -...
卡方独立性检验
参考书《白话统计学》非参数检验卡方独立性检验适用于数据来自两个分类的变量，样本对象通过分类变量分成不同类型，...
卡方优度检测 (Python 实现) --基于jupyter
卡方独立性检验（2）参数说明【输入】：[http://localhost:8888/notebooks/%E5...
R语言卡方检验大全
本文首发于公众号：医学和生信笔记医学和生信笔记，专注R语言在临床医学中的使用，R语言数据分析和可视化。主要分享R...