R之四格表与列联表数据的统计处理

四格表与列联表数据的统计处理

卡方检验用途：

• 推断两个或两个以上总体率或总体构成比间有无差别；
• 推断两种属性或两个变量(指标)之间有无关联性；
• 频数分布的拟合优度检验（检验某实际分布与理论分布的吻合情况）

卡方检验类型：

• 普通四格表资料卡方检验
• 配对四格表资料卡方检验
• 行×列（R×C）表资料卡方检验
• 行×列（R×C）列联表卡方检验

四格表卡方检验应用条件：

• 当n≥40且所有T≥5,用普通X2检验
• 当n≥40，但1≤T<5时, 用校正的X2检验
• 当n<40 或 T<1时，用四格表确切概率法
• 若p ≈ α ，改用四格表确切概率法（直接计算概率法）

一、基础介绍

1、二项分布、Poisson分布

二项分布(binomial distribution)：临床试验结果时常只有相互对立的结果，比如：生或死；有效或无效；阳性或阴性等。这种资料也叫二分类资料。其阳性率用π表示，阴性率用（1－ π）表示，其对应分布是二项分布

二项分布资料的假设检验： u检验或卡方检验，二者是等价的

Poisson分布：当π很小而n很大时，二项分布逼近Poisson分布，因此， Poisson分布是二项分布的特例，

Poisson分布主要应用： 单位容积、面积、时间或人群内某事件的发生数。比如：单位体积水中细菌数，单位体积空气中粉尘数，单位时间内放射性物质放射出的质点数，单位空间中某些昆虫数，一定人群中某种恶性肿瘤发生或罕见非传染性疾病的患病数或死亡数

2、卡方值的计算方法（通用公式）

例如：

**研究目的：**判断A与B两药的有效率有无差别？

表 5-1 两种药物治疗脑血管病有效率比较

分组	有效	无效	合计	有效率 %
实验组A药	73（a）	9（b）	82	89.02 %
对照组B药	52（c）	22（d）	74	70.27 %
合 计	125	31	156	80.13 %

检验假设

H0: A、 B两药的有效率无差别？
H1: A、 B两药的有效率有差别？

计算理论频数

表 5-1 两种药物治疗脑血管病有效率比较

分组	有效	无效	合计	有效率 %
实验组A药	73（65.7）	9（16.3）	82	89.02 %
对照组B药	52（59.3）	22（14.7）	74	70.27 %
合 计	125	31	156	80.13 %

3、卡方检验基本思想

从卡方值计算公式可以看出，卡方检验是检验实际分布和理论分布的吻合程度。若H0假设成立，则实际分布(A)和理论分布(T)相差不大， 卡方值应较小；若H0假设不成立，则实际分布(A)和理论分布(T)相差较大，卡方值应较大。另外卡方值的大小尚与格子数（自由度）有关，格子数越多（自由度越大），卡方值越大。可以根据卡方分布原理，由卡方值确定P值，从而作出推论。

自由度： V=(行数-1)(列数-1)

4、解题思路

1. 判断基本数字abcd是否给出？
2. 总数是否大于40？（大于40，可考虑用卡方检验，小于40改用Fisher确切概率法）
3. 判断是否需要校正？（计算最小理论数，若最小理论数>5，不需要校正。否则，用校正公式）
4. 按假设检验基本步骤进行卡方检验（基本公式或专用公式均可）

二、样本率与总体率的比较

样本率与总体率比较，判断样本与总体有无差别，是统计推断内容，因计算的是率，因而介绍一下：

p <- 0.2
n_case <- 48
n <- 152
stde <- sqrt (p *(1- p)/n)
u <- (n_case/n-p)/stde
u
pvalue <- 2*(1-pnorm(u))
pvalue 

三、 两个样本率的比较

x1 <-23

x2 <-13
n1 <- 80
n2 <-85
p <- (x1+x2)/(n1+n2)
std<-sqrt(p*(1-p)*(1/n1+1/n2))
u <- (x1/n1-x2/n2)/std
pvalue <- 2*(1-pnorm(u))
pvalue 

四、普通四格表代码实现

1、案例1

Example13_2  <- read.table ("example13_2.csv", header=TRUE, sep=",")
attach(Example13_2)
mytable  <-  xtabs(~a + b)
library(gmodels)
CrossTable(a, b)
chisq.test(mytable)
detach (Example13_2)

2、案例2

Example13_3  <- read.table ("example13_3.csv", header=TRUE, sep=",")
attach(Example13_3)
mytable  <-  xtabs(~a + b)
library(gmodels)
CrossTable(a, b)
# CrossTable(a, b, expected=TRUE,fisher=TRUE)这样也行
chisq.test(mytable)
chisq.test(mytable,correct=False) #默认是校正卡方
detach (Example13_3)

五、配对四格表

配对资料的展示形式：

1、配对资料的几种形式

