白话统计学—卡方检验
基本原理
当我们要进行分类资料和分类资料之间的比较的时候可以使用卡方检验。
卡方检验的基本思想其实就是，在我们形成四格表的时候，比较如果没有差异的各个格的格数和现实的格数是否存在差异。如果存在差异的就是有意义的。
R语言实现
R语言要进行卡方检验的时候，可以用到chisq.test。在进行卡方检验之前，我们可以需要先形成一个表格，然后才可以进行卡方检验
library(MASS)
car.data <- data.frame(Cars93$AirBags, Cars93$Type)
target <- table(car.data)
chisq.test(target)
## Warning in chisq.test(target): Chi-squared approximation may be incorrect
##
## Pearson's Chi-squared test
##
## data: target
## X-squared = 33.001, df = 10, p-value = 0.0002723
卡方检验的替换
当总例数> 40且理论频数 < 5的时候，建议使用Yates较正的卡方检验。当总例数 < 40 或理论频数 < 1的时候，使用fisher精准检验更加稳妥一些。
R当中的chisq.test函数实现了Yates较正。通过fisher.test可以实现fisher检验
ca = matrix(c(20,2,14,3), ncol = 2)
chisq.test(ca, correct = T)
## Warning in chisq.test(ca, correct = T): Chi-squared approximation may be
## incorrect
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: ca
## X-squared = 0.095844, df = 1, p-value = 0.7569
fisher.test(ca)
##
## Fisher's Exact Test for Count Data
##
## data: ca
## p-value = 0.6362
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.2112916 28.2603619
## sample estimates:
## odds ratio
## 2.100606
组内两两比较
如果分组有很多种需要多组比较的话，可以使用rcompanion::pairwiseNominalIndependence
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")
## Comparison p.Chisq p.adj.Chisq
## 1 物理疗法组 : 药物治疗组 1.68e-02 0.025200
## 2 物理疗法组 : 外用膏药组 9.34e-06 0.000028
## 3 药物治疗组 : 外用膏药组 4.78e-02 0.047800
等级资料的比较
当分组其实含有等级的变量的时候，如果我们想要把结局当作一个有序变量的话，需要使用秩和检验。如果是不想当作有序变量的话，则使用chisq.test
单向R×C列联表分析——列有序
如果列联表是单项有序的。可以使用Ridit::ridit进行检验
data54 <- c(33,18,24,31,1,9,2,2)
data54m <- matrix(data54,nrow=2,
dimnames=list(c("对照组","实验组"),
c("无效","有效","显效","痊愈")))
data54m
## 无效 有效 显效 痊愈
## 对照组 33 24 1 2
## 实验组 18 31 9 2
#从数据来看，列变量，即无效、有效、显效、痊愈这四个等级是有序的，
# 即越来越好，而行变量，对照组与实验组，这是无效的，适合采用Ridit分析：
library(Ridit)
ridit(data54m,1)
##
## Ridit Analysis:
##
## Group Label Mean Ridit
## ----- ----- ----------
## 1 对照组 0.4267
## 2 实验组 0.5733
##
## Reference: Total of all groups
## chi-squared = 9.2758, df = 1, p-value = 0.002322
# 函数ridit的参数中1表示分组信息是行，即行是无序的，列是有序的
# 如果是2，则表示分组信息在列，即列是无序的，行是有序的
双向有序R×C列联表分析-CMH检验
对于双向有序的RC列联表进行分析，常用于CMH检验，卡方检验常用于分类变量资料的一种假设检验，该法主要用于两个或多个样本率的比较，也可以用于两变量间的关联分析、频属分析的拟合优度检验等。CMH检验是在1959年提出的mh统计方法的基础上提出来的，1988年才发展完善，现在习惯上称之为扩展的MH卡方检验。他们的主要区别在于，如果用同样的方法做同样的实验，由于实验的地点不同，若将这些不同地点的实验数据简单合并后有卡方检验，会参杂很多混杂因素，这会使检验结果产生很大的偏差—主要是由不同实验的点的软、硬条件不同造成的。采用CMH检验可以解决这类问题，采用的是多中心或分层分析方法，多用R×C列联表的分层统计处理1。此外，CMH还应用于双向有序的RC列联表的分析。
R中可以使用vcdExtra::CMHtest进行检验
library(vcdExtra)
## Loading required package: vcd
## Loading required package: grid
## Loading required package: gnm
# 载入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) # 分析数据
## Cochran-Mantel-Haenszel Statistics for age by degree
##
## AltHypothesis Chisq Df Prob
## cor Nonzero correlation 63.389 1 1.6962e-15
## rmeans Row mean scores differ 63.450 3 1.0758e-13
## cmeans Col mean scores differ 65.669 3 3.6072e-14
## general General association 71.176 9 8.9519e-12
结果提过了四种检验结果。第一行代表行和列都是有序变量的结果；第二行代表行是有序变量的结果；第三行代表列是有序变量的结果；第四行代表行和列都是无序变量的结果
CMH检验对多层分类资料的分析
由于CMH检验主要还是用于多分层的检验。因此可以继续检验结果
data56 <- c(8,4,22,26,11,9,19,21)
data56a <- array(data56,dim=c(2,2,2))
data56a
## , , 1
##
## [,1] [,2]
## [1,] 8 22
## [2,] 4 26
##
## , , 2
##
## [,1] [,2]
## [1,] 11 19
## [2,] 9 21
dimnames(data56a) <- list(
c("对照组","试验组"),
c("无效","有效"),
c("中心1","中心2"))
CMHtest(data56a,overall=TRUE)
## $`:中心1`
## Cochran-Mantel-Haenszel Statistics
## in stratum :中心1
##
## AltHypothesis Chisq Df Prob
## cor Nonzero correlation 1.6389 1 0.20048
## rmeans Row mean scores differ 1.6389 1 0.20048
## cmeans Col mean scores differ 1.6389 1 0.20048
## general General association 1.6389 1 0.20048
##
##
## $`:中心2`
## Cochran-Mantel-Haenszel Statistics
## in stratum :中心2
##
## AltHypothesis Chisq Df Prob
## cor Nonzero correlation 0.295 1 0.58703
## rmeans Row mean scores differ 0.295 1 0.58703
## cmeans Col mean scores differ 0.295 1 0.58703
## general General association 0.295 1 0.58703
##
##
## $ALL
## Cochran-Mantel-Haenszel Statistics
## Overall tests, controlling for all strata
##
## AltHypothesis Chisq Df Prob
## cor Nonzero correlation 1.5436 1 0.21408
## rmeans Row mean scores differ 1.5436 1 0.21408
## cmeans Col mean scores differ 1.5436 1 0.21408
## general General association 1.5436 1 0.21408
结果现实多个中心各自的结果，然后会现实进行较正后的结果。这个其实类似于去除批次效应。
Cochran-Armitage趋势检验
前提条件
- 分组变量必须是三组或者以上才能使用的。两组无所谓是否有趋势
- 分组必须是有序变量，无序变量使用意义不大。例如：随着血型的变化，阳性率升高。这样的结果没有意义
- 结局变量必须是二分类的。
- CA检验只能检验是否存在线性关系，对于不是线性的，不能检验出来
R语言实现
R中可以通过prop.trend.test来实现。 其中P < 0.05代表存在线性关系。
smokers <- c( 83, 90, 129, 70 )
patients <- c( 86, 93, 136, 82 ) ###提供每组的总例数
prop.trend.test(smokers, patients)
##
## Chi-squared Test for Trend in Proportions
##
## data: smokers out of patients ,
## using scores: 1 2 3 4
## X-squared = 8.2249, df = 1, p-value = 0.004132
但是上述的分析没办法看出是怎么一个趋势。这里我们可以使用DescTools::CochranArmitageTest结果当中，如果Z是正的则是上升趋势，如果Z是负的则是下降趋势。
lungtumor <- data.frame(dose = rep(c(0, 1, 2), c(40, 50, 48)),
tumor = c(rep(c(0, 1), c(38, 2)),
rep(c(0, 1), c(43, 7)),
rep(c(0, 1), c(33, 15))))
DescTools::CochranArmitageTest(table(lungtumor))
##
## Cochran-Armitage test for trend
##
## data: table(lungtumor)
## Z = -3.2735, dim = 3, p-value = 0.001062
## alternative hypothesis: two.sided
分类变量的赋值对结果的影响
- 如果是无序资料的话，分类变量的赋值是不影响结果的。如果是有序资料的分析则不然，分类变量的赋值是会影响其变化的。例如:MH卡方检验以及CA趋势检验。
- 如果进行合理的赋值呢？最好还是利用图形简单的探索一下各个部分如果。然后根据分布来进行相对应的赋值。一般来说直接看柱状图的分布即可
prop.trend.test中提供了赋值的选项。可以进行不同的赋值
prop.trend.test(smokers, patients, c(0,0,0,1))
##
## Chi-squared Test for Trend in Proportions
##
## data: smokers out of patients ,
## using scores: 0 0 0 1
## X-squared = 12.173, df = 1, p-value = 0.0004848
讲频数表转换为原始数据
有时候我们需要把频数表转换为原始数据。这个时候我们可以使用NCStats::expandTable PS：这个包不能通过install.packages安装。需要通过source("http://www.rforge.net/NCStats/InstallNCStats.R")安装
library(NCStats)
## Loading required package: FSA
## ## FSA v0.8.22\. See citation('FSA') if used in publication.
## ## Run fishR() for related website and fishR('IFAR') for related book.
##
## Attaching package: 'FSA'
## The following object is masked from 'package:gnm':
##
## se
## ## NCStats v0.4.7 by Derek H. Ogle, Northland College.
##
## Attaching package: 'NCStats'
## The following object is masked from 'package:FSA':
##
## rSquared
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
frc <- expandTable(rc,c("group","somker status"))
head(frc)
## group somker status
## 1 g1 somker
## 2 g1 somker
## 3 g1 somker
## 4 g2 somker
## 5 g2 somker
## 6 g2 somker
网友评论