分割表〜遺伝情報活用のための遺伝統計学

  • ケース・コントロールスタディ
    • (患者、非患者) x (要因あり、なし)
---
title: "遺伝情報活用のための遺伝統計学"
author: "ryamada"
date: "2016年10月2日"
output: html_document
---

# 概論

遺伝的バリアントのリスクに関するデータを読み、それがもたらす遺伝的リスクを理解することを目標とする。

取得データのタイプにより推定指標は異なるが、その点推定と区間推定を行う。また事前確率と組み合わせたベイズ流の事後分布を推定することを原則とする。

複数のパターンを取り上げ、

* その事例検討
* Rを用いた統計解析の実践
* モデル設定した上でのデータのシミュレーション生成
* 生成データからの統計量・推定量分布の確認

パターンとして以下を行う。

* 1SNPのケース・コントロール解析
    + その性別影響の加味
    + その年齢影響の加味
* 1SNPの量的表現型解析
    + その性別影響の加味
    + その年齢影響の加味

# 検定と推定の基礎

## 基本確認〜これはやりませんが、解らない人は2x2表の独立性検定とp値について復習すること

2x2表を用いて、2要素の間の独立性の検定を行った。

帰無仮説が成り立っているときにサンプリングを繰り返して検定すると、p値は一様分布となる。
一様分布からの乱数は、度数分布にすると平坦になり、ソートしてプロットすると直線状になる。

```{r}
n.iter <- 10^4
N <- 1000

ps <- rep(0,n.iter)
for(i in 1:n.iter){
  x1 <- sample(0:1,N,replace=TRUE)
  x2 <- sample(0:1,N,replace=TRUE)
  ps[i] <- chisq.test(table(x1,x2),correct=FALSE)$p.value
}
hist(ps)
plot(sort(ps))
```


## 例1 2x2表データの取得と推定:ケース・コントロールスタディ 

ある集団があり、人口がNであるという。

今、ある調査をして、Ne人をランダムに調べたところ、二値型表現型Yが、ne0,ne1 (ne0+ne1=Ne)だったという。

また、この集団において、表現型Yのケース・コントロールスタディを実施し、ケース(Y=1)、コントロールをそれぞれ、mcase,mcontrol人集め、ある要因Xの有無について測定した。
ケース、コントロールのXの要因別の内訳は、mca0 + mca1=mcase,mco0+mco1=mcontrolとなったという。

この結果を用いて、ある個人について表現型Yが1となる確率について考えたい。

```{r}
ne0 <- 1000 - 70
ne1 <- 70
mca0 <- 500
mca1 <- 400
mco0 <- 700
mco1 <- 400
```

要因Xの情報がないのであれば、ne0,ne1の値を基に推定するしかない。

おおよそ
```{r}
ne1/(ne0+ne1)
```
がY=1の確率。

区間推定をすれば
```{r}
ps <- seq(from=0,to=1,length=10000)
pry1 <- dbeta(ps,ne1+1,ne0+1)
plot(ps,pry1,type="l")
plot(ps,pry1,type="l",xlim=c(0.005,0.008))
abline(v=ne1/(ne0+ne1),col=2)
library(binom)
binom.confint(ne1,ne0+ne1)
ci.bayes.central <- binom.confint(ne1,ne0+ne1,method="bayes",type="central",prior.shape1=1,prior.shape2=1)
abline(v=ci.bayes.central,col=3)
abline(v=ne1/(ne0+ne1),col=2)
```

関連検定をする。
```{r}
tab <- matrix(c(mca0,mca1,mco0,mco1),byrow=TRUE,2,2)
tab
chisq.test(tab)
```

関連はあるのだろう。


要因あり、要因なし、それぞれでのY=1の確率を推定する。

```{r}
pry1_x0 <- dbeta(ps,mco1+1,mco0+1)
pry1_x1 <- dbeta(ps,mca1+1,mca0+1)

matplot(ps,cbind(pry1_x0,pry1_x1),type="l")
```


今、Xの有無にかかわらず、Y=1である確率が、uであるとし、
Y=0のときのX=1である確率が、vであるとし、X=1のときのX=1である確率がwとすると、集団全体での、X=1であるのは、u*w + (1-u)*vであるから、

X=1であることがわかっているときにY=1であるのは
u*w/(u*w+(1-u)*v)となる。

今、u,v,wのそれぞれの値が、分布として推定されているから

```{r}
n.iter <- 10^5
u <- rbeta(n.iter,ne1+1,ne0+1)
v <- rbeta(n.iter,mca1+1,mca0+1)
w <- rbeta(n.iter,mco1+1,mco0+1)
ret <- u*v/(u*v+(1-u)*w)

ret. <- u*(1-v)/((u*(1-v)+(1-u)*(1-w)))
hist(ret)
```
```{r}
hist(ret,density=33)
hist(u,add=TRUE,col=2,density=25)
hist(ret.,add=TRUE,col=3,density=19)
summary(u)
summary(ret)
```
```{r}
matplot(log10(cbind(sort(u),sort(ret),sort(ret.))),type="l")
```

