卡方检验用途:
卡方检验类型:
四格表卡方检验应用条件:
二项分布(binomial distribution):临床试验结果时常只有相互对立的结果,比如:生或死;有效或无效;阳性或阴性等。这种资料也叫二分类资料。其阳性率用π表示,阴性率用(1- π)表示,其对应分布是二项分布
二项分布资料的假设检验: u检验或卡方检验,二者是等价的
Poisson分布:当π很小而n很大时,二项分布逼近Poisson分布,因此, Poisson分布是二项分布的特例,
Poisson分布主要应用: 单位容积、面积、时间或人群内某事件的发生数。比如:单位体积水中细菌数,单位体积空气中粉尘数,单位时间内放射性物质放射出的质点数,单位空间中某些昆虫数,一定人群中某种恶性肿瘤发生或罕见非传染性疾病的患病数或死亡数
例如:
**研究目的:**判断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 % |
从卡方值计算公式可以看出,卡方检验是检验实际分布和理论分布的吻合程度。若H0假设成立,则实际分布(A)和理论分布(T)相差不大, 卡方值应较小;若H0假设不成立,则实际分布(A)和理论分布(T)相差较大,卡方值应较大。另外卡方值的大小尚与格子数(自由度)有关,格子数越多(自由度越大),卡方值越大。可以根据卡方分布原理,由卡方值确定P值,从而作出推论。
自由度: V=(行数-1)(列数-1)
样本率与总体率比较,判断样本与总体有无差别,是统计推断内容,因计算的是率,因而介绍一下:
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
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)
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)
配对资料的展示形式:
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 表可以分为:双向无序、单向有序、双向有序属性相同和双向有序属性不同等4 类。
双向无序R X C 表, 表中两个分类变量皆为无序分类变量: ①若研究目的为多个样本率(或构成比) 比较,可用行×列表资料的卡方检验: ②若研究目的为分析两个分类变量之间有无关联性以及关系的密切程度,可以用行×列表资料的卡方检验,以及Pearson 列联系数分析
单向有序R X C 表有两种形式:
①RXC 表中的分组变量是有序的,而结局变量是无序的,此种单向有序RXC 表资料可以用行×列表资料的卡方检验进行分析
②RXC 表中的分组变量是无序的,而结局变量是有序的,此种单向有序RXC表资料直用秩和检验分析
双向有序属性相同R X C 表,表中的两个分类变量皆为有序且属性相同。实际上是2 X 2配对设计的扩展,此时宜用一致性检验(或称Kappa 检验)
双向有序属性不同R X C 表,表中两个分类变量皆为有序且属性不同。对于该类资料,需要分析两个有序分类变量间是否存在线性变化趋势, 可用有序分组资料线性趋势检验
治愈赋值为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
例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)
例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)
例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)
例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)
例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)
例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)
联系客服