原文链接：http://blog.sina.com.cn/s/blog_5de124240101rcn1.html

卡方分析用于比较不同组之间的构成比，它的零假设是假定各组之间的构成是相同的，计算出理论每组的理论构成比，再计算理论值与实际值的差别，如果差别大的话，就拒绝零假设。它的扩展分析方法有Fisher精确分析，Ridit分析，CMH分析。与卡方分析有关的数据形式主要有四格表检验，无序RC表的卡方检测，单向有序表的卡方检验，双向有序表的卡方检测。

    卡方分布是一种连续分布，只有一个参数，即自由度ν，当自由度ν≤2时,卡方分布的曲线呈L形，随着ν的增加，曲线逐渐趋于对称，当自由度无穷大时，接近于正态分布

# 绘制此图的R命令：

x <- c(seq(0,10,length = 1000)) #生成从0到10之间的1000个数字

y1 <- dchisq(x,1) # 生成对x的卡方密度值，自由度为2

y6 <- dchisq(x,6) # 生成对x的卡方密度值，自由度为3

y10 <- dchisq(x,10) # 生成对x的卡方密度值，自由度为5

plot(x,y1,type = "l", ylim = c(0,0.5),xlim=c(0,15),xlab="", ylab="", lty = 1, main = expression(paste(chi^2,"Distribution")))

# type ="l"表示绘制线，如果是n则不绘线；lty=1表示绘制实线

lines(x,y6, type = "l", xlab = "", ylab = "", lty =2) # lty=2表示绘制虚线

lines(x,y10, type = "l", xlab = "", ylab = "", lty= 3) # lty=3表示绘制点线

text(locator(3), c(expression(chi^2*(1)),expression(chi^2*(6)),expression(chi^2*(10))))

# locator()读取鼠标的坐标，expression_r()用以添加数学公式

 

 

目录：

普通四格表卡方检验
	
配对四格表卡方检验与Kappa检验
	
Fisher精确检验
	
无序RC表检验
	
卡方分析的两两比较
	
单向有序RC表检验（Ridit分析）
	
双向有序表的CMH检验
 

一、普通四格表卡方检验

卡方分布是一种连续型分布，而四格表资料是分类资料，属离散型分布，由此计算的卡方值的抽样分布也应当是不连续的，当样本量较小时，两者间的差异不可忽略，应进行连续性校正（在每个单元格的残差中都减去0.5）。

若n > 40，此时有1< T <5时，需计算Yates连续性校正（Yates' correction）卡方值；将原始公式分子部分变为（|Ai-Ti|-0.5）即可。在例7-1中，期望值没有在1< T <5这个区间的，因此不用进行Yates连续性校正，需要在参数中设为chisq.test(x,correct=F)；但 例7-2中需要进行连续性校正，R中默认的有，不用设定。

例7-1:某医院欲比较异梨醇口服液（试验组）和氢氯噻嗪+地塞米松（对照组）降低颅内压的疗效。将200例颅内压增高症患者随机分为两组，结果见表7-1。问两组降低颅内压的总体有效率有无差别？（《医学统计学》第三版，孙振球）


			
组别
			
			

			
有效
			
			

			
无效
			
		

			
试验组
			
			

			
99 
			
			

			
5 
			
		

			
对照组
			
			

			
75 
			
			

			
21 
			
		
# 源代码例7-1：
data71<-matrix(c(99,75,5,21),nr = 2,dimnames = list(c("试验组", "对照组"),c("有效", "无效")))

data71

chi <- chisq.test(data71)

chi$expected #　显示理论数

chisq.test(data71,correct=F) # correct = F表示不进行Yates连续性校正

 

# 运行结果 
> data71<-matrix(c(99,75,5,21),nr = 2,dimnames = list(c("试验组", "对照组"),c("有效", "无效")))

> data71

       有效 无效

试验组   99    5

对照组   75   21

> chi <- chisq.test(data71)
> chi$expected #　显示理论数
        有效  无效

试验组 90.48 13.52

对照组 83.52 12.48
> chisq.test(data71,correct=F) # correct = F表示不进行Yates连续性校正

        Pearson's Chi-squared test

 data:  data71

X-squared = 12.8571, df = 1, p-value = 0.0003362

   