## 例1の解説

```{r}
ne0 <- 1000 - 70
ne1 <- 70
mca0 <- 500
mca1 <- 400
mco0 <- 700
mco1 <- 400
```

#### 割合の推定

標本割合を計算する。
```{r}
ne1/(ne0+ne1)
```

確率pでおきる事象をNe回観察し、成否がne1,ne0となる確率は二項分布。

Ne回観察し、その成否がne1,ne0となったときに、成否確率がpである尤度は、ベータ分布。

二項分布は離散的(折れ線グラフ)。ベータ分布は連続的(曲線グラフ)。

二項分布とベータ分布は、2次元プロットを見る方向の違い。

```{r}
N <- 10 # 総標本数
ne1s <- 0:N # 可能な成功回数
ps <- seq(from=0,to=1,length=100) # 成功率
pr <- matrix(0,length(ne1s),length(ps)) # 確率・尤度を納める行列
for(i in 1:length(ne1s)){
  pr[i,] <- dbinom(ne1s[i],N,ps)
}
par(mfcol=c(1,2))
plot(pr[3,],type="l")
plot(pr[,30],type="h")
par(mfcol=c(1,1))
persp(pr,phi=30,theta=5)
```

二項分布で算出したが、同様の分布はベータ分布としても計算できる
```{r}
out.binom <- dbinom(3,10,ps)
out.beta <- dbeta(ps,3+1,7+1) # 下駄をはかせることに注意
par(mfcol=c(2,2))
plot(ps,out.binom,type="l")
plot(ps,out.beta,type="l")
plot(out.binom,out.beta)
par(mfcol=c(1,1))
```

ベータ分布の最頻値(最尤推定値)、期待値(推定の期待値)、区間

```{r}
library(binom)
a <- 3
b <- 7
pr <- dbeta(ps,a+1,b+1)
plot(ps,pr,type="l")
# 最頻値〜標本平均
abline(v=a/(a+b),col=2)
# 期待値
abline(v=(a+1)/(a+1+b+1),col=3)
ci.bayes <- binom.bayes(a,a+b,type="highest",prior.shape1=1,prior.shape2=1)
abline(v=ci.bayes[6:8],col=4)
abline(h=dbeta(ci.bayes$lower,a+1,b+1),col=5)
dbeta(ci.bayes$lower,a+1,b+1)
dbeta(ci.bayes$upper,a+1,b+1)
```

## 何の分布を知りたいかは、場合によりけり

2つの分布を比べる

ケースとコントロールのそれぞれの群について、X=1の確率分布が得られた。

X=1である確率がケース・コントロールのどちらの方で高いのかは、2つの分布を比べる必要がある。
2つの分布を単純に比較するときは、同時分布というものを使う。
簡単に言えば、横軸にPr(X=1|Y=0)をとり、縦軸にPr(X=1|Y=1)を取って描ける、2次元平面の分布のこと。

Pr(X=1|Y=0),Pr(X=1|Y=1)はどちらもベータ分布であって、

```{r}
mca0 <- 5
mca1 <- 4
mco0 <- 70
mco1 <- 40
ps1 <- ps2 <- seq(from=0,to=1,length=100)
prx0 <- dbeta(ps1,mco1+1,mco0+1)
prx1 <- dbeta(ps2,mca1+1,mca0+1)
ps12 <- expand.grid(ps1,ps2)
tmp <- expand.grid(prx0,prx1)
joint.p <- tmp[,1]*tmp[,2]
persp(matrix(joint.p,length(ps1),length(ps2)))
image(matrix(joint.p,length(ps1),length(ps2)))
```

同時分布から何がわかるか?

Pr(X=1|Y=1) > Pr(X=1|Y=0)の確率
```{r}
image(matrix(joint.p,length(ps1),length(ps2)))
abline(0,1,col=1)
```

Pr(X=1|Y=1) > Pr(X=1|Y=0)*2 の確率
```{r}
image(matrix(joint.p,length(ps1),length(ps2)))
abline(0,2,col=1)
```

Pr(X=1|Y=1) > Pr(X=1|Y=0) + 0.2 の確率
```{r}
image(matrix(joint.p,length(ps1),length(ps2)))
abline(0.2,1,col=1)
```


一番知りたいところの、X=0,1ごとのY=1の確率は、有病率に依存するので、少し工夫が居る。

