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

作者: 小程的学习笔记 | 来源:发表于2023-03-28 23:53 被阅读0次

基本原理：卡方检验就是统计样本的实际观测值与理论推断值之间的 $\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语言统计3：独立性检验（卡方检验）

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读