例7-2 某医师欲比较胞磷胆碱与神经节苷酯治疗脑血管疾病的疗效，将78例脑血管疾病患者随机分为两组，结果见表7-2。问两种药物治疗脑血管疾病的有效率是否有差别？（四格表卡方分析，《医学统计学》第三版，孙振球）


			
组别
			
			

			
有效
			
			

			
无效
			
		

			
胞磷胆碱组
			
			

			
46 
			
			

			
6 
			
		

			
神经节苷脂组
			
			

			
18 
			
			

			
8 
			
		
 

# 源代码例7-2：
data72<-matrix(c(46,18,6,8),nr = 2,dimnames = list(c("胞磷胆碱组", "神经节苷脂组"),c("有效", "无效")))

data72

chisq.test(data72) # 其中有理论数是小于5的，需要进行Yates校正，R自带这种校正

 

# 运行结果

> chisq.test(data72)

 

        Pearson's Chi-squared test with Yates' continuity correction

 

data:  data72

X-squared = 3.1448, df = 1, p-value = 0.07617
 

Warning message:
In chisq.test(compare_2) : Chi-squared approximation may be incorrect

> chisq.test(data72) # 查看理论数，里面有值小于5的，所以需要Yates校正，但R中默认有这种校正

                 有效     无效

胞磷胆碱组   42.66667 9.333333

神经节苷脂组 21.33333 4.666667

Warning message:

In chisq.test(data72) : Chi-squared approximation may be incorrect

   

二、配对四格表卡方检验与Kappa检验

配对四格表常用于两种检验方法、培养方法、诊断方法的比较。特点是对样本中各观察单位分别用两种方法处理，然后观察两种处理方法的某两分量的计数结果。这种资料可以进行两种分析，一种是分析A法和B法的阳性率是否有差异，另一种是分析A法和B法之间是否有关联(association)。对于第一种目的，如果b＋c>40，则用McNemar配对卡方检验;如果b＋c<40，则用校正的配对卡方检验。对于第二种目的，可进行一致性检验，即 Kappa检验(Kappa test)，这种检验是分析A法与B法在实现目的上，是否吻合。Kappa值的范围在-1至+1之间，取值越大，说明吻合程度越高。Kappa≥0.75一致性好；0.75>Kappa ≥0.4一致性一般；Kappa<0.4一致性较差。

 

例7-3 某实验室分别用乳胶凝集法和免疫荧光法对58名可疑系统红斑狼疮患者血清中抗核抗体进行测定，结果见表7-3。问两种方法的检测结果有无差别？（配对四格表的卡方分析，《医学统计学》第三版，孙振球）

 
			

			
乳胶凝集法
			
		

			
免疫荧光法
			
			

			
+ 
			
			

			
- 
			
		

			
+ 
			
			

			
11 
			
			

			
12 
			
		

			
- 
			
			

			
2 
			
			

			
33 
			
		
对于配对的四格表有四种情况：(1)两种检测方法都为阳性(a)；(2)两种方法都为阴性(d)；(3)免疫荧光法为阳性，乳胶凝集法为阴性(b)；(4)乳胶凝集法为阳性，免疫荧光法为阴性数(c)。

# 源代码 例7-3：
data73<-matrix(c(11,2,12,33),nr = 2,dimnames = list(c("+", "-"),c("+", "-")))

install.packages("epiR")

library(epiR) # 载入epiR包，使用其中的epi.kappa函数

mcnemar.test(data73) # 配对检验

epi.kappa(data73) # 一致性检验，即考察这两个方法在检验方面是否吻合

 

# 运行结果

> mcnemar.test(data73)

 

        McNemar's Chi-squared test with continuity correction

 

data:  data73
McNemar's chi-squared = 5.7857, df = 1, p-value = 0.01616

> epi.kappa(data73)

$kappa

est se lower upper

1 0.4550336 0.1209499 0.2179761 0.6920911

 

$z

test.statistic p.value

1 3.762165 0.000168449

结果分析：McNemar配对卡方检验的p值为0.01616,在0.05水平上，有显著性差异，免疫荧光法的总体阳性检出率要高于乳胶凝集法。

Kappa检验结果的p值小于0.01，估计值为0.455，说明两种方法的吻合度有统计学意义，但吻合度一般。