```{r}
n.iter <- 10^4
u <- rbeta(n.iter,ne1+1,ne0+1) # 有病率の分布に従う乱数
v <- rbeta(n.iter,mca1+1,mca0+1) # ケースのX=1確率
w <- rbeta(n.iter,mco1+1,mco0+1) # コントロールのX=1確率
# さまざまな有病率における、X=1の場合のY=1率
ret <- u*v/(u*v+(1-u)*w)
# 対応する有病率における、X=0の場合のY=1率
ret. <- u*(1-v)/((u*(1-v)+(1-u)*(1-w)))
# Pr(Y=1|X=0),Pr(Y=1|X=1)の散布図
plot(ret,ret.,pch=20,cex=0.1,xlim=c(0,1),ylim=c(0,1))
```

何がわかるか。

Pr(Y=1|X=1) > Pr(Y=1|X=0)の確率
```{r}
plot(ret,ret.,pch=20,cex=0.1,xlim=c(0,1),ylim=c(0,1))
abline(0,1,col=2)
```
Pr(Y=1|X=1) > Pr(Y=1|X=0)*2の確率
```{r}
plot(ret,ret.,pch=20,cex=0.1,xlim=c(0,1),ylim=c(0,1))
abline(0,2,col=2)
```

Pr(Y=1|X=1) > Pr(Y=1|X=0)*+0.2の確率
```{r}
plot(ret,ret.,pch=20,cex=0.1,xlim=c(0,1),ylim=c(0,1))
abline(0.2,1,col=2)
```

Pr(Y=1|X=1)/Pr(Y=1|X=0)の分布
```{r}
hist(ret/ret.)
```
Pr(Y=1|X=1)-Pr(Y=1|X=0)の分布
```{r}
hist(ret-ret.)
```

Xの情報を得る前と後とのY=1確率の違い
Pr(Y=1|X=1) / Pr(Y=1|X=0 or 1)
Pr(Y=1|X=0) / Pr(Y=1|X=0 or 1)
```{r}
hist(ret/u)
hist(ret./u)
```

## 発展

### 共役事前分布

二項分布に従う事象について、ベータ分布を事後分布として取った。
この二項分布とベータ分布はベイズ推定をするときにとてもうまくできた関係にある。

他にも観察事象をA分布、それのモデルとなる確率分布がB分布となっているような場合がいくつもある。

このようなAとBの関係になっているときに、BをAの共役事前分布と呼ぶ。

### 無情報事前分布

dbeta(p,a+1,b+1)+1について。
式を確認。

この+1は、尤度の計算のためのもの。

事後分布 = 尤度 x 事前分布

と言う関係があり、事前分布を定める必要がある。

成功率pが一様分布であるという事前仮定のときには、観察が(0,0)のときに、dbeta(p,0+1,0+1)としてやると、これが一様分布となるから、
pが一様であるという事前仮定を体現しているのが+1。

ただし、pが一様〜どんな値pも同じくらいの程度で「ありそう〜likely」とは言うものの、pというパラメタの取り方をこうしなくてはいけないのか!という一見変な突っ込みをし始めると、「何も情報がないとき」に、どういう事前分布を仮定するべきなのか、ということが解らなくなってくる。

そんな「何も情報がないときの事前分布」としてある特徴を持っていて、それなりに意味がありそうなものにJeffrenys priorというものがある。

ベータ分布の場合には、+0.5することに対応する。

```{r}
ps <- seq(from=0,to=1,length=100)
jeffreys <- dbeta(ps,0.5,0.5)
plot(ps,jeffreys,type="l")
```

全く成功しないか、ほぼ確実に成功するという場合を大きく見積もっている仮説である。

この事前分布を仮定すると、同じ件数の観察をすると、真の成功率がどれであっても、観察後の推定の曖昧さが同程度に縮む、という性質もある。

情報は隠された真によらない、という立場での無情報事前分布とも言える。

# 1SNPのケースコントロールスタディの場合


区間推定をすれば
```{r}
ps <- seq(from=0,to=1,length=10000)
pry1 <- dbeta(ps,ne1+1,ne0+1)
plot(ps,pry1,type="l")
plot(ps,pry1,type="l",xlim=c(0.005,0.008))
abline(v=ne1/(ne0+ne1),col=2)
library(binom)
binom.confint(ne1,ne0+ne1)
ci.bayes.central <- binom.confint(ne1,ne0+ne1,method="bayes",type="central",prior.shape1=1,prior.shape2=1)
abline(v=ci.bayes.central,col=3)
abline(v=ne1/(ne0+ne1),col=2)
```

関連検定をする。
```{r}
tab <- matrix(c(mca0,mca1,mco0,mco1),byrow=TRUE,2,2)
tab
chisq.test(tab)
```

関連はあるのだろう。


要因あり、要因なし、それぞれでのY=1の確率を推定する。

