title: DALS007-統(tǒng)計(jì)推斷(Inference)06-關(guān)聯(lián)檢驗(yàn)(Association Tests)
date: 2019-08-07 12:0:00
type: "tags"
tags:
- 統(tǒng)計(jì)推斷
categories: - 生物統(tǒng)計(jì)
前言
前面所提到的案例都是普通的統(tǒng)計(jì)學(xué)案例，這些普通案例指的是键袱，它們的數(shù)據(jù)是二元的(binary)俊抵，分類的(categorical)或有序的(ordinal)。現(xiàn)在我們考慮一種特殊的案例淋肾，例如基因型數(shù)據(jù)柳洋，現(xiàn)在我們有兩組基因型數(shù)據(jù)毕箍，一組基因是AA/Aa，另外一組是aa（某個(gè)特定疾病的對(duì)照）。現(xiàn)在我們的問題是胰耗，這個(gè)基因型(genotype)是否與疾病相關(guān)。這個(gè)案例與我們前面提到的案例就有點(diǎn)像了查描，我們有兩個(gè)總體，即AA/Aa與aa匠襟，我們使用1或0來表示它們對(duì)應(yīng)的疾病狀態(tài)啃勉，其中0是不患病，1是患病的。此時(shí)我們不能使用t檢驗(yàn)先改，因?yàn)檫@種數(shù)據(jù)明顯不符合正態(tài)分布，無法使用t檢驗(yàn)蒸走。如果樣本足夠大的話仇奶，可以使用CLT，否則比驻，可以使用關(guān)聯(lián)檢驗(yàn)(association test)该溯。
女士品茶
關(guān)于假設(shè)檢驗(yàn)最有名的一個(gè)案例就是Fisher提出的。這個(gè)案例大概是這個(gè)樣子的别惦，一名女士聲稱狈茉，如果牛奶與紅茶加入的先后順序不同，得到的奶茶口味也不同掸掸，并且她自己能品嘗出來氯庆。Fisher就給了4杯奶茶，這4杯奶茶加入的花與牛奶順序是不同的扰付，是隨機(jī)的堤撵。后來這個(gè)女士品嘗對(duì)了3個(gè)，此時(shí)羽莺，我們能認(rèn)為這個(gè)女士確實(shí)有這種能力么实昨，她品嘗出的這3杯奶茶是否是出于偶然因素？假設(shè)檢驗(yàn)就能回答這個(gè)問題盐固，并且能夠?qū)@種偶然因素進(jìn)行定量荒给。
現(xiàn)在我們提出我們的基礎(chǔ)問題：假如這個(gè)女士品嘗出的這3杯奶茶是她自己瞎猜（瞎猜說明她沒有品嘗出奶茶這種能力）丈挟，并且猜對(duì)了，那么她瞎猜對(duì)3杯锐墙，或者是比3杯更難發(fā)生的事件（例如猜對(duì)4杯）的概率是多少礁哄？現(xiàn)在我們做一個(gè)零假設(shè)，假設(shè)她能猜對(duì)4杯溪北，這個(gè)事件的概率是多少桐绒？
現(xiàn)在我們就把這個(gè)事件再簡(jiǎn)化一下，我們做一個(gè)零假設(shè)之拨，我們可以把女士品嘗這個(gè)特殊案例看作是從一個(gè)含有4個(gè)綠球（正確答案）和4個(gè)紅球（錯(cuò)誤答案）的甕中挑出4個(gè)小球的試驗(yàn)茉继。
零假設(shè)就是：這個(gè)女士品嘗出正確奶茶的結(jié)果純粹是她自己瞎猜的，也就是說蚀乔，每個(gè)小球有著相同的機(jī)會(huì)被挑中∷附撸現(xiàn)在計(jì)算一下每個(gè)事件的概率，那么挑中3個(gè)綠球的概率就是：
挑中4個(gè)綠球的概率是：
現(xiàn)在我們把這兩個(gè)數(shù)值相應(yīng)吉挣，就是約等于0.24派撕，這個(gè)數(shù)字表示，這個(gè)女士能夠品嘗對(duì)3杯奶茶睬魂，甚至比3杯奶茶這個(gè)事件更加極端的概率（例如4杯奶茶）终吼，這就是p值。這個(gè)p值產(chǎn)生的過程就是Fisher精確檢驗(yàn)(Fisher's exact test)氯哮，這類數(shù)據(jù)服從超幾何分布(hypergeometric distribution)际跪。
二聯(lián)表
前面提到的數(shù)據(jù)我們可以使用二聯(lián)表來創(chuàng)建，如下所示：
tab <- matrix(c(3,1,1,3),2,2)
rownames(tab)<-c("Poured Before","Poured After")
colnames(tab)<-c("Guessed before","Guessed after")
tab
數(shù)據(jù)如下所示：
> tab
Guessed before Guessed after
Poured Before 3 1
Poured After 1 3
使用fisher.test()這個(gè)函數(shù)就能計(jì)算出p值喉钢，如下所示：
fisher.test(tab,alternative="greater")
計(jì)算結(jié)果如下所示：
> 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
卡方檢驗(yàn)
GWAS全稱是Genome-wide association studies姆打，即全基因組關(guān)聯(lián)分析。這種分析的統(tǒng)計(jì)結(jié)果常用曼哈頓圖(Manhattan)圖來表示肠虽。曼哈頓圖的y軸是p值的負(fù)的log值（底為10）幔戏，這個(gè)p值則是通過計(jì)算數(shù)百萬個(gè)SNP的關(guān)聯(lián)檢驗(yàn)得到的。x軸通常是染色體的位置（例如chromosome 1 to 22, X,Y等）舔痕。這引起p值的計(jì)算類似于前面提到的女士品茶的計(jì)算過程评抚。但是，在女士品茶那個(gè)案例中伯复，綠球與紅球的數(shù)目是固定的慨代，針對(duì)每種事件的發(fā)生次數(shù)也是固定的，換句話講啸如，在列聯(lián)表中侍匙，行與列的中的數(shù)字之和是固定的。這種限定條件就對(duì)二聯(lián)表中所有的可能性進(jìn)行了約束，并且還可以使用超幾何分布來計(jì)算p值想暗。但是妇汗，實(shí)際情況并非如此，這時(shí)说莫，可以使用另外一種方法杨箭，即卡方檢驗(yàn)(Chi-squared test)。
假如我們有250個(gè)人储狭，他們中有一些人患有某種疾病互婿，剩下的則不患。在這些人中辽狈，純合子(aa)有50人慈参，這些純合子中又有20%的人患病，即有10名患有該病刮萌，而非純合子(AA/Aa)的人中則有10%患有該病驮配。如果我們?cè)偬?50個(gè)人，那么又是什么情況呢着茸？
現(xiàn)在我們先根據(jù)既有的數(shù)據(jù)來過計(jì)算一下壮锻，如下所示：
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)
結(jié)果如下所示：
> 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()函數(shù)可以創(chuàng)建一個(gè)列聯(lián)表，如下所示：
table(genotype)
結(jié)果如下所示：
> table(genotype)
genotype
AA/Aa aa
200 50
查看一下疾病的列聯(lián)表涮阔，如下所示：
table(disease)
結(jié)果如下：
> table(disease)
disease
control cases
220 30
現(xiàn)在把基因型與疾病關(guān)聯(lián)起來躯保，如下所示：
tab <- table(genotype,disease)
tab
結(jié)果如下所示：
> tab
disease
genotype control cases
AA/Aa 180 20
aa 40 10
從上面結(jié)果可以看出來，我們使用了2個(gè)因素澎语，即基因型，疾病验懊，函數(shù)table()返回了2X2列聯(lián)表擅羞。
對(duì)于上述的這種數(shù)據(jù)，我們通常會(huì)使用比數(shù)比(OR, odds ratio)來表示∫逋迹現(xiàn)在我們計(jì)算一下aa的比數(shù)(odd)减俏，即10/40，也就是說在50個(gè)純合子(aa)中碱工，有10名患病娃承，40名不患病，這個(gè)數(shù)值是患病與不患病的比值怕篷。那么AA/Aa這種基因型患病的比數(shù)為20/180历筝，這兩個(gè)比值再進(jìn)行比，就是比數(shù)比(odd ratio)廊谓，即：(10/40)/(20/180)梳猪，如下所示：
(tab[2,2]/tab[2,1]) / (tab[1,2]/tab[1,1])
結(jié)果如下所示：
> (tab[2,2]/tab[2,1]) / (tab[1,2]/tab[1,1])
[1] 2.25
我們并不使用OR來直接計(jì)算p值。相反蒸痹，我們會(huì)先假設(shè)基因型與疾病沒有關(guān)聯(lián)春弥，然后計(jì)算出表格中每個(gè)單元格的期望數(shù)值呛哟。在零假設(shè)下，含有200個(gè)人的AA/Aa和含有50個(gè)人的aa都被計(jì)算出相應(yīng)的患病數(shù)匿沛，把他們都當(dāng)成一個(gè)總體扫责，那么患病率是：
p=mean(disease=="cases")
p
結(jié)果如下所示：
> p
[1] 0.12
那么這個(gè)表格的期望值是（期望值是使用每組的人數(shù)乘以發(fā)病率）：
expected <- rbind(c(1-p,p)*sum(genotype=="AA/Aa"),
c(1-p,p)*sum(genotype=="aa"))
dimnames(expected)<-dimnames(tab)
expected
結(jié)果如下所示：
> expected
disease
genotype control cases
AA/Aa 176 24
aa 44 6
卡方檢驗(yàn)的p值計(jì)算結(jié)果如下：
chisq.test(tab)$p.value
結(jié)果為：
> chisq.test(tab)$p.value
[1] 0.08857435
大樣本，小p值
前面提到逃呼，在報(bào)道一項(xiàng)發(fā)現(xiàn)時(shí)鳖孤，僅報(bào)道p值是不合適。許多遺傳學(xué)關(guān)聯(lián)研究似乎過于強(qiáng)調(diào)p值蜘渣，在這些研究中淌铐，研究者們使用了大量的樣本，因此他們報(bào)道的p值非常人低蔫缸。但是腿准，如果我們仔細(xì)看一下他們的結(jié)果就會(huì)發(fā)現(xiàn)，比數(shù)比(odd ratio)則表現(xiàn)得很一般：很少有大于1的拾碌。在這種情況下吐葱，基因型為AA/Aa或aa的差異也許并不改變個(gè)體的疾病風(fēng)險(xiǎn)，雖然這有著實(shí)際意義(practically significant)校翔，但是弟跑，一個(gè)人不可能根據(jù)這些基因型與疾病的關(guān)系而對(duì)自己的行為進(jìn)行改變。
比數(shù)比(odds ration)與p值之間并非是一對(duì)一的關(guān)系防症。為了說明這一點(diǎn)孟辑，我們?cè)诒ＷC所有比例相同的前提下來重新計(jì)算一下p值。現(xiàn)在我們把樣本數(shù)目擴(kuò)大10倍蔫敲，這就會(huì)明顯地降低p值饲嗽，如下所示：
tab <- tab*10
chisq.test(tab)$p.value
計(jì)算結(jié)果如下所示：
> tab <- tab*10
> chisq.test(tab)$p.value
[1] 1.219624e-09
比數(shù)比的置信區(qū)間(Confidence Intervals For The Odd Ratio)
在數(shù)學(xué)上計(jì)算OR的置信區(qū)間并不簡(jiǎn)單明了。在計(jì)算其它統(tǒng)計(jì)量時(shí)奈嘿，我們可以推導(dǎo)出那些統(tǒng)計(jì)量的近似貌虾，但是OR的計(jì)算與之不同，OR不僅是一個(gè)比值裙犹，而是一個(gè)比值的比值尽狠。因此，沒有簡(jiǎn)單的方法（例如CLT）來計(jì)算OR叶圃。
一種方法是使用廣義線性模型(generalized linear models )的理論袄膏，這種理論可以對(duì)比數(shù)比的log值(log odds ratio)進(jìn)行估計(jì)。經(jīng)過log轉(zhuǎn)換的比數(shù)比就可以被證明是接近正態(tài)分布的盗似，如下所示（具體的理論可以參考文獻(xiàn)）：
fit <- glm(disease~genotype,family="binomial",data=dat)
coeftab<- summary(fit)$coef
coeftab
計(jì)算結(jié)果如下所示：
> 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
數(shù)據(jù)的表格中的第2行就是估計(jì)值與log(odds ratio)的值哩陕。數(shù)學(xué)理論告訴我們，這個(gè)估計(jì)值是接近于正態(tài)分布的。因此我們可以構(gòu)建一個(gè)置信區(qū)間悍及，然后再化為正常值闽瓢，就能計(jì)算出OR的置信區(qū)間，如下所示：
ci <- coeftab[2,1] + c(-2,2)*coeftab[2,2]
exp(ci)
計(jì)算結(jié)果為：
> ci <- coeftab[2,1] + c(-2,2)*coeftab[2,2]
> exp(ci)
[1] 0.9618616 5.2632310
這個(gè)置信區(qū)間包括1心赶，這與p值大于于0.05的結(jié)果是一致的扣讼。需要注意的是，這里的p值是基于另一種不同方法計(jì)算的缨叫，它以前面使用卡方檢驗(yàn)得到的p值計(jì)算過程不一樣椭符。
練習(xí)
P91
