R语言和医学统计学（3）：卡方检验

本文首发于公众号：医学和生信笔记，完美观看体验请至公众号查看本文。

医学和生信笔记，专注R语言在临床医学中的使用，R语言数据分析和可视化。

前言

这是R语言和医学统计学的第3篇内容。

主要是用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

结果也是和课本一致的！

本文首发于公众号：医学和生信笔记，完美观看体验请至公众号查看本文。

医学和生信笔记，专注R语言在临床医学中的使用，R语言数据分析和可视化。