PH525x series - Association Test
2019-11-19 本文已影响0人
3between7
- Chi-square Test（卡方检验）
假设有250个个体，有些患有某种疾病：在基因型为aa的人群中，有20%的人患有疾病，而剩余的人中有10%患有疾病。那么如果另选250个个体，我们是否还能看到这样的结果？
为了解答这一问题，首先构建如下一组数据：
disease=factor(c(rep(0,180),rep(1,20),rep(0,40),rep(1,10)),
labels=c("control","cases"))
genotype=factor(c(rep("AA/Aa",200),rep("aa",50)),
levels=c("AA/Aa","aa"))
dat <- data.frame(disease, genotype)
dat <- dat[sample(nrow(dat)),] #打乱顺序
head(dat)
## disease genotype
## 67 control AA/Aa
## 93 control AA/Aa
## 143 control AA/Aa
## 225 control aa
## 50 control AA/Aa
## 221 control aa
## 构建2X2的表格
tab <- table(genotype,disease)
tab
## disease
## genotype control cases
## AA/Aa 180 20
## aa 40 10
概括上述结果的一个经典统计手段是计算odd ratio（OR）,即比值比。对于发病率很低的疾病来说，它是相对危险度的精确估计值。计算方法如下：
##就是先计算基因型分别为aa、AA/Aa的群体中，患病与未患病人数的比值，然后再计算两个比值的比值。
(tab[2,2]/tab[2,1]) / (tab[1,2]/tab[1,1])
## [1] 2.25
但是为了计算p值，我们并不会直接使用OR值。相反，我们首先会假设基因型与疾病之间无关联（即零假设），然后计算2x2表格中每一个cell的值。那么，在零假设下，患病的概率是：
p=mean(disease=="cases")
p
## [1] 0.12
那么零假设下，期望的table应该是：
#按照患病与未患病概率去计算相应cell的数值
expected <- rbind(c(1-p,p)*sum(genotype=="AA/Aa"),
c(1-p,p)*sum(genotype=="aa"))
dimnames(expected)<-dimnames(tab)
expected
## disease
## genotype control cases
## AA/Aa 176 24
## aa 44 6
而后通过比较理论table与实际table之间的偏差，即：
再通过查表便可以获得p值大小。当然也可以直接使用chisq.test函数获得p值：
chisq.test(tab)$p.value
## [1] 0.08857435
#p值大于0.05，所以无法拒绝原假设，无法得出基因型与患病与否之间有关联的结论
- "大样本，小p值"问题
在进行统计推断时，仅报道p值并不可取，很多遗传相关方面的研究都过度强调了p值。在这种研究中，往往样本量非常大，p值非常小，但是OR值往往是相当小的：很少有大于1的。以本例为例，我们将样本量扩大10倍，但不改变比值，重新计算p值后会发现，p值会变得非常小：
tab<-tab*10
chisq.test(tab)$p.value
## [1] 1.219624e-09
- OR值的置信区间
由于OR值是比值的比值，所以不能用诸如CLT等理论简单的计算置信区间，但是可以通过广义线性模型理论去计算（使用OR的log值）：
fit <- glm(disease~genotype,family="binomial",data=dat)
coeftab<- summary(fit)$coef
coeftab
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.1972246 0.2356828 -9.322803 1.133070e-20
## genotypeaa 0.8109302 0.4249074 1.908487 5.632834e-02
coeftab第二行提供了OR的log值的预期值与SE，由于预期值大致为正态分布，故可以进行置信区间的计算:
ci <- coeftab[2,1] + c(-2,2)*coeftab[2,2]
#取幂值
exp(ci)
## [1] 0.9618616 5.2632310
