title: DALS007-统计推断(Inference)06-关联检验(Association Tests)
date: 2019-08-07 12:0:00
type: "tags"
tags:
- 统计推断
categories: - 生物统计
前言
前面所提到的案例都是普通的统计学案例，这些普通案例指的是，它们的数据是二元的(binary)，分类的(categorical)或有序的(ordinal)。现在我们考虑一种特殊的案例，例如基因型数据，现在我们有两组基因型数据，一组基因是AA/Aa，另外一组是aa（某个特定疾病的对照）。现在我们的问题是，这个基因型(genotype)是否与疾病相关。这个案例与我们前面提到的案例就有点像了，我们有两个总体，即AA/Aa与aa，我们使用1或0来表示它们对应的疾病状态，其中0是不患病，1是患病的。此时我们不能使用t检验，因为这种数据明显不符合正态分布，无法使用t检验。如果样本足够大的话，可以使用CLT，否则，可以使用关联检验(association test)。
女士品茶
关于假设检验最有名的一个案例就是Fisher提出的。这个案例大概是这个样子的，一名女士声称，如果牛奶与红茶加入的先后顺序不同，得到的奶茶口味也不同，并且她自己能品尝出来。Fisher就给了4杯奶茶，这4杯奶茶加入的花与牛奶顺序是不同的，是随机的。后来这个女士品尝对了3个，此时，我们能认为这个女士确实有这种能力么，她品尝出的这3杯奶茶是否是出于偶然因素？假设检验就能回答这个问题，并且能够对这种偶然因素进行定量。
现在我们提出我们的基础问题：假如这个女士品尝出的这3杯奶茶是她自己瞎猜（瞎猜说明她没有品尝出奶茶这种能力），并且猜对了，那么她瞎猜对3杯，或者是比3杯更难发生的事件（例如猜对4杯）的概率是多少？现在我们做一个零假设，假设她能猜对4杯，这个事件的概率是多少？
现在我们就把这个事件再简化一下，我们做一个零假设，我们可以把女士品尝这个特殊案例看作是从一个含有4个绿球（正确答案）和4个红球（错误答案）的瓮中挑出4个小球的试验。
零假设就是：这个女士品尝出正确奶茶的结果纯粹是她自己瞎猜的，也就是说，每个小球有着相同的机会被挑中。现在计算一下每个事件的概率，那么挑中3个绿球的概率就是：
挑中4个绿球的概率是：
现在我们把这两个数值相应，就是约等于0.24，这个数字表示，这个女士能够品尝对3杯奶茶，甚至比3杯奶茶这个事件更加极端的概率（例如4杯奶茶），这就是p值。这个p值产生的过程就是Fisher精确检验(Fisher's exact test)，这类数据服从超几何分布(hypergeometric distribution)。
二联表
前面提到的数据我们可以使用二联表来创建，如下所示：
tab <- matrix(c(3,1,1,3),2,2)
rownames(tab)<-c("Poured Before","Poured After")
colnames(tab)<-c("Guessed before","Guessed after")
tab
数据如下所示：
> tab
Guessed before Guessed after
Poured Before 3 1
Poured After 1 3
使用fisher.test()这个函数就能计算出p值，如下所示：
fisher.test(tab,alternative="greater")
计算结果如下所示：
> fisher.test(tab,alternative="greater")
Fisher's Exact Test for Count Data
data: tab
p-value = 0.2429
alternative hypothesis: true odds ratio is greater than 1
95 percent confidence interval:
0.3135693 Inf
sample estimates:
odds ratio
6.408309
卡方检验
GWAS全称是Genome-wide association studies，即全基因组关联分析。这种分析的统计结果常用曼哈顿图(Manhattan)图来表示。曼哈顿图的y轴是p值的负的log值（底为10），这个p值则是通过计算数百万个SNP的关联检验得到的。x轴通常是染色体的位置（例如chromosome 1 to 22, X,Y等）。这引起p值的计算类似于前面提到的女士品茶的计算过程。但是，在女士品茶那个案例中，绿球与红球的数目是固定的，针对每种事件的发生次数也是固定的，换句话讲，在列联表中，行与列的中的数字之和是固定的。这种限定条件就对二联表中所有的可能性进行了约束，并且还可以使用超几何分布来计算p值。但是，实际情况并非如此，这时，可以使用另外一种方法，即卡方检验(Chi-squared test)。
假如我们有250个人，他们中有一些人患有某种疾病，剩下的则不患。在这些人中，纯合子(aa)有50人，这些纯合子中又有20%的人患病，即有10名患有该病，而非纯合子(AA/Aa)的人中则有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)),] #shuffle them up
head(dat)
结果如下所示：
> head(dat)
disease genotype
130 control AA/Aa
191 cases AA/Aa
187 cases AA/Aa
233 control aa
144 control AA/Aa
91 control AA/Aa
使用table()函数可以创建一个列联表，如下所示：
table(genotype)
结果如下所示：
> table(genotype)
genotype
AA/Aa aa
200 50
查看一下疾病的列联表，如下所示：
table(disease)
结果如下：
> table(disease)
disease
control cases
220 30
现在把基因型与疾病关联起来，如下所示：
tab <- table(genotype,disease)
tab
结果如下所示：
> tab
disease
genotype control cases
AA/Aa 180 20
aa 40 10
从上面结果可以看出来，我们使用了2个因素，即基因型，疾病，函数table()返回了2X2列联表。
对于上述的这种数据，我们通常会使用比数比(OR, odds ratio)来表示。现在我们计算一下aa的比数(odd)，即10/40，也就是说在50个纯合子(aa)中，有10名患病，40名不患病，这个数值是患病与不患病的比值。那么AA/Aa这种基因型患病的比数为20/180，这两个比值再进行比，就是比数比(odd ratio)，即：(10/40)/(20/180)，如下所示：
(tab[2,2]/tab[2,1]) / (tab[1,2]/tab[1,1])
结果如下所示：
> (tab[2,2]/tab[2,1]) / (tab[1,2]/tab[1,1])
[1] 2.25
我们并不使用OR来直接计算p值。相反，我们会先假设基因型与疾病没有关联，然后计算出表格中每个单元格的期望数值。在零假设下，含有200个人的AA/Aa和含有50个人的aa都被计算出相应的患病数，把他们都当成一个总体，那么患病率是：
p=mean(disease=="cases")
p
结果如下所示：
> p
[1] 0.12
那么这个表格的期望值是（期望值是使用每组的人数乘以发病率）：
expected <- rbind(c(1-p,p)*sum(genotype=="AA/Aa"),
c(1-p,p)*sum(genotype=="aa"))
dimnames(expected)<-dimnames(tab)
expected
结果如下所示：
> expected
disease
genotype control cases
AA/Aa 176 24
aa 44 6
卡方检验的p值计算结果如下：
chisq.test(tab)$p.value
结果为：
> chisq.test(tab)$p.value
[1] 0.08857435
大样本，小p值
前面提到，在报道一项发现时，仅报道p值是不合适。许多遗传学关联研究似乎过于强调p值，在这些研究中，研究者们使用了大量的样本，因此他们报道的p值非常人低。但是，如果我们仔细看一下他们的结果就会发现，比数比(odd ratio)则表现得很一般：很少有大于1的。在这种情况下，基因型为AA/Aa或aa的差异也许并不改变个体的疾病风险，虽然这有着实际意义(practically significant)，但是，一个人不可能根据这些基因型与疾病的关系而对自己的行为进行改变。
比数比(odds ration)与p值之间并非是一对一的关系。为了说明这一点，我们在保证所有比例相同的前提下来重新计算一下p值。现在我们把样本数目扩大10倍，这就会明显地降低p值，如下所示：
tab <- tab*10
chisq.test(tab)$p.value
计算结果如下所示：
> tab <- tab*10
> chisq.test(tab)$p.value
[1] 1.219624e-09
比数比的置信区间(Confidence Intervals For The Odd Ratio)
在数学上计算OR的置信区间并不简单明了。在计算其它统计量时，我们可以推导出那些统计量的近似，但是OR的计算与之不同，OR不仅是一个比值，而是一个比值的比值。因此，没有简单的方法（例如CLT）来计算OR。
一种方法是使用广义线性模型(generalized linear models )的理论，这种理论可以对比数比的log值(log odds ratio)进行估计。经过log转换的比数比就可以被证明是接近正态分布的，如下所示（具体的理论可以参考文献）：
fit <- glm(disease~genotype,family="binomial",data=dat)
coeftab<- summary(fit)$coef
coeftab
计算结果如下所示：
> 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
数据的表格中的第2行就是估计值与log(odds ratio)的值。数学理论告诉我们，这个估计值是接近于正态分布的。因此我们可以构建一个置信区间，然后再化为正常值，就能计算出OR的置信区间，如下所示：
ci <- coeftab[2,1] + c(-2,2)*coeftab[2,2]
exp(ci)
计算结果为：
> ci <- coeftab[2,1] + c(-2,2)*coeftab[2,2]
> exp(ci)
[1] 0.9618616 5.2632310
这个置信区间包括1，这与p值大于于0.05的结果是一致的。需要注意的是，这里的p值是基于另一种不同方法计算的，它以前面使用卡方检验得到的p值计算过程不一样。
练习
P91