简单来说，McNemar检验是比较两种方法A与B对同一个对象的检验是否有差异，判断哪种方法比较好，比较容易检验出来结果。而Kappa检验就是比较A与B这两种方法检验出来的结果是不是一致的。在有些配对案例中可能会出现A的方法优于B，但A与B检验的结果不一致的情况。

 

三、精确Fisher检验

当四格表资料中出现n<40或T<1，或用一般公式计算出的Chi值所得到的概率与alpha接近时，需要用四格表资料的Fisher确切概率法，即对所有可能出现的情况进行排列组合，计算每一种情况的概率。

例7-4-某医师为研究乙肝免疫球蛋白预防胎儿宫内感染HBV的效果，将33例HBsAg阳性孕妇随机分为预防注射组和非预防组，结果见表7-4。问两组新生儿的HBV总体感染率有无差别？（精确Fisher检验，《医学统计学》第三版，孙振球）


			
组别
			
			

			
阴性
			
			

			
阳性
			
		

			
预防注射组
			
			

			
4 
			
			

			
18 
			
		

			
非预防组
			
			

			
5 
			
			

			
6 
			
		
 

# 源代码 例7-4： 
data74<-matrix(c(4,5,18,6),nr = 2,dimnames = list(c("预防注射组", "非预防组"),c("阳性", "阴性")))

chi<-chisq.test(data74)  

chi$expected # 查看理论频数       

fisher.test(data74)

 

# 运行结果

> data74<-matrix(c(4,5,18,6),nr = 2,dimnames = list(c("预防注射组", "非预防组"),c("阳性", "阴性")))
> chi<-chisq.test(data74) 

Warning message:
In chisq.test(data74) : Chi-squared approximation may be incorrect
> chi$expected # 查看理论频数      
           阳性 阴性

预防注射组    6   16

非预防组      3    8
> fisher.test(data74)

 

        Fisher's Exact Test for Count Data

 

data:  data74

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

 

例7-5-某单位研究胆囊腺癌、腺瘤的P53基因表达，对同期手术切除的胆囊腺癌、腺瘤标本各10份，用免疫组化法检测P53基因，资料见表7-6。问胆囊腺癌和胆囊腺瘤的P53基因表达阳性率有无差别？（精确Fisher检验，《医学统计学》第三版，孙振球）


			
组别
			
			

			
阴性
			
			

			
阳性
			
		

			
胆囊腺癌
			
			

			
6 
			
			

			
4 
			
		

			
胆囊腺瘤
			
			

			
1 
			
			

			
9 
			
		
 

# 源代码例7-5
data75<-matrix(c(6,1,4,9),nr = 2,dimnames = list(c("胆囊腺癌", "胆囊腺瘤"),c("阳性", "阴性")))

fisher.test(data75)

 

# 运行结果：

> fisher.test(data75)

         Fisher's Exact Test for Count Data

 data:  compare_5

p-value = 0.05728

alternative hypothesis: true odds ratio is not equal to 1

95 percent confidence interval:

   0.9487882 684.4235629

sample estimates:

odds ratio

   11.6367

   

四、无序R×C列表的卡方检验

例7-6-某医师研究物理疗法、药物治疗和外用膏药三种疗法治疗周围性面神经麻痹的疗效，资料见表7-8。问三种疗法的有效率有无差别（RC都无序，《医学统计学》第三版，孙振球）？


			
疗法
			
			

			
有效
			
			

			
无效
			
		

			
物理疗法
			
			

			
199 
			
			

			
7 
			
		

			
药物治疗组
			
			

			
164 
			
			

			
18 
			
		

			
外用膏药组
			
			

			
118 
			
			

			
26 
			
		
 

# 源代码：
data76<-matrix(c(199,164,118,7,18,26),nr = 3,dimnames = list(c("物理疗法组", "药物治疗组","外用膏药组"),c("有效", "无效")))

chi.result <- chisq.test(data76)

chi.result

chi.result$expected

 

 # 运行结果：
> data76<-matrix(c(199,164,118,7,18,26),nr = 3,dimnames = list(c("物理疗法组", "药物治疗组","外用膏药组"),c("有效", "无效")))
> chi.result <- chisq.test(data76)
> chi.result

         Pearson's Chi-squared test

 data:  data76

X-squared = 21.0377, df = 2, p-value = 2.702e-05

 