• 同一个病例/样本用两种不同方法处理的结果
• 同一个病例/样本用两种方法诊断或检查结果
• 同一个病例/样品治疗前后的疗效等指标比较
• 结成对子的两个对象分别接受不同处理的结果

2、案例

Example13_5  <- read.table ("example13_5.csv", header=TRUE, sep=",")
attach(Example13_5)
mytable  <-  xtabs(~a + b)
library(gmodels)
CrossTable(a, b)
mcnemar.test(mytable)

使用kappa值评估两种检验方法是否一致

# Kappa.test(mytable, conf.level=0.95) # wrong 需提前载入 fmsb 包
library(vcd)
K<-Kappa(mytable)
K
confint(K)
summary(K)
print(K, CI = TRUE,level=0.95)
print(K, CI = TRUE,level=0.99)
detach (Example13_5)

六、R X C 列表卡方

1、列表分类

R X C 表可以分为：双向无序、单向有序、双向有序属性相同和双向有序属性不同等4 类。

1. 双向无序R X C 表， 表中两个分类变量皆为无序分类变量: ①若研究目的为多个样本率(或构成比) 比较，可用行×列表资料的卡方检验: ②若研究目的为分析两个分类变量之间有无关联性以及关系的密切程度，可以用行×列表资料的卡方检验，以及Pearson 列联系数分析

2. 单向有序R X C 表有两种形式：

   ①RXC 表中的分组变量是有序的，而结局变量是无序的，此种单向有序RXC 表资料可以用行×列表资料的卡方检验进行分析

   ②RXC 表中的分组变量是无序的，而结局变量是有序的，此种单向有序RXC表资料直用秩和检验分析

3. 双向有序属性相同R X C 表，表中的两个分类变量皆为有序且属性相同。实际上是2 X 2配对设计的扩展，此时宜用一致性检验(或称Kappa 检验)

4. 双向有序属性不同R X C 表，表中两个分类变量皆为有序且属性不同。对于该类资料，需要分析两个有序分类变量间是否存在线性变化趋势， 可用有序分组资料线性趋势检验

2、列变量为有序变量的行均分检验（Cochran-Mantel-Haenszel）

治愈赋值为4，显效赋值为3，好转赋值为2，无效赋值为1

ABC三种药物分布别用123表示

Example13_7  <- read.table ("example13_7.csv", header=TRUE, sep=",")
attach(Example13_7)
mytable  <-  xtabs(~a + b)
library(gmodels)
CrossTable(a, b)
library(vcdExtra)
CMHtest(mytable)     
detach (Example13_7)

输出为
Cochran-Mantel-Haenszel Statistics for a by b 

                 AltHypothesis  Chisq Df
cor        Nonzero correlation 46.297  1
rmeans  Row mean scores differ 58.678  2 #读取这一行
cmeans  Col mean scores differ 48.838  3
general    General association 66.958  6
              Prob
cor     1.0160e-11
rmeans  1.8125e-13
cmeans  1.4122e-10
general 1.7167e-12

也可以用回归，等级回归模型：

# 用Ordinal logistic回归模型
Example13_7  <- read.table ("example13_7.csv", header=TRUE, sep=",")
attach(Example13_7)
Example13_7$x1  <- ifelse (a==1, 1, 0)#对于ABC三种药物设置哑变量
Example13_7$x2  <- ifelse (a==2, 1, 0)
Example13_7$x3  <- ifelse (a==3, 1, 0)
library(rms)
fit1 <- lrm(b~ x1 + x2 ,  data=Example13_7, model=FALSE, x=FALSE, y=FALSE) #以药物C为参照
fit1
coefficients(fit1)
exp(coefficients(fit1))
detach (Example13_7)
#对于累积比数因变量模型，平行性假设决定了每个自变量的OR值对于前g-1 个模型是相同的。
#例如，自变量xl 的OR=8.044 ，表示使用A 药物治愈的可能性是C药物的8.044 倍;
#也表示使用A 药物显效或治愈的可能性是C药的8.044倍;
#同时也表示使用A 药物至少好转的可能性是C药的8.044 倍。

备注：假设治愈与显效相比，显效与好转相比，好转与无效相比的风险是相同的，因而会有三个回归系数

结果输出为：

fit1
Logistic Regression Model
 
 lrm(formula = b ~ x1 + x2, data = Example13_7, model = FALSE, 
     x = FALSE, y = FALSE)
 
 
 Frequencies of Responses
 
   1   2   3   4 
  51 126  73  20 
 
                        Model Likelihood      Discrimination    Rank Discrim.    
                              Ratio Test             Indexes          Indexes    
 Obs           270    LR chi2      66.98      R2       0.241    C       0.686    
 max |deriv| 3e-07    d.f.             2      R2(2,270)0.214    Dxy     0.372    
                      Pr(> chi2) <0.0001    R2(2,235.3)0.241    gamma   0.536    
                                              Brier    0.174    tau-a   0.249    
 
      Coef    S.E.   Wald Z Pr(>|Z|)
 y>=2  0.9782 0.2248   4.35 <0.0001 
 y>=3 -1.5521 0.2432  -6.38 <0.0001 
 y>=4 -3.7483 0.3375 -11.11 <0.0001 
 x1    2.0850 0.3088   6.75 <0.0001 
 x2    0.0028 0.2955   0.01 0.9926  
 
