漸近近似ベイズ因子 - ryamadaの遺伝学・遺伝統計学メモ

こちらにGWASでのSNP関連検定に漸近近似ベイズ因子を用いる論文について書いた。
以下、20080821と重複になる部分もあるが、漸近近似に中心を移して書く。
事前確率分布のおき方によっては、解析的にベイズ因子を算出するのは難しい。
適当な事前確率分布を仮定して、計算を楽にしよう。
必要な事前情報
- リスクローカスがあったとして、そのローカスの「相対危険度」の強さに関する分布
- データのサンプリングバイアス
リスクの事前分布
- リスクローカスのリスクの分布は、相対危険度 $r=e^{\theta}$ としたときの $\theta$ が平均 $0$ の正規分布になると仮定する(こうすると、データのサンプリングバイアスの分布も、正規分布近似をすることができるので、楽だし、そのように仮定することに、『生物学的実感上、無理がないから』)。
- $\theta \sim N(0,W)$ であるから、 $W$ を与える必要がある。リスクアレルのf(たとえば95%)が $\frac{1}{g} \le r \le g$ ( $-log(g) \le \theta log(g)$ )にあると仮定すると、Rで言えば、以下のような関係になる、そんなWを仮定していることになる。

pnorm(log(g),mean=0,sd=sqrt(W),lower.tail=FALSE)*2=1-f

$W=\frac{log(g)}{qnorm(\frac{1+f}{2},mean=0,sd=1)}$

getW<-function(g,f){
ret<-log(g)/qnorm((1+f)/2,mean=0,sd=1)
return(ret)
}

データのサンプリングバイアス
- データのサンプリングによる分布は、２ｘ２表の場合には、一般に、以下の式で与えることになる。 $V=\frac{1}{n1}(\frac{1}{n11}+\frac{1}{n12})+\frac{1}{n2}(\frac{1}{n21}+\frac{1}{n22})$
- SNPの２ｘ３表に対する、ロジスティック回帰(これだと、リスクをとしているので、このモデルによる扱いが適当)のときには、統計パッケージで得られる。
  - 例として、以下の、Coefficients: の G のEstimated の 0.5073 が $\hat{\theta}$ 、次のカラムのStd. Errorが $\sqrt{V}$ 、次のz valueが $Z$ 。

data<-read.table(file="infile.txt")
P<-data[,1];G<-data[,2];
result<-glm(formula = P ~ G, family = binomial)
summary(result)

Call:
glm(formula = P ~ G, family = binomial)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.36225  -1.14255   0.03196   1.21268   1.43593  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  -0.5898     0.2627  -2.246  0.02473 * 
G            0.5073     0.1880   2.699  0.00696 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 277.26  on 199  degrees of freedom
Residual deviance: 269.73  on 198  degrees of freedom
AIC: 273.73

Number of Fisher Scoring iterations: 4