> chi.result$expected

               有效     无效

物理疗法组 186.2519 19.74812

药物治疗组 164.5526 17.44737

外用膏药组 130.1955 13.80451

 

例7-7-某医师在研究血管紧张素I转化酶(ACE)基因I/D多态与2型糖尿病肾病(DN)的关系时，将249例2型糖尿病患者按有无糖尿病肾病分为两组，资料见表7-9，问两组2型糖尿病患者的ACE基因型总体分布有无差别（RC都无序，《医学统计学》第三版，孙振球）？


			
组别
			
			

			
DD 
			
			

			
ID 
			
			

			
II 
			
		

			
DN组
			
			

			
42 
			
			

			
48 
			
			

			
21 
			
		

			
无DN组
			
			

			
30 
			
			

			
72 
			
			

			
36 
			
		
 

# 源代码7-7：
data77 <- matrix(c(42,30,48,72,21,36),nrow=2,dimnames=list(c("DN组","无DN组"),c("DD","ID","II")))

chi.result77 <- chisq.test(data77)

chi.result77$expected

chi.result77

 

# 运行结果：
> data77 <- matrix(c(42,30,48,72,21,36),nrow=2,dimnames=list(c("DN组","无DN组"),c("DD","ID","II")))
> chi.result77 <- chisq.test(data77)
> chi.result77$expected

             DD       ID       II

DN组   32.09639 53.49398 25.40964

无DN组 39.90361 66.50602 31.59036
> chi.result77

 

        Pearson's Chi-squared test

 

data:  data77

X-squared = 7.9127, df = 2, p-value = 0.01913

 

例7-8：测得某地5801人的ABO血型和MN血型结果如表7-10，问两种血型系统之间是否有关联？（RC都无序，《医学统计学》第三版，孙振球）


			
ABO血型
			
			

			
MN血型
			
		

			
M 
			
			

			
N 
			
			

			
MN 
			
		

			
O 
			
			

			
431 
			
			

			
490 
			
			

			
902 
			
		

			
A 
			
			

			
388 
			
			

			
410 
			
			

			
800 
			
		

			
B 
			
			

			
495 
			
			

			
587 
			
			

			
950 
			
		

			
AB 
			
			

			
137 
			
			

			
179 
			
			

			
32 
			
		
 

# 源代码例7-8：
data78 <- matrix(c(431,388,495,137,490,410,587,179,902,800,950,32),nrow=4,dimnames=list(c("O","A","B","AB"),c("M","N","MN")))

chi.result <- chisq.test(data78)

chi.result

sqrt(chi.result$statistic/(sum(data78)+chi.result$statistic)) # 列联系数

 

# 运行结果：
> data78 <- matrix(c(431,388,495,137,490,410,587,179,902,800,950,32),nrow=4,dimnames=list(c("O","A","B","AB"),c("M","N","MN")))
> chi.result <- chisq.test(data78)
> chi.result

 

        Pearson's Chi-squared test

 

data:  data78

X-squared = 213.1616, df = 6, p-value < 2.2e-16

 

> sqrt(chi.result$statistic/(sum(data78)+chi.result$statistic)) # 列联系数
X-squared

0.1882638 

 

五、卡方分析的两两比较：在对多个样本的构成比之间进行两两比较时，需要对检验水平进行校正。校正方法如下：

(1)多组间的两两比较：alpha'=alpha/N，N=n(n-1)/2 n为参加检验的组数，alpha为显著性水平，如例7-9。

(2)多个实验组与统一对照组的比较：alpha=alpha/(K-1)，K为实验组与对照组的组数之和，如例7-10。

 

例7-9  对例7-6中的资料进行两两比较，以推断是否任两种疗法治疗周围性面神经麻痹的有效率均有差别？

#源代码7-9

data76<-matrix(c(199,164,118,7,18,26),nr = 3,dimnames = list(c("物理疗法组", "药物治疗组","外用膏药组"),c("有效", "无效")))

chisq.test(data76[1:2,],correct=F) #比较物理疗法组与药物治疗组，R中默认四格表进行了Yates校正，此处设为F，不进行校正

chisq.test(data76[2:3,], correct=F) #比较药物治疗组与外用膏药组

chisq.test(data76[-2,],correct=F) # 比较物理疗法组与外用膏药组

 

