R语言和医学统计学系列（3）：卡方检验
前言
这是R语言和医学统计学的第3篇内容。
第1篇请见：R语言和医学统计学系列（1）：t检验
第2篇请见：R语言和医学统计学系列（2）：方差分析
使用R语言进行统计学是我学习R语言最开始的初衷，没想到从此一发不可收拾，打开了新世界的大门。这个系列也是我最开始学习R语言时的笔记。希望对大家有帮助。
主要是用R语言复现课本中的例子。我使用的课本是孙振球主编的《医学统计学》第4版，封面如下：
我的研究生课程并没有把整本书的全部学完，只学习了其中的一部分，因此本系列也是只针对其中学过的部分进行复现。另外对于统计描述部分也不在这里探讨。
四格表资料的卡方检验
使用课本例7-1的数据。
首先是构造数据，本次数据自己从书上摘录。。
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)
str(data1)
## 'data.frame': 200 obs. of 3 variables:
## $ ID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ treat: Factor w/ 2 levels "placebo","treated": 2 2 2 2 2 2 2 2 2 2 ...
## $ impro: Factor w/ 2 levels "marked","none": 1 1 1 1 1 1 1 1 1 1 ...
数据一共3列，第一列是id，第二列是治疗方法，第三列是等级（有效和无效）。
简单看下各组情况：
table(data1$treat,data1$impro)
##
## marked none
## placebo 75 21
## treated 99 5
做卡方检验有2种方法，分别演示：
方法1
直接使用gmodels里面的CrossTable()函数，非常强大，直接给出所有结果，和SPSS差不多。
library(gmodels)
CrossTable(data1$treat, data1$impro, digits = 4, expected = T, chisq = T, fisher = T, mcnemar = T, format = "SPSS")
##
## Cell Contents
## |-------------------------|
## | Count |
## | Expected Values |
## | Chi-square contribution |
## | Row Percent |
## | Column Percent |
## | Total Percent |
## |-------------------------|
##
## Total Observations in Table: 200
##
## | data1$impro
## data1$treat | marked | none | Row Total |
## -------------|-----------|-----------|-----------|
## placebo | 75 | 21 | 96 |
## | 83.5200 | 12.4800 | |
## | 0.8691 | 5.8165 | |
## | 78.1250% | 21.8750% | 48.0000% |
## | 43.1034% | 80.7692% | |
## | 37.5000% | 10.5000% | |
## -------------|-----------|-----------|-----------|
## treated | 99 | 5 | 104 |
## | 90.4800 | 13.5200 | |
## | 0.8023 | 5.3691 | |
## | 95.1923% | 4.8077% | 52.0000% |
## | 56.8966% | 19.2308% | |
## | 49.5000% | 2.5000% | |
## -------------|-----------|-----------|-----------|
## Column Total | 174 | 26 | 200 |
## | 87.0000% | 13.0000% | |
## -------------|-----------|-----------|-----------|
##
##
## Statistics for All Table Factors
##
##
## Pearson's Chi-squared test
## ------------------------------------------------------------
## Chi^2 = 12.85707 d.f. = 1 p = 0.0003362066
##
## Pearson's Chi-squared test with Yates' continuity correction
## ------------------------------------------------------------
## Chi^2 = 11.3923 d.f. = 1 p = 0.0007374901
##
##
## McNemar's Chi-squared test
## ------------------------------------------------------------
## Chi^2 = 50.7 d.f. = 1 p = 1.076196e-12
##
## McNemar's Chi-squared test with continuity correction
## ------------------------------------------------------------
## Chi^2 = 49.40833 d.f. = 1 p = 2.078608e-12
##
##
## Fisher's Exact Test for Count Data
## ------------------------------------------------------------
## Sample estimate odds ratio: 0.1818332
##
## Alternative hypothesis: true odds ratio is not equal to 1
## p = 0.0005286933
## 95% confidence interval: 0.05117986 0.5256375
##
## Alternative hypothesis: true odds ratio is less than 1
## p = 0.0002823226
## 95% confidence interval: 0 0.4569031
##
## Alternative hypothesis: true odds ratio is greater than 1
## p = 0.9999541
## 95% confidence interval: 0.06281418 Inf
##
##
##
## Minimum expected frequency: 12.48
#可以通过设定参数直接给出理论频数，卡方检验，校正的卡方检验及Fisher精确概率法的结果
可以看到这个函数直接给出所有结果，根据需要自己选择合适的即可。
本例符合pearson卡方，卡方值为12.85707，p<0.01，和课本一致。
方法2
先把数据变成2x2列联表，然后用chisq.test函数做
mytable <- table(data1$treat,data1$impro)
mytable
##
## marked none
## placebo 75 21
## treated 99 5
chisq.test(mytable,correct = F) # 和SPSS一样
##
## Pearson's Chi-squared test
##
## data: mytable
## X-squared = 12.857, df = 1, p-value = 0.0003362
这个结果和课本也是一致的，和SPSS算出来的也是一样的。
“2种方法都演示过了，很明显方法1更简单直接，无需考虑各种条件，根据结果自己选择即可！
关于四格表资料卡方检验的专用公式/四格表资料卡方检验的校正公式/配对四格表资料的卡方检验/四格表资料的Fisher精确概率法，只要用方法1可直接解决，就不在继续演示了，大家可以自己尝试！
连续校正的卡方检验也是使用chisq.test()函数，需加上correct = T参数；Fisher精确概率法，也可使用fisher.test()函数完成。
行 x 列表资料的卡方检验
行 x 列表资料的卡方检验有很多种情况，不是所有的列联表资料都可以直接用卡方检验，大家要注意甄别！
使用课本例7-6的数据。
首先是构造数据，本次数据直接读取，也可以自己手动摘录。
df <- foreign::read.spss("E:/各科资料/医学统计学/研究生课程/5计数资料统计分析18-9研/例07-06物理药物外用膏用.sav", to.data.frame = T)
str(df)
## 'data.frame': 6 obs. of 3 variables:
## $ 疗法: Factor w/ 3 levels "物理疗法组","药物治疗组",..: 1 1 2 2 3 3
## $ 疗效: Factor w/ 2 levels "有效","无效": 1 2 1 2 1 2
## $ f : num 199 7 164 18 118 26
数据一共3列，第1列是疗法，第2列是有效无效，第3列是频数.
进行 行 x 列表资料的卡方检验，首先要对数据格式转换一下，变成table或者矩阵：
M <- as.table(rbind(c(199, 7),
c(164, 18),
c(118, 26)))
dimnames(M) <- list(trt = c("物理", "药物", "外用"),
effect = c("有效","无效"))
进行 行 x 列表资料的卡方检验：
kf <- chisq.test(M, correct = F)
kf
##
## Pearson's Chi-squared test
##
## data: M
## X-squared = 21.038, df = 2, p-value = 2.702e-05
结果也是和课本一致的！
欢迎大家交流，可直接评论区留言或添加我的微信。
欢迎关注我的公众号：医学和生信笔记
“医学和生信笔记 公众号主要分享：1.医学小知识、肛肠科小知识；2.R语言和Python相关的数据分析、可视化、机器学习等；3.生物信息学学习资料和自己的学习笔记！
往期精彩内容：
R语言画多时间点ROC和多指标ROC曲线
R语言ggplot2画相关性热图
R语言画好看的聚类树