```{r}
pry1_x0 <- dbeta(ps,mco1+1,mco0+1)
pry1_x1 <- dbeta(ps,mca1+1,mca0+1)

matplot(ps,cbind(pry1_x0,pry1_x1),type="l")
```


今、Xの有無にかかわらず、Y=1である確率が、uであるとし、
Y=0のときのX=1である確率が、vであるとし、X=1のときのX=1である確率がwとすると、集団全体での、X=1であるのは、u*w + (1-u)*vであるから、

X=1であることがわかっているときにY=1であるのは
u*w/(u*w+(1-u)*v)となる。

今、u,v,wのそれぞれの値が、分布として推定されているから

```{r}
n.iter <- 10^5
u <- rbeta(n.iter,ne1+1,ne0+1)
v <- rbeta(n.iter,mca1+1,mca0+1)
w <- rbeta(n.iter,mco1+1,mco0+1)
ret <- u*v/(u*v+(1-u)*w)

ret. <- u*(1-v)/((u*(1-v)+(1-u)*(1-w)))
hist(ret)
```
```{r}
hist(ret,density=33)
hist(u,add=TRUE,col=2,density=25)
hist(ret.,add=TRUE,col=3,density=19)
summary(u)
summary(ret)
```
```{r}
matplot(log10(cbind(sort(u),sort(ret),sort(ret.))),type="l")
```
N人のうち、実際には、N0+N1=Nであるが、この真の値は分からない。


要因Xを有するとき、表現型Yが1となる確率がp1であり、要因Xを有しないとき、表現型Yが1となる確率がp0とする。

今、集団全体で、要因Xを有する割合qとする。

集団全体では
```{r}
p1 <- 0.2
p0 <- 0.1
q <- 0.2

tab <- matrix(c((1-q)*(1-p0),(1-q)*p0,q*(1-p1),q*p1),byrow=TRUE,2,2)
tab
apply(tab,1,sum)
```

集団全体では、(X,Y) = (0,0)0.891(0,1)0.099(1,0)0.008(1,1)0.002の割合でいることになる。

Y=(0,1)の割合(非病気・病気の割合)は
```{r}
y01 <- apply(tab,2,sum)
y01
```
となっている。

ケース・コントロールスタディをして、ケース500人、コントロール1000人を集めたとする。

```{r}
tab.casecontrol <- cbind(tab[,1]/y01[1]*1000,tab[,2]/y01[2]*500)

tab.casecontrol
apply(tab.casecontrol,2,sum)
```


```{r}
chisq.test(tab)
```

要因xの保持と表現型Yの保持に関係があるらしいことが解った。

## 推定

p0,p1が不明。
p0,p1の値がを推定したい。

さらには、Xが0だったとき(Xが1だったとき)、Yが1である確率も推定したい。

ベイズ推定。

p0=tであるときに、1000人中100人でY=1となる確率は

```{r}
t <- 0.5
dbinom(100,1000,t) # すごく低い確率
t <- 0.1
dbinom(100,1000,t) # 高い確率
```
1000人中100人でY=0という観察が得られたときに、p0=tである尤度、とも言う。

実際には、tの値は0から1までどんな値かわからないので
```{r}
ts <- seq(from=0,to=1,length=10000)
likes <- dbinom(100,1000,ts)
plot(ts,likes,type="l")
abline(v=0.1,col=2)
```
このプロットの山のもっとも高くするtの値がp0の最尤推定値。今回の場合は0.1。

拡大してみる。

尤度の高さを基準に区間に組み込むかどうかを判断して、95%が組み込まれるようにしてみる。

```{r}
plot(ts,likes,type="l")
plot(ts,likes,type="l",xlim=c(0.08,0.13),ylim=c(0,0.01))
abline(v=0.1,col=2)
library(binom)
binom.confint(100,1000)
ci.bayes.central <- binom.confint(100,1000,methods="bayes",type="central")
ci.bayes.central
#ci.bayes.highest <- binom.confint(100,1000,methods="bayes",type="highest")
#ci.bayes.highest
abline(v=ci.bayes.central,col=3)
#abline(v=ci.bayes.highest,col=4)
```


# 分割表
```{r}
tab <- matrix(c(900,4200,4900,1100,4500,4700),byrow=TRUE,2,3)
```

## 複数の検定
```{r}
chisq.test(tab,correct=FALSE)
prop.trend.test(tab[1,],tab[1,]+tab[2,])
prop.trend.test(tab[1,],tab[1,]+tab[2,],score=c(1,0.5,0))
prop.trend.test(tab[1,],tab[1,]+tab[2,],score=c(1,1,0))
prop.trend.test(tab[1,],tab[1,]+tab[2,],score=c(1,0,0))
```

ジェノタイプとフェノタイプとが独立でないことが解った。

## 推定する

推定しよう。

さて、何を



# 回帰

# マルチプルテスティング