> chisq.test(data76[1:2,],correct=F) #比较物理疗法组与药物治疗组

 

    Pearson's Chi-squared test

 

data: data76[1:2, ]

X-squared = 6.756, df = 1, p-value = 0.009343

 

> chisq.test(data76[2:3,], correct=F) #比较药物治疗组与外用膏药组

 

    Pearson's Chi-squared test

 

data: data76[2:3, ]

X-squared = 4.591, df = 1, p-value = 0.03214

 

> chisq.test(data76[-2,],correct=F) # 比较物理疗法组与外用膏药组

 

    Pearson's Chi-squared test

 

data: data76[-2, ]

X-squared = 21.3228, df = 1, p-value = 3.881e-06

由上述结果可以看出，在0.0125显著性水平上：此例中两两比较的次数为4,0.05/4=0.125，因此可以说：

物理疗法组与药物治疗组有差别

物治疗组与外用膏药组没有差别

物理疗法组与外用膏药组有差别

 

例 7-10  以表7-8资料中的药物治疗组为对照组，物理疗法组与外用膏药组为试验组，试分析两试验组与对照组的总体有效率有无差别？

# 源代码7-10：

chisq.test(data76[1:2,],correct=F) #比较物理疗法组与药物治疗组

chisq.test(data76[2:3,], correct=F) #比较药物治疗组与外用膏药组

 

> chisq.test(data76[1:2,],correct=F) #比较物理疗法组与药物治疗组

 

    Pearson's Chi-squared test

 

data: data76[1:2, ]

X-squared = 6.756, df = 1, p-value = 0.009343

 

>

> chisq.test(data76[2:3,], correct=F) #比较药物治疗组与外用膏药组

 

    Pearson's Chi-squared test

 

data: data76[2:3, ]

X-squared = 4.591, df = 1, p-value = 0.03214

此例中有3组，另外两组分别与药物治疗组比。2×(k-1)=2×(3-1)=4,0.05/4=0.125，虽然与上题的检验水平数值一样，但计算的过程是不同的。

物理疗法组与药物治疗组的p值为0.009343<0.0125，有显著性差异；外用膏药组与药物治疗组相比，p值为0.03214>0.0125，没有显著性差异。

 

六、单向有序RC表检验（Ridit分析）

Ridit分析全称为relative to an identified distribution unit，即按等级分布的频率转化为与某特定分布相对应的单位进行分析。属于半定量的分析方法，适合等级资料或数量分组不很确切的资料。常用于行有序或者是列有序的数据分析，属于卡方分析的一种扩展。但从查到的文献来看，这种分析方法，多数是中国人采用，国外采用的不多。

 

例6-1：比较三个药物的疗效，如下所示，在疗效一列中，治愈，显效，好转，无效这些疗效都是有顺序的，而药物A、B、C则是无序的，可以用Ridit分析。


			
疗效
			
			

			
药物
			
		

			
A 
			
			

			
B 
			
			

			
C 
			
		

			
治愈
			
			

			
17 
			
			

			
5 
			
			

			
3 
			
		

			
显效
			
			

			
51 
			
			

			
11 
			
			

			
17 
			
		

			
好转
			
			

			
33 
			
			

			
52 
			
			

			
47 
			
		

			
无效
			
			

			
7 
			
			

			
24 
			
			

			
26 
			
		
 