> coefficients(fit1)
        y>=2         y>=3         y>=4 
 0.978150920 -1.552142851 -3.748324155 
          x1           x2 
 2.084976695  0.002753265 
> exp(coefficients(fit1))
      y>=2       y>=3       y>=4         x1 
2.65953400 0.21179365 0.02355719 8.04440400 
        x2 
1.00275706 

3、行列变量均为有序变量检验方法

3.1 双向有序属性不同R X C表

例13-8 为了研究晶状体混浊程度是否与年龄相关，将资料整理为表13 -6 的形式，试编写R程序，分析年龄与晶状体混浊程度的相关关系。

# 双向有序属性不同 秩相关
Example13_8  <- read.table ("example13_8.csv", header=TRUE, sep=",")
attach(Example13_8)
cor(Example13_8, method="spearman")
cor.test(a,  b, method="spearman")
detach (Example13_8)

还可进行线性趋势检验

# 双向有序属性不同 线性趋势检验
Example13_8  <- read.table ("example13_8.csv", header=TRUE, sep=",")
attach(Example13_8)
library(gmodels)
CrossTable(a, b)
mytable <-  xtabs(~a + b)
chisq.test(mytable)
fit <- lm(a~b)
summary(fit)
coefficients(fit)
confint(fit)
detach (Example13_8)

3.2 双向有序属性相同R X C表

例13-11 某学校学生的文化课成绩和体育课成绩整理如表13-7 所示，试对学生文化课和体育课成绩进行一致性检验

计算kappa值：

Example13_10  <- read.table ("example13_10.csv", header=TRUE, sep=",")
attach(Example13_10)
library(gmodels)
CrossTable(a, b)
mytable  <-  xtabs(~a + b)
mcnemar.test(mytable)#配对卡方，意义在于如优秀率与不及格率有差别
library(fmsb)
Kappa.test(mytable, conf.level=0.95)
detach (Example13_10)

4、分层行列表分析 Mantel-Haenszel 检验

例13-12 为研究心肌梗塞与近期使用避孕药之间的关系，在5所医院中采用病例-对照研究方法调查了234名心肌梗塞病人与1742 名对照者使用口服避孕药状况，资料见表13-8 。请在排除了研究医院影响后，分析使用口服避孕药与否对是否患心肌梗塞病的影响情况。

Example13_12  <- read.table ("example13_12.csv", header=TRUE, sep=",")
attach(Example13_12)
mytable  <-  xtabs(~drug + case + hos) #分层输出各家医院四格表频数
mytable
prop.table(mytable,3) #分层输出各家医院四格表百分数
addmargins(mytable) #分层输出各家医院四格表边际频数
mantelhaen.test(mytable)
detach (Example13_12)

6、线性趋势卡方检验

例13-14 为了研究晶状体混浊程度是否与年龄相关，将资料整理成表13-10 的形式，试编写趋势卡方检验的R程序，分析年龄与晶状体混浊程度的相关关系

Example13_13  <- read.table ("example13_13.csv", header=TRUE, sep=",")
attach(Example13_13)
library(gmodels)
CrossTable(a, b)
mytable  <-  xtabs(~a + b)
chisq.test(mytable)
fit  <-  lm(a~b)
summary(fit)
coefficients(fit)#回归系数
confint(fit)#回归系数95%CI
detach (Example13_13)

七、秩和检验处理四格表及列联表资料

1、等级分组资料的非参数检验

例14- 13 用某药治疗不同病情( A 型和B 型)的老年慢性支气管炎病人， 疗效见表14-10 ，试比较该药对两种病情的疗效

1为控制，2为显效，3为有效，4为治愈

example14_15  <- read.table ("example14_15.csv", header=TRUE, sep=",")
attach(example14_15)
wilcox.test(x ~ g)
detach(example14_15)

2、多个组的等级指标的非参数检验

例14-14 4 种疾病患者痰液内嗜酸性粒细胞的检查结果见表14-11。试分析4 种疾病患者痰液内嗜酸性粒细胞有无差别

# 多个组的等级指标的非参数检验
example14_16  <- read.table ("example14_16.csv", header=TRUE, sep=",")
attach(example14_16)
kruskal.test(x~ g)
library(nparcomp)
nparcomp(x ~ g, data=example14_16, alternative = "two.sided")
detach(example14_16)
cox.test(x ~ g)
detach(example14_15)