if(!require("Ridit"){

install.packages("Ridit")

library(Ridit)

} # 载入Ridit包

 

data1 <- c(17,51,33,7,5,11,52,24,3,17,47,26)

data1 <- matrix(data1,nrow=4,dimnames=list(c("治愈","显效","好转","无效"),c("A","B","C")))

ridit(data1,2)

# 参数2表示data1的分组信息在列，即药物A，药物B，药物C，有序的是行，即治愈，显效，好转，无效

 

# 运行结果：

> ridit(data1,2)

Ridit Analysis:

Group Label Mean Ridit

----- ----- ----------

1 A 0.3465

2 B 0.5921

3 C 0.5872

Reference: Total of all groups

chi-squared = 54.7507, df = 2, p-value = 1.291e-12

 

分析结果：卡方值为54.7507,p值小于0.01，表示药物的疗效不同，MeanRidit的意思为数值越大，疗效越差。但这个结果没有给出两两之间的比较。有人有SPSS进行Ridit分析，在进行两两比较的时候，是将Ridit转化为ANOVA进行比较。

 

例子6-2：甲中药与乙中药疗效的表格（来自于百度文库的课件）。


			
组别
			
			

			
无效
			
			

			
好转
			
			

			
显效
			
			

			
控制
			
		

			
甲中药
			
			

			
800 
			
			

			
1920 
			
			

			
680 
			
			

			
60 
			
		

			
乙中药
			
			

			
10 
			
			

			
60 
			
			

			
26 
			
			

			
18 
			
		
 

# 源代码6-2：

data2 <- c(800,10,1920,60,680,26,60,18)

data2 <- matrix(data2, nrow=2,dimnames=list(c("甲中药","乙中药"),c("无效","好转","显效","控制")));data2

ridit(data2,1)

 

# 运行结果：

> data2 <- matrix(data2, nrow=2,dimnames=list(c("甲中药","乙中药"),c("无效","好转","显效","控制")));data2

无效 好转 显效 控制

甲中药 800 1920 680 60

乙中药 10 60 26 18

> ridit(data2,1)

 

Ridit Analysis:

 

Group Label Mean Ridit

----- ----- ----------

1 甲中药 0.4957

2 乙中药 0.6318

 

Reference: Total of all groups

chi-squared = 30.2509, df = 1, p-value = 3.796e-08

 

 

例6-3：用中药复方治疗3330支气管炎病人，同时考察了另一种中药制剂治疗94例病人。结果如下

 


			
组别
			
			

			
处方1
			
			

			
处方2
			
			

			
处方3
			
			

			
处方4
			
		

			
痊愈
			
			

			
45 
			
			

			
11 
			
			

			
11 
			
			

			
4 
			
		

			
显著好转
			
			

			
38 
			
			

			
24 
			
			

			
8 
			
			

			
8 
			
		

			
好转
			
			

			
51 
			
			

			
34 
			
			

			
16 
			
			

			
25 
			
		

			
无效
			
			

			
38 
			
			

			
31 
			
			

			
18 
			
			

			
25 
			
		
 

data3 <- c(45,38,51,38,11,24,34,31,11,8,16,18,4,8,25,25)

data3 <- matrix(data3, nrow=4,dimnames=list(c("痊愈","显著好转","好转","无效"),c("处方1","处方2","处方3","处方4")));data3

ridit(data3,2)

 

> data3 <- matrix(data3, nrow=4,dimnames=list(c("痊愈","显著好转","好转","无效"),c("处方1","处方2","处方3","处方4")));data3

处方1 处方2 处方3 处方4

痊愈 45 11 11 4

显著好转 38 24 8 8

好转 51 34 16 25

无效 38 31 18 25

> ridit(data3,2)

 

Ridit Analysis:

 

Group Label Mean Ridit

----- ----- ----------

1 处方1 0.4382

2 处方2 0.5297

3 处方3 0.5178

4 处方4 0.6084

 

Reference: Total of all groups

chi-squared = 19.243, df = 3, p-value = 0.0002435

 

七、双向有序表的CMH检验

RC有序表并不适合用卡方分析，如果列或行有序，则无论怎么调换顺序，得到的卡方值是一样的。因此在面临行或列有序分析的时候通常采用Ridit分析；在面临行与列都有序的时候则需要CMH(Cochran-Mantel-Haenszel)分析方法，例如研究年龄与冠状动脉硬化之间的关系等。《医学统计学》（第三版，孙振球）中没有提到RC有序表中的统计方法，只是把算法列了出来，查过之后，发现比较符合CMH检验。CMH检验在新药中多用于多中心试验的检验。在书中，只用到了它对于两个有序变量的线性相关分析。在R中可以使用vcdExtra包中的CMHtest函数进行这种分析。

 

CMHtest的参数如下：

CMHtest(x, strata = NULL,rscores = 1:R, cscores = 1:C,types = c("cor", "cmeans", "rmeans", "general"), overall=FALSE, details=overall, ...)

 

x:两列或两列以上的列联表，数据为数组(array)或表(table)，可以用分类标签标注行与列。

formula:用于指定数据列联表中变量的公式。可以单边公式(array数据)，也可以是双边公式(data.frame)。可以为分层数据指定一个条件公式。

data:数据框，或table，或ftable。

strata:三维或高维表，其名称与因子的数目可以作为一个分层依据。默认情况下，前两个因子作为主要的表变量，其它的视为分层因子。

rscores：行得分。可以为一组数字(通常是1:R)，或者是用字符串"midrank"来进行标准中秩得分，设置为NULL则对行不进行得分统计。

cscores ：同rscores.

types：CMH的计算数据：一个或多个c("cor", "cmeans", "rmeans", "general")检验，或者用"ALL"指全部检验。

 

例7-11 某研究者欲研究年龄与冠状动脉粥样硬化等级之间的关系，将278例尸解资料整理成表7-13，问年龄与冠状动脉粥样硬化等级之间是否存在线性变化趋势？


			
表7-13 年龄与冠状动脉硬化的关系
			
		

			
年龄(岁)

			
(X) 
			
			

			
冠状动脉硬化等级(Y)
			
			

			
合计
			
		

			
— 
			
			

			
+ 
			
			

			
++ 
			
			

			
+++ 
			
		

			
20～
			
			

			
70 
			
			

			
22 
			
			

			
4 
			
			

			
2 
			
			

			
98 
			
		

			
30～
			
			

			
27 
			
			

			
24 
			
			

			
9 
			
			

			
3 
			
			

			
63 
			
		

			
40～
			
			

			
16 
			
			

			
23 
			
			

			
13 
			
			

			
7 
			
			

			
59 
			
		

			
≥50 
			
			

			
9 
			
			

			
20 
			
			

			
15 
			
			

			
14 
			
			

			
58 
			
		

			
合计
			
			

			
122 
			
			

			
89 
			
			

			
41 
			
			

			
26 
			
			

			
278 
			
		
# 源代码

if(!require("vcdExtra "){

install.packages("vcdExtra ")

library(vcdExtra)

} # 载入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) # 分析数据

 

# 运行结果
> CMHtest(data711)

Cochran-Mantel-Haenszel Statistics for age by degree


			
   
			
			

			
AltHypothesis
			
			

			
Chisq
			
			

			
Df
			
			

			
Pro
			
		

			
Cor
			
			

			
Nonzero correlation
			
			

			
63.389
			
			

			
1
			
			

			
1.6962e-15
			
		

			
Cmeans
			
			

			
Col mean scores differ
			
			

			
63.450
			
			

			
3
			
			

			
1.0758e-13
			
		

			
Rmeans
			
			

			
Row mean scores differ
			
			

			
65.669
			
			

			
3
			
			

			
3.6072e-14
			
		

			
general
			
			

			
General association
			
			

			
71.176
			
			

			
9
			
			

			
8.9519e-12
			
		

> chisq.test(data711)

 

        Pearson's Chi-squared test

 

data:  data711

X-squared = 71.4325, df = 9, p-value = 7.97e-12

注：卡方检验的结果值为71.4325，表明年龄与冠状动脉粥样硬化之间有关系。CMH检验的结果中有三项，主要是看Nonzero correlation这一项，它表示，列与行有序的时候，数据的检验结果，它的H0是年龄与冠状动脉粥样硬化之间没有线性关系，它的p值小于0.01，因此可以判断年龄与冠状动脉粥样硬化等级之间还存在线线性趋势。

 

八、拟合优度与卡方检验

例7-12 观察某克山病患者的空间分布情况，调查者将该地区划分为279个取样单位，统计各取样单位历年累计病例数，问是否服从Poisson分布

# 源代码：
x<-0:6
y<-c(26,51,75,63,38,17,9)

mean<-mean(rep(x,y))

q<-ppois(x,mean)

n<-length(y)

p<-c()

p[1]<-q[1]

p[n]<-1-q[n-1]
for(i in 2:(n-1))

p[i]<-q[i]-q[i-1]

chisq.test(y, p=p,correct=F)

 

# 运行结果：
> chisq.test(y, p=p,correct=F)

 

        Chi-squared test for given probabilities

 

data:  y

X-squared = 2.0569, df = 6, p-value = 0.9144

统计学简介之十四——列联表介绍

一、列联表构造






二、卡方检验

数据来源： 

图灵社区 
http://www.ituring.com.cn/book/1809 

1、创建列联表



dat <- read.csv("sample.csv")
 table(dat)


管道处理，用左边的数据生成表

install.packages("dplyr")
library(dplyr)
 dat  %>% table //等价于 table(dat)

独立性检验



dat2 <- dat  %>% table
chisq.test(dat2)
结果是
Pearson's Chi-squared test with Yates' continuity correction

data:  dat2
X-squared = 0.15476, df = 1, p-value = 0.694


皮尔逊卡方检验的结果显示，p大于0.05.故不是小概率事件，不拒绝原假设，从而两因素是相互独立的。

用python和Excel进行列联表分析
卡方检验回顾（独立性检验）
行(\(r_i\))列(\(c_j\))
\(j=1\)
\(j=2\)
\(\cdots\)
合计
\(i=1\)
\(f_{11}\)
\(f_{12}\)
\(r_1\)
\(i=2\)
\(f_{21}\)
\(f_{22}\)
\(r_2\)
\(\vdots\)
\(\vdots\)
\(\vdots\)
\(\vdots\)
合计
\(c_1\)
\(c_2\)
\(\cdots\)
\(n\)
建立原假设和备择假设
\[H_0:\text{两变量独立 }vs\text{ }H_1:\text{两变量不独立 }
\]
计算自由度和理论频数
自由度：\((r-1)*(c-1)\)，其中\(r,c\)分别是行数和列数
我们将实际频数记作\(f_{ij}\)，期望（理论）频数记作\(e_{ij}\),且\(e_{ij}=\frac{r_i*c_j}{n}\)
计算卡方统计量
\[\chi^{2}=\sum\limits_{i=1}^{r} \sum\limits_{j=1}^{c} \frac{\left(f_{i j}-\mathrm{e}_{i j}\right)^{2}}{\mathrm{e}_{i j}} \sim \chi^{2}(d f)
\]
实际观察次数与理论次数之差的平方再除以理论次数得到的统计量近似服从卡方分布
查χ2方分布临界值表，确定接受域
\[\chi_{1-\frac{\alpha}{2}}^{2}(\mathrm{df})<\chi^{2}<\chi_{\frac{\alpha}{2}}^{2}(\mathrm{df})
\]
数据示例
真实频数
男
女
总计
未流失
29699
8954
38653
流失
9660
2762
12422
总计
39359
11716
51075
Excel操作
计算期望频数
期望频数
男
女
总计
未流失
29786.46
9572.54
39359
流失
8866.54
2849.46
11716
总计
38653
12422
51075
以及真实频数和期望频数的差值
0.257
0.799
0.863
2.684
p值：excel用=CHISQ.TEST(真实频数,期望频数)，wps用=CHITEST(真实频数,期望频数)来计算。
卡方统计量：用=SUM()合计差值区域即可
自由度：如前文所述
1-置信水平：一般选0.1,0.05,0.01
接受域上下界：excel用=CHISQ.INV.RT(置信水平/2,自由度)，wps用=CHISQ.INV(置信水平/2,自由度)
所给数据计算结果如下：
p值
0.031916
卡方统计量
4.603016363
自由度
1
\(\alpha\)
0.01
上界
3.92704E-05
下界
6.634896601
Python代码解决方案
使用scipy.stats.chi2_contingency来进行卡方检验。
import numpy as np
from scipy.stats import chi2_contingency
#性别
d = np.array([[29699,9660], [8954,2762]])
chi2_contingency(d, correction = False)
#默认correction = True，这里是为了和excel计算结果一致
关于correction参数，如果为True，且自由度为1，则应用Yates校正以保持连续性。 校正的效果是将每个观察值向相应的期望值调整0.5。
Yates correction资料点击此处，总之它会使得卡方统计量更接近卡方分布。
输出结果如下
(4.60301636324016,
0.03191575584773069,
1,
array([[29786.45965737,  8866.54034263],
[ 9572.54034263,  2849.45965737]]))
输出结果分别为
卡方统计量
p值
自由度
期望频数
这题p值约为0.032，也就是说当显著性水平\(\alpha\)大于0.032时我们会拒绝原假设，当\(\alpha<0.032\)时我们会接受原假设。如果选择\(\alpha=0.05\)，那么拒绝原假设，认为性别与流失是不独立的。