ぱらぱらめくる『Causal Inference for Statistics, Social, and Biomedical Sciences』

Causal Inference for Statistics, Social, and Biomedical Sciences: An Introduction

Causal Inference for Statistics, Social, and Biomedical Sciences: An Introduction

  • 目次
  • Part I Introduction
  • Part II Classical Randomized Experiments
  • Part III Regular Assignment Mechanisms: Design
  • Part IV Regular Assignment Mechanisms: Analysis
  • Part V Regular Assignment Mechanisms: Supplementary Analyses
  • Part VI Regular Assignment Mechanisms with Noncompliance: Analysis
  • Part VII Conclusion
  • まず、「アサインメント(割り付け)をどうするか」が問題であることがわかる。割り付けのバリエーションの中にRandomizationがあり、Classicalと表現されている。割り付けの色々と"Causal Inference"との間柄を説明するための、道具立て(割り付けの異同、Causal Inferenceの異同の表現手段とその間柄)をわかることがこの本のテーマらしい

Parts III,IV,V,VI 各論ぱらぱらめくる『Causal Inference for Statistics, Social, and Biomedical Sciences』

  • 割り付けがどうなされているか・観察データの得られ方がどうなっているかによって、スタディを分類し、それぞれに対する手法があるので、それについての説明になっている
  • 分類と帰属手法名とを目次を見て列挙するとよいだろう
  • Part II Randomized Assignment
    • フィッシャーの正確確率、ノイマンの繰り返しサンプリング、回帰、モデルあてはめ/推定、Randomizedだけど層別する場合、ペアワイズ
  • Parts III & IV & VI Regular Assignment mechanisms: Design & Analysis
    • Propensity Score(Prop.Scについて書かれた『読み物的理解用の記事』)、Covariate Distributionの扱い、マッチングする、Covariatesを扱いやすくするためのTrimming、小分けにする(subclassification)、Causal effect値を推定すること、一般的にCausal Estimandsの点推定/区間推定のためのこと、その他、confoundingについてわかることがあるなら取り出そう
  • Part VI Regular Assignment mechanisms with noncompliance: analysis
    • ここは「個別方法〜non-compliance対策〜」的な章

Part I Introductionぱらぱらめくる『Causal Inference for Statistics, Social, and Biomedical Sciences』

  • 1 Causality: The Basic Framework
  • 2 A Brief History of the Potential Outcomes Approach to Causal Inference
  • 3 A Classification of Assignment Mechanisms
  • この本で扱う"Causality"について定義
  • "Potential Outcomes" Approach と言うのがこの本の肝
  • 割り付けとは、割り付けがしていることとは、もこの本の構成上の基本事項
  • この1,2,3章がわかれば、残りのParts II, III, IV, V, VIは各論でPart VIIが、Introduction3章と対を成すまとめなので、1,2,3章がつかめてPart VIIがわかれば、この本の大意はつかめたことになりそう
  • 1 Causality
    • 1.1 Introduction
    • 1.2 Potential Outcomes
      • "Unit" に"Action"が作用して、別の"Unit"に変わる。変わりうる"Unit"状態は複数あり得るが実現するのは1つ。この変わりうる"Unit"状態がPotential Outcome
    • 1.3 Definition of Causal Effects
      • Causal effectの定義と、(定義に沿った)Causal effectについての推定の違いについて
      • 複数のPotential outcomesがあるときにあるoutcomeが出やすく他のoutcomeが出にくいことがActionのCausal effect
    • 1.4 Causal Effects in Common Usage
      • 1Unitの観察ではActionとoutcome(potential outcomesのうちの1つ)のペアしかないので、「そー、それで」となる。なので複数のUnitについてActionの具合とOutcomeの具合とのペアを取ることとする
    • 1.5 Learning about Causal Effects: Multiple Units
      • 複数のAction-outcomeのペアは、別物なので、それを「合わせる」には取決め・仮定が必要
    • 1.6 The Stable Unit Treatment Value Assumption(SUTVA)
      • 単純に確率・組み合わせ計算を用いることができる条件仮定
    • 1.7 The Assignment Mechanism: An Introduction
      • 割り付けによってAction選択のぶれをコントロールする
    • 1.8 Attributes, Pre-Treatment Variables, or Covariates
      • Unitsの違いを取り扱うために投入する情報たち
    • 1.9 Potential Outcomes and Lord's Paradox
    • 1.10 Causal Estimanads
      • Causalityに関して推定する統計量について
    • 1.11 Structure of the Book
      • Part I がイントロ
      • Part II がRandomized Assignmentで、解釈が簡単な例。割り付けが確率的にきちんと扱える
      • Parts III, IV 割り付けをコントロールできず、「こうわりつけちゃったけど(割り付け確率がわからないので、推定するとか、推定的に考えるとかしないといけない)、Action-Outcomeのペアはあって、Potential outcomesについて比較可能にはなっている」という状況の例
        • 観察研究とかがこれ
      • Parts V, VIでは、割り付けが共変量に依存しており、outcomeと結びついている場合(割り付けがどういう確率過程かも推定するとか推定的に考えるとかしないといけないけれど、outcome情報の得られ方にもそれを考える必要がある)
        • 観察研究の多くはこれ
    • 1.12 Samples, Populations, and Super-Populations
    • 1.13 Conclusion
  • 第2,3章の目次は以下のようになるが、第1章に関する上記の知識で十分(な人も多いはず)。
  • 2 A Brief History of the Potential Outcomes Approach to Causal Inference
    • 2.1 Introduction
    • 2.2 Potential Outcomes and the Assignment Mechanism before Neyman
    • 2.3 Neyman's (1923) Potential Outcome Notation in Randomized Experiments
    • 2.4 Earlier Hints for Physical Randomizing
    • 2.5 Fisher's (1925) Proposal to Randomize Treatments to Units
    • 2.6 The Observed Outcome Notation in Observational Studies for Causal Effects
    • 2.7 Early Uses of Potential Outcomes in Observational Studies in Social Sciences
    • 2.8 Potential Outcomes and the Assignment Mechanism in Observational Studies: Rubin (1974)
  • 3 A Classification of Assignment Mechanisms
    • 3.1 Introduction
    • 3.2 Notation
    • 3.3 Assignment Probabilities
    • 3.4 Restrictions on the Assignment Mechanism
    • 3.5 Assignment Mechanisms and Super-Populations
    • 3.6 Randomized Experiments
    • 3.7 Observational Studies: Regular Assignment Mechanisms
    • 3.8 Observational Studies: Irregular Assignment Mechanisms
    • 3.9 Conclusion

マルチプルテスティングとFDR〜オミックス統計学入門2014

  • 2014年度講義資料
  • FDRとマルチプルテスティング
  • 90分で1-2回分相当
  • Rmdファイルです。html化、epub化できます(やり方はこちら)
  • html化、epub化が面倒くさければ、kindleで1米ドルでも(こちら)

# マルチプルテスティングとFalse Discovery Rate

# 検定とp値と、どれくらい信じるか(帰無仮説の棄却)

## 一様分布とp値

一様分布というのは、どんな値も同じ確率で出現する分布である。
横軸が0から1、縦軸も0から1の正方形になっている。
x軸の値によらず、出現確率が同じ高さである。
検定で用いられるp値はこのような性質を持っている。
小さいp値も大きいp値も中くらいのp値も同じ確率で出現する。


```{r}
x <- seq(from=0,to=1,length=10000)
d <- dunif(x)
plot(x,d,type="l",ylim=c(0,1.1),col=2)
abline(v=c(0,1))
abline(h=0)
```

たとえば0.05より小さいp値は次の図の赤い部分であって、全体の0.05を占める。


```{r}
x <- seq(from=0,to=1,length=10000)
col <- rep(1,length(x))
col[which(x<0.05)] <- 2
plot(x,d,type="h",ylim=c(0,1.1),col=col)
abline(v=c(0,1))
abline(h=0)
```

従って、p値が0.05であるということは、「p値が0.05以下になる確率が0.05である」ということを意味している。
これを別の方法で示す。


```{r}
p <- punif(x)
plot(x,p,type="p",cex=0.01,ylim=c(0,1.1),col=col)
abline(v=c(0,0.05,1))
abline(h=c(0,0.05,1))
```

p値を小さい方から累積してプロットすると、グラフは$y=x$の直線となり、原点と(1,1)の2点を結んだものとなる。
p値が0.05であるということは、この直線の左側0.05の部分を登ったところであることを示している。

## 帰無仮説を棄却する

ある試料を一斉に集め、適当に分配して、A・B2社に外注したところ、A社とB社とのそれぞれからの納品濃度の分布が次のようになっていたとする。
A社とB社とで納品濃度に違いがないと考えるかどうかを試してみる。

2社の納品濃度の分布が同じであると仮定してデータを作成し、t検定をしてみる。
同じことを仮定しているので、帰無仮説は棄却されないはずである。
p値は0.05より小さくならないだろうと予想されると言うこともできる。

```{r}
ma <- 120
sa <- 10
mb <- 120
sb <- 10
na <- 20
nb <- 15
A <- rnorm(na,ma,sa)
B <- rnorm(nb,mb,sb)
h <- hist(c(A,B),plot=FALSE)
hist(A,breaks=h$breaks,ylim=h$ylim,density=10)
hist(B,breaks=h$breaks,density=7,col=2,add=TRUE)
t.test(A,B)$p.value
```

小さいp値が出ないだろうという予想は、高い確率で裏切られませんが、
低い確率では裏切られる。
それを試してみる。
```{r}
n.iter <- 10000
ps <- rep(0,n.iter)
for(i in 1:n.iter){
  A <- rnorm(na,ma,sa)
  B <- rnorm(nb,mb,sb)
  ps[i] <- t.test(A,B)$p.value
}
hist(ps)
plot(sort(ps))
```
たしかに、p値の分布は一様分布らしくなっており、0.05より小さいp値は確率0.05で出現する。

## たくさんの検定〜GWAS〜

ゲノムワイドに50万個のSNPと誕生月が偶数か奇数かとの関係を調べるとする。
無関係な50万の検定が繰り返される。


```{r}
n.snp <- 500000
p <- runif(n.snp)
plot(p,cex=0.1)
```
点の数が多すぎて何のことだかわからないので、対数表示にしてやる。
こうすると、いわゆるGWASのマンハッタンプロットになる

```{r}
plot(-log10(p),cex=0.1,ylim=c(0,10))
```

一番小さいp値はいくつだろうか。
プロット上で強調しておく。
```{r}
min(p)
plot(-log10(p),cex=0.1,ylim=c(0,10))
minp <- which(p==min(p))
points(minp,-log10(p[minp]),col=2,pch=20)
```

今、p値はすべて帰無仮説が成り立つ状態で発生させたのであるから、
この小さいp値を見ても、

> 帰無仮説が成り立っているときの、ごく普通の値であって、珍しくはない

と感じるべきである。

ではどれくらい小さいp値を見たら、帰無仮説を疑うほど、小さいことになるだろか。

珍しさは、「何度も繰り返してもめったに起きない」ことで決まるので、50万SNPの検定を何度も繰り返し、そのもっとも小さいp値がどのような分布をするかを調べてみればよい。

その分布に照らしてめったに起きないほど小さければ、そのとき初めて帰無仮説を疑えばよい。

```{r}
n.iter <- 1000
minp <- apply(matrix(runif(n.iter*n.snp),nrow=n.iter),1,min)
hist(minp)
```

図に示すように、大きい方に裾を引いた分布である。
```{r}
summary(minp)
quantile(minp,c(0.01,0.05,0.1))
```

50万個のSNPによるGWASを1000回実施したときの各GWASでの最小p値の分布である。
1000回のGWASのそれぞれは、全く遺伝因子が関与していないと思われるようなフェノタイプに対して実施してある。
小さい方から1パーセンタイル、5パーセンタイル、10パーセンタイルを示したが、この値に比べて、小さい最小p値を観測したときにはじめて、そのp値をもたらしたSNPとフェノタイプとの間では帰無仮説が成り立たないのではないかと考えることとなる。

この数値は大まかには$0.01/500000=2\times 10^{-8}$でとして算出できる。
珍しさ0.01を検定の個数50万で割った値である。

このようにして、珍しさの基準を決める方法をボンフェロニ補正法と呼ぶ。

最も小さいp値の大小を問題にしたが、すべてのp値についても情報を得ることにする。
たくさんのp値のすべてが帰無仮説に由来しているのかどうかを図で示すには、小さい順にならべて、$y=x$の直線に乗るかどうかでわかるのだった。

実際にやってみる。
横軸が0から1となるように工夫しておく。

```{r}
plot(ppoints(n.snp,a=0),sort(p),cex=0.01)
```

ただし、これでは、p値が小さい部分が、原点付近の狭い範囲に固まっていて、よくわからないので、その部分を拡大してみるために、対数をとれば
```{r}
plot(log10(ppoints(n.snp,a=0)),log10(sort(p)),cex=0.01)
```

また、場合によって(論文ではそちらの方が多いが)、p値ではなく、統計量($\chi^2$統計量)の観測値とその期待値とで、直線になっているかを示す場合もある。

```{r}
plot(qchisq(ppoints(n.snp,a=0),1),qchisq(sort(p),1),cex=0.01)
```

# 対立仮説と推定

## たくさんの治療法報告

2種類の治療薬の効果を比較した研究が次から次へと発表される。
小さなp値とともに有効だ、とする報告もあれば、あと一歩で帰無仮説が棄却できるレベルだったけれど、ちょっと足らなかったが、でも、惜しかった、という感じで報告されることもある。

単純に考えるためにp値が0.05であったら、必ず発表され、0.05から0.07の間のp値の場合には、10%が発表されるとする。

そんな報告をたくさん集めると、p値はどんな分布になるか。

```{r}
n.t <- 100000
p <- runif(n.t)
ok <- which(p<0.05)
sub.ok <-which(p>=0.05 & p<0.07)
published.sub.ok <- sample(sub.ok,length(sub.ok)*0.1)
hist(p[c(ok,published.sub.ok)])
```

0.05より小さい部分と0.05から0.07の部分とで段差はできるが、それぞれの範囲で一様分布になる。

## 検出力

実際の治療法の報告のp値の分布がこうなっているかというとそんなことはない。

カットオフの閾値0.05を下回ったあたりにたくさんの報告が集中する。

それはなぜかというと、

> 効果があるならこれくらい、という見込みの下、必要な標本数を計算して、スタディをデザインする

からである。

一般的な標本数の計算をしてみる。

治療法Aの有効率が50%、治療法Bのそれが60%であると見積もられるとき、p値の閾値として0.05で両者に違いがないという帰無仮説を棄却することにして、デザインする。
きちんとp値が0.05を下回る確率が80%となるには、検体数を2群に何人ずつとればよいかを考える。

```{r}
pw.out <- power.prop.test(p1 = .50, p2 = .60, power = .80)
pw.out$n
n. <- round(pw.out$n)
n.
```
で得られる数値が、2群に同数の検体数を取る場合に確保するべき検体数である。

実際、この条件で、このスタディを繰り返してp値がどうなるかを見てみる。
```{r}
n.iter <- 1000
ps <- rep(0,n.iter)
pa <- 0.5
pb <- 0.6
for(i in 1:n.iter){
  A <- sample(0:1,n.,replace=TRUE,prob=c(1-pa,pa))
  B <- sample(0:1,n.,replace=TRUE,prob=c(1-pb,pb))
  tbl <- matrix(c(table(A),table(B)),ncol=2)
  ps[i] <- chisq.test(tbl,correct=FALSE)$p.value
}
hist(ps)
```
実際p値が0.05を下回った割合は
```{r}
ok <- which(ps<0.05)
length(ok)/n.iter
```
と、80% 程度になっている。


では、このようなスタディがたくさんあって、p値が0.05を下回ったら報告され、
若干、上回ったときには、その10%が「惜しい」として報告されるとすると、p値はどんな分布になるのだろうか。

```{r}
sub.ok <- which(ps>=0.05 & ps<0.07)
published.sub.ok <- sample(sub.ok,length(sub.ok)*0.1)
hist(ps[c(ok,published.sub.ok)])
```

0.05ぎりぎりの報告は少なく、かなり小さいp値の報告が目白押しになるはずであることがわかる。

実際の報告のp値の分布は、もっと、0.05ぎりぎりのものが多い印象がないだろうか。

## 色々な治療効果のスタディの混成

ここで示した分布は、2群の成功率が0.5,0.6と「きちんとした差がある」場合のp値の分布であり、
実際の報告の集まりの中には、このように違いのはっきりしたスタディの報告もあれば、まったく2群に違いがないのにもかかわらず、確率0.05で生じてしまう小さなp値のせいで報告されるものも含まれること、また、若干の差はあるが、2群の成功率が0.50.55だったり、0.50.51だったりするものも含まれるからである。

では、A群の成功率は0.5であって、B群の成功率は0.5の場合と0.6の場合とが半々だとどういうことになるだろうか。

```{r}
n.iter <- 1000
ps <- rep(0,n.iter)
pa <- 0.5
pb <- 0.6
for(i in 1:n.iter){
  A <- sample(0:1,n.,replace=TRUE,prob=c(1-pa,pa))
  if(i < n.iter/2){
    pb. <- pa
  }else{
    pb. <- pb
  }
  B <- sample(0:1,n.,replace=TRUE,prob=c(1-pb.,pb.))
  tbl <- matrix(c(table(A),table(B)),ncol=2)
  ps[i] <- chisq.test(tbl,correct=FALSE)$p.value
}
hist(ps)
```

p値のうち半分は、上で描いた0寄りの分布であり、残りの半分は0から1の一様分布であって、それを併せた分布となる。

このようなスタディから、報告がなされるとき、p値の分布はどうなるだろうか。

今回は、p値が0.05の場合にのみ報告し、それ以外は厳密に報告しないものとする。
```{r}
ok <- which(ps<0.05)
ok.false <- which(ok<(n.iter/2))
ok.true <- which(!(ok<(n.iter/2)))
h <- hist(ps[ok],plot=FALSE)
hist(ps[ok[ok.true]],breaks=h$breaks,ylim=h$ylim,density=20)
hist(ps[ok[ok.false]],breaks=h$breaks,ylim=h$ylim,density=15,col=2,add=TRUE)
```

たくさんの黒(A,Bの2群に違いがあってp値が0.05より小さくなったスタディ)があり、少数の赤(A,Bの2群に違いがないのにp値が0.05より小さくなったスタディ)があることがわかる。

黒のバーはp値が小さい方で長いが赤のバーは高さが一様である。

そのことから、p値が小さいスタディがあれば、それは本当に2群間に違いがあることを反映している割合が高いし、p値が0.05に近ければ、本当は2群に違いがないのに報告されてしまっている可能性も低くないことがわかる。

では、2群に違いがある(有用な治療法である)と判断したもののうち、どれくらいの割合が、「本当に有用」なのだろうか。

```{r}
length(ok.true)/length(ok)
```

そして、本当は有用でないのに、有用であると誤って判断されたのはどれくらいの割合だろうか。
```{r}
length(ok.false)/length(ok)
```

## 感度・特異度・事前確率・Positive Predictive Value・Negative Predictive Value・False Discovery Rate

これが、診断のための検査であったとすると、どう考えればよいだろうか。

まず、2群に違いがあると仮定したスタディ、違いがないと仮定したスタディが50%ずつだったので、検査の前に真に「有効〜病気である」ことの事前確率は50%。

ついで、本当に「差がある〜病気である」場合に、p値が0.05未満になる確率は80%であるから、これは、検査で言えば、感度が80%。

また、本当は「差がない〜病気ではない」場合にp値が0.05未満になる確率は5%であるから、特異度は95%。

このとき、p値が0.05未満になったときに、本当に病気である確率が、上で計算した
```{r}
length(ok.true)/length(ok)
```
であり、これはPPVに相当し、

```{r}
length(ok.false)/length(ok)
```
は本当は病気でないのに病気であると判断されてしまった確率で、これはFalse Discovery Rate(FDR)である。

# FDR法 とq値

すべての検定が帰無仮説を満足することを前提にしてp値を補正する例をGWASで見た。

かなりの割合で対立仮説が成り立つ場合に、検定をすると、PPVが問題になり、FDR=1-PPVも気になることを、有用な治療法の研究を例として見た。

たくさんの検定を実施し、それらのうちのいくらかが対立仮説を満足するような状況の典型例は、2(疾患群vs.健康群、正常組織 vs.異常組織)の間での遺伝子発現量を全遺伝子について比較するトランスクリプトーム解析である。

2つの群は明らかに異なっており、その2群での遺伝子発現量が違っているのも確実であり、興味があるのは、「どの遺伝子がどれくらい」2群間で違うのか、である。
そのようなとき、すべての遺伝子に関する2群間の違いついて帰無仮説検定をし、そのすべての遺伝子に関して帰無仮説が成り立つことを想定したうえでマルチプルテスティング補正をするのは適当ではない。

ある閾値で違いがありそうな遺伝子を選び出したときに、選ばれた遺伝子にどれくらい本物があり、どれくらい偽物があるかの比率が、納得の行くものであることが大事であり、その比率がわかっていること、わかった上でその結果を解釈することが大事である。

このような場合に、FDRを制御して全遺伝子が「本当に2群間で違う」かどうかを評価する手法を「FDR法」と言う。

FDR法では、たくさんのp値が得られたときに、そのpの値をあるルールで変換して0からの1の値に変える。

その値を補正後p値とかFDR法によるq値と言う。

p値が小さい方が「本物」らしいわけだが、補正後p値も小さい方が「本物らしい。また、2つのp値があったら、その大小関係は、補正後p値でも同じである。

FDR法にはいくつかの方法があるが、まずはBH法と呼ばれる方法でやってみる。

p値と補正後p値との関係は以下のようになる。
値が変わっているので、一直線ではなくなっている。
また、pの値よりも補正後p値の方が大きくなっていることも見て取れる。
```{r}
p.corr <- p.adjust(ps,method="BH")
plot(ps,p.corr,cex=0.1)
```

このようにして計算した補正後p値の度数分布を補正前のそれと比べてみる。
全体的に少し大きな値に補正するということが度数分布ではどのように見えるかがわかる。

```{r}
hist(ps,density=20)
hist(p.corr,col=2,density=15,add=TRUE)
```

このような補正をした結果、PPVとFDR=1-PPVがいくつになったかを調べてみる。

```{r}
ok.corr <- which(p.corr<0.05)
ok.corr.false <- which(ok.corr<(n.iter/2))
ok.corr.true <- which(!(ok.corr<(n.iter/2)))
length(ok.corr.true)/length(ok.corr)
length(ok.corr.false)/length(ok.corr)
# 補正前の値は
length(ok.false)/length(ok)
```

FDRが補正により小さくなっていることがわかる。
そうは言うものの、FDRの値がきっちり0.05になるわけではないこともある。

では、きっちりFDRがコントロールできていないのならば、どんな基準で補正がなされているのだろうか。

すべて帰無仮説が成り立っているときにはどうなるかを見る。

```{r}
p.null <- runif(1000)
p.corr.null <- p.adjust(p.null,method="BH")
plot(p.null,p.corr.null)
hist(p.corr.null)

```

ほとんどすべての補正値が1付近になっており、補正後に0.05になるようなものはないことがわかる。
このことから、

>たくさんのp値が帰無仮説に由来していると考えられる場合には、そのすべてが「ありきたり〜補正後p値が1」となるような補正法であることがわかる。

考え方としてはこのような補正法であるが、さまざまな改良がなされている。
その基本とは、補正前のp値の分布が、すべて帰無仮説の場合と比べてどのようにずれているかを評価し、それに応じて、どのくらいの割合で「対立仮説が真である」かを見積もり、補正後のp値(q値)が、「偽の陽性率」とみなせるように補正するものである。

その一つを以下に示す。
```{r}
library(fdrtool)
fdr <- fdrtool(ps, statistic="pvalue",verbose=FALSE,plot=FALSE)
plot(ps,fdr$qval) # estimated Fdr values
hist(fdr$qval)
```

次世代シークエンサーを使ったデータ解析〜オミックス統計学入門2014

  • ひとまず、次世代データ、Deep sequencing データの資料はこんなもの
  • 「こうすれば、よい」という段階ではないので、コンセプト説明が重くなり、また、いざというときの汎用性重視になった…結局、「入門」としては役に立てにくいのだが…
  • まだ、1度も講義使用していないのでベータ版扱いです
  • 90分で1-2回分相当
  • Rmdファイルです。html化、epub化できます(やり方はこちら)
  • html化、epub化が面倒くさければ、kindleで1米ドルでも(こちら)

# 次世代シークエンサーによるDeep sequencing

# はじめに
次世代シークエンサーによるDeep sequencingの概要に触れた後、
そのデータ産出特性に合わせて行われる確率的データ解析とはどんなことかを簡単にまとめ、
それに引き続いて、ゲノム・レファレンス配列を用いてバリアント検出をすることを対象として、データ処理・データ解析をするプロセスを確認する。
そのプロセスを確認しながら、そのデータ解釈やデータ処理のどういうところに確率や尤度やモデルが使われるのかに関して触れることで、確率的データ解析とその解釈の勘所を理解することを目標とする。
実務とは少しずれた話を中心にするが、それは、この分野のデータ解析はまだ実験系自体が発展途上であり、それに合わせたデータ解析手法も変化しており、その細部の知識が役に立つ期間は長くないと予想されるからである。逆に、どうして確率的にデータを扱うのか、その背景の考え方は何かを押さえておくことは、Deep sequencingにも有用でありそれ以外のデータ解析にも重要だからである。

# 次世代シークエンサーによるDeep Sequencingの基礎事項

## ゲノム・エピゲノム・トランスクリプトーム
次世代シークエンサーは核酸配列を読む。
従って、ゲノム(DNA)、トランスクリプトーム(RNA)の両方が対象となり、DNAの修飾状態によって対象を絞る技術と組み合わせることでエピゲノムの解析もできる()<img src="DeepSequencing.jpg" />

## 短いリードが大量に
次世代シークエンサーで読まれるリードは長くする方向での進展もあるが、現時点では短くてもよいから大量に高速に読むことでアドバンテージを稼ぐ技術となっている。
大量に読むと、ある特定の塩基箇所・範囲について複数回読むことになる。その複数の読み重ねを「深さ Depth」と言うのでDeep sequencingと言う。

そのため、図に示すように、データを活用するにあたって小規模実験データと異なる手順が必要となる
### 生データ処理

バルクで出てくるデータなので、それをひとまとめで扱うための処理がある。

### データのクオリティ・チェックとリードごとの処理

> データは全体としてどのようなクオリティなのかを、クオリティ指標の代表値と分布に照らして評価する。

> 同様に1本1本のリードごとに、そのリードのクオリティの評価をする。

> さらに個々のリードの個々の塩基のコールについても評価する。

### アセンブルかマッピングか

ヒトの全ゲノム配列のような、ある特定の配列(レファレンス配列)を地図として、その上のどこかしらにリード配列を割り当てることをマッピングと呼び、リード配列同士をつなぎ合わせて長い配列(contig)にすることをアセンブルと言う。

レファレンス配列がないときには図の左端のようにリード配列をアセンブルし、レファレンス配列があるときには、それ以外の場合のようにリード配列をマッピングする。

ジグソーパズルで言えば、完成図がない状態でピースをつなぐのがアセンブル、完成図に照らして絵を作るのがマッピングである。

## Deep sequensingの苦労
### 短いことの苦労
リードが短いと配列の特異性が低いので、区別のつかないリードがたくさんあり過ぎて、アセンブルもマッピングも難しい。
ジグソーパズルで言えば似過ぎたピースが多いのと、出来上がる絵のあちこちに似たような部分が多いのとで、完成が難しいのと同じようなものである。

### リードの精度による苦労
たくさんのリードが得られるがそれらの11つが実験結果であり、きれいに読めているものとそうでないものとがある。高速に大量に読むので、個々のリードの精度には必然的にばらつきがある。不確か情報・偽情報もあるということである。
それらを用いることはジグソーパズルで言えばピースの絵柄に間違いがあったり汚れていて見えないところがあったりすること・ピースの輪郭の作り粗雑でピース同士(リード同市)をつなぐかつながないかの判断が難しかったりすること、に相当する。
これはアセンブル・マッピングの段階のことであるが、なんとか作り上げたアセンブリ・マッピング後情報を読み取る際にも苦労がある。同じくパズルに照らせば、2つのサンプルから作り上げた地図の違いを見つけるのは、間違いさがしであり、出来たパズルの絵に不確かさや汚れがあればその作業は難しい。
また、Deep sequencing では、リードというピースを同じ場所に積み上げ、その高さを場所ごとに比較したり(RNA-seqやChIP-seqの定量評価)、積み上がったところのアレルの比率を検討したり(バリアント検出、頻度推定)する作業も、同じく難しくなる。

### 場所ごとの読まれる回数のばらつきによる苦労
サンプルのどこがどれくらい読まれるかは、読まれる配列がどれくらいあるか、その配列がどれくらい読みやすいか、という要素に左右され、ばらつく。また、実験過程で配列を指数関数的に増幅する過程(PCR)を使う場合には、量のばらつきが増幅される。

元のサンプルがゲノム配列であるときは、2コピーずつあることが原則である。
この場合には、2コピーが均等に読まれていないことを考慮してバリアントの有無・ジェノタイプの推定をする必要が出る。また、部位ごとに得られる情報の量に違いが出るため、部位ごとにバリアント検出・アレル頻度推定の感度や特異度がばらつくという問題が生じる。
一方、読んだリードのデプスによって元のサンプルの量を定量的に評価したい場合(RNA-seq,ChIP-seq,プールしたゲノムを読んでアレル頻度定量をする)には、読まれる回数にばらつきがあることは、定量結果に直結する問題である。

# 確率的データ解析とは
## Deep sequencing における確率的データ解析の枠組み
以上のように次世代シークエンサーによるDeep sequencingでは複数の苦労がある。
これらはいずれも

> 不確かさのある情報が得られる

> 大量に得られるので分布が取れる

> 不確かさを分布によって補う

> 補うにあたって確率的判断をする

という枠組みで対処しようとしている。

## 確率・尤度、事前確率・事後確率
この確率的データ解析とは次のようにまとめられる。

> 知らないことがあって、それを知りたい(ある場所にバリアントがあるか、ある場所のアレル頻度はいくつか、ある座位のトランスクリプトはどれくらいあるか、など)

> 知らないことがもし、こうなっていたらと仮定する(バリアントがなかったら、バリアントがあったら、ある座位のトランスクリプトが別の座位に比べて●倍だったら、など)

> 仮定があったときに、ある結果がどれくらいの確率で生じるかを計算する

> 実際に得られた結果に照らして、どの仮定ではどれくらいの確率だったかを確認する。この確率を尤度と言う

> 仮定ごとに尤度の大小が出るので、それに応じて、どの仮定を信じるのがよいのかを判断する

これが、観察データに基づいて仮説を確率・尤度的に判断するプロセスである。

また、この確率・尤度の計算にあたって、色々と条件を想定する。
どのような配列は読まれやすいか(GC含量は読まれやすさに影響する要素)、読み間違いエラー率はどれくらいか・リードの個々の塩基はどれくらい信用できるか、などがそれにあたる。

このような条件を加味するときには、どのような値で取り込むかが問題になる。
ここで用いる条件の正確な値は知りえないからである。
条件に用いる値が変われば、それに基づいて計算する確率・尤度ももちろん変わり、最終結果ももちろん変わる。
従って、確率的データ解析の結果は、その真贋についても確率的に判断しなくてはならない。

また、大量のデータを取ってその分布を利用すると書いた。
この分布の利用にあたって、その分布がどのような形をしているのかにモデルとなる分布式を使うかもしれない。
モデルを入れると真実がそれとどれくらい合ってるのか、という点で、データ解析結果は影響を受ける。
モデルが複数考えられる場合には、どのモデルを使うかで結果が変わってくる。
従って、Deep sequencingデータの解析結果は、モデルによって結果が変わることも念頭に置いて解釈しなくてはならない。

# リシークエンシングによるバリアント検出

1人のDNAをレファレンス配列にマッピングし、バリアント箇所を検出する場合を考える。
検出結果の解釈に影響を与える点を中心にその過程をたどる。

## サンプル調整

読む個人を決めてDNAを採取・調整する。
全ゲノムにするか・エクソン部分だけにするか(エクソーム・シークエンシング)の選択肢があるが、

> エクソン部分だけにする場合には、その部分だけを取り出すというプロセスが入り、どのセグメントのどちらのアレル(母方・父方の2アレル)がどれだけsequencingに回るかに関してばらつきが入る。

> セグメントごとの取り出され方の差は、セグメントごとのDepthに影響が及ぶ。Depthが小さいセグメントは情報量が少なくなり、Depthが大きければ情報量が増えるので、セグメントごとにバリアントを見つけられるかどうかのパワーが異なり、バリアントであると判断した場合の信憑性が変わる

> バリアントは同一座位に1つのアレルがあるのか2つのアレルがあるのかを判断することで検出される。2アレルがサンプル調整の段階でどれだけ取り込まれるかにばらつきが出るので、そのばらつきがバリアントの検出感度に影響する。50:50で取り込まれた場合が最も感度が高くなる

> とりこみに差が出るのは、次のように分類できる

>> 片方がもう片方よりたまたま多く取り込まれることによって生じる

>> 片方が実験原理的に多く取り込まれる理由があることによって生じる

>>> バリアントのアレル自体が取り込まれに効率に影響を与える場合

>>> バリアントの周辺の配列差(近傍バリアント)が取り込みまれ効率に影響を与える場合


## マッピング前の評価

準備したサンプルをDeep sequencingする。

実験全体を評価する。

実験は生ものなので、その実験全体がうまく行っているのかどうかの評価が必要である。
その実験手法で「成功」とされる標準的な値、との範囲が量と質についてあるはずなので、それと比較する。

また、複数の実験をしているときには、実験単位で量的評価・質的評価をして、その異同を確認し、突出しているものについてはその原因を確認する。

### 量的評価

どれくらい読めたかの総量に関する評価が必要である。

> 全部で何リードか $n.read$

> それはリード長を加味して何塩基対分か $\sum_{i=1}^{n.read} L_i$ (ただし$L_i$は第i番リードの長さ)

> 読む対象の全長(ゲノム全長はいくつか、エクソームなら対象の全長はいくつか) $L_{total}$に対して、平均して何回ずつ読まれているか $\frac{\sum_{i=1}^{n.read} L_i}{L_{total}}$

最後の指標はある一定以上のDepthがあって初めてバリアントの検出が可能であることから、大まかな指標としてわかりやすい。
Depthとバリアントの検出とに関する確率的実習は、後述の『確率的データハンドリング:Depthとバリアントの検出』を参照。




### 質的評価

リード属性によって実験全体の質的評価をする

リードの属性としては、量的属性として長さ、質的属性として塩基単位のクオリティ値がある。
多数のリードについて、その量的・質的属性の評価をする。

また、リードの属性には、読まれ方・データのクオリティに影響する因子として、
リード内の相対的位置(読み始め・読み終わり・その中間)と、塩基の構成比(GC割合)などがあるので、
その条件での評価も有用である。

> リード長の評価をする

> リードの各塩基のクオリティ値の評価をする

> リードの相対的位置(読み始め・読み終わり・その中間)とクオリティ値には関係があるので、リードの相対的位置ごとのクオリティ値の評価をする

> リードのGC割合を条件として、リードのクオリティ代表値を評価する。

なお、属性の評価をするとは、その分布をとることと、代表値を確認することである。

その分布・代表値は、実験手法の標準的なそれと比較する。
また、自前の実験群で同様に比較し、突出があれば原因を確認する。

塩基別のクオリティは、実験上の塩基が実際の塩基と一致しているかどうかを表している。
つまり、リードの塩基は誤っているかもしれないわけである。
この誤った塩基を持ったリードがマッピングされれば、そこに本来はないはずのバリアントが存在することになる。
従って、クオリティ値が低めなリードのデータを採用すると、バリアントの検出にあたってFalse positiveを生じることになる。

この問題の解決法としては、間違いの可能性が非常に低い(全くない、ということはない、という立場に立つのが確率的データ解析である)リードのみ、そのような塩基のみを解析に用いるという方法がある。
しかしながら、これは、活用する総リード数・総塩基数を減らすことであり、Depthを小さくする。
Depthが小さくなればFalse negativeが増えるわけで、むやみに不採用のリード・塩基を増やすことは適当ではない。

コールのエラー率とFalse positiveの関係については、後述の『確率的データハンドリング:コールエラーによるFalse positive』を参照


## マッピングの評価

ここまでは、実験全体、リードの集まり全体に関する評価について述べてきた。
リードをレファレンス配列にマッピングするステップを考える。

読まれるべきリードがあまり読まれなければ、完璧にマッピングが成功してもDepthは読まれた総本数以上にはならないし、読まれたリードのクオリティが完璧でなければ、マッピング自体がうまく行っても、リードの配列自体に偽情報があるわけであるから、やはりマッピングはうまくいかない。後者の場合は、

> マッピングできる位置が見つからない

> 間違ったところにマッピングされる

> マッピングが複数の場所に対応づく

という場合が多くなる、という形で現れる。

この節では、それらについては、いったん横に置き、リード自体は完璧であるが、マッピングに支障が生じる場合について検討する。

本数とクオリティはそれに重畳する形で影響を与えると考えればよい。

### レファレンス配列との不一致

レファレンスゲノム配列は、ある代表的な配列である。
従って、ある特定の個人のゲノム配列が、それと一致しないことは珍しいことではない。
逆に、バリアント検出のためにDeep sequencingをしている場合には、一致しない箇所を探しているのであるから、一致しない〜マッピングに困難がある〜場合にこそ興味があるわけである。

問題は、レファレンス配列と一致しない座位を有するリードと一致するリードに比べてマッピングされないかもしれないことと、もし、そうであるとすると、それはどういう場合であって、どの程度、どういう形で影響を及ぼすかである。

レファレンス配列とサンプル配列との違い

> リード上にバリアントが一つあれば、ない場合に比べて、リードがレファレンス配列とマッチするスコアは下がる

> 近傍にバリアントがあれば、さらに下がり、場合によっては、レファレンス配列にマッチしたとのスコアはさらに下がる

> ゲノム配列には民族特有なバリアントがあるので、サンプルを提供した個人とレファレンス配列の提供者とに民族的な違いが大きければ、リード配列とレファレンス配列との不一致確率はさらに上昇する

> ゲノム上には似た配列がある。いわゆるリピート配列もそれに相当し、コピー数多型もこれに該当する。また、pseudogeneのような相同性の高い配列もある。その場合には、より複雑な状況が出現する

### バリアントを含むリードのマッピング

#### バリアントの密度、ハプロタイプ

ある部分が長さ$N$のリードとして読まれるときに、レファレンス配列と同じ配列を読めば、完璧なリードの場合には$n_A=N$塩基が一致する。
$k$個の1塩基バリアントが含まれれば、$n_B=N-k$塩基が一致する。

単純に、$t$箇所までの不一致は許容し、$t+1$箇所以上の不一致は許容しないとすれば、リードに含まれるほどの近傍にバリアントが$t-1$個までしかないバリアントは、100%マッピングされ、リードに含まれるほどの近傍に$t$個以上のバリアントがあって、自身を含めて$t+1$箇所以上が不一致となるバリアントは、まったくマッピングされない。

このような単純な割り切り方をする場合には、一定以上にバリアント密度が高い部位では、バリアントの検出感度が下がることがわかる。

また、この「近傍にバリアントがある」というときには、レファレンス塩基と異なる塩基が同一DNA分子上に並んで、複数のバリアントのハプロタイプを作っていることを指している。

従って、個人のゲノムにおいて、バリアントが同じように近傍に存在しているとしても、母方・父方で異なるDNA分子上にあるバリアントの場合より、バリアントが同一のDNA上にあってハプロタイプを形成している場合の方が、マッピングし損ねる率が上がる。

バリアントのゲノム上の密度と、ハプロタイプ構成とがマッピングに影響を与える。

バリアントの密度が高くなったり、民族性を反映して、そもそもレファレンス配列がマッピング先を含んでいない場合には、レファレンス配列を複数にしたり、レファレンス配列を三省せずにリード同士をつなぐ(アセンブル)の要素を取り込むなどを考慮する必要が生じる。

複数のレファレンス配列を用いる場合には、人工的に配列を重複させている側面もあるので、そのことに配慮したマッピング判定が必要となり、その判定は確率・尤度的な評価となる。
また、アセンブルを加味した場合には、それによって見出されるバリアントの感度と特異度は、アセンブルを加味しないそれとどのような関係にあるのかも、確率・尤度的に評価する必要がある。

#### リードのクオリティ

ここでは、簡単のために、塩基の一致数・不一致数による単純なマッピングを想定して述べたが、塩基コールのエラー率を加味してマッピングする場合にも、影響は及ぶ。

一致数・不一致数に関する閾値を設定する場合には、アレル別のリード数がバリアントで少なくなる傾向がある、という形で現れ、そのリード数は整数である。

他方、エラー率を加味すると、個々のリードがどこにどれくらいの確率でマッピングされるかの数字が細かい数字になる。
この確率は、エラー率を加味しても、バリアントを有する配列に由来するリードは、バリアントを有さない配列に由来するリードよりも低くなる。

マッピングにあたって、この確率を閾値に照らして、マッピングされるか否かを分けるのであれば、アレル別のリード数はやはり整数となり、バリアントを有する配列に由来するリード本数は少なくなる。
また、閾値を定めずに、マッピングされるかどうかの確率を用いて、重みづけ加算をするとしても、やはり、バリアントを有する配列に由来するリードの重みづけ加算値は小さくなる。

これは、バリアントの検出感度に影響する。
そして、その感度への影響は、個々の箇所についてはたいしたものではないが、多くの箇所をゲノムワイドに対象とする場合には、影響が無視しえない可能性が出てくる。

#### リードのクオリティとバリアントの存在

バリアント塩基を完璧なクオリティで読んだとき、そのリード配列とレファレンス配列との違いは、1塩基分になる。
同じ部分を読んだリードのバリアント塩基のクオリティが完璧でないときには、そのリードとレファレンス配列との違いは、1塩基分より小さいとみなせる。
単純に言うと、レファレンス配列との違いが小さい方が、マッピングの成功度は高いのであるから、クオリティが完璧でないリードの方が、マッピング成績が良くなる。

近傍にバリアントがあるときも同様である。
あるバリアントは完璧なクオリティで読まれているとする。
その近傍にあってハプロタイプを構成しているバリアントも完璧なクオリティで読まれれば、そのリードはかなりレファレンス配列から遠くなる。
一方、近傍バリアントの塩基クオリティが完璧ではないときには、かえってマッピング成績が上がるから、ここでも、クオリティが完璧でないリードの方が、マッピング成績が良くなる。

このように、「そもそも違っている」何かを検出するときに、「違いに曖昧さが入っている」方が、好成績でマッピングされたり、「曖昧さを持つ」リードを取り込んだ方が、検出されやすくなったりする可能性があることには注意が必要である。

クオリティ等の情報をすべて取り込んで、すべての可能性を確率的に扱えば、この問題のある部分は解消するかもしれないが、どこかしらのステップで閾値設定をすれば、その段差の影響が発生し、感度が下がったり上がったりすることになる。

また、すべての情報を取り込む、と簡単に書いたが、情報を取り込んだ上で、スコアリングするにあたっては、事前確率設定やモデル設定が影響を及ぼしてくるから、やはり、完璧ではない。

ここから言えることは、素朴な方法で臨んでも、工夫を凝らして情報を駆使しても、いずれにも何がしかの曖昧さや捉えきれない幅、false positive, false negativeが発生するということである。

その曖昧さの程度・大雑把さの程度についての感覚を得るには、実験系・処理フローとそのオプションが変わるごとに何がしかのシミュレーションをしてみることも役に立つ。


### リピートの存在、pseudogeneの存在

レファレンス配列には、多くの似通った配列が含まれる。
非コーディング遺伝子領域のリピート配列もそれに含まれるし、
遺伝子とそのpseudogeneもそうである。
また、複数のコーディング遺伝子が良く似た配列を共有している場合もある。
似通った配列が複数箇所にある場合の、マッピングの一般的な問題は、マッピング箇所が一意に決めにくいということである。

バリアントが存在する場合には、事情が込み入ってくる。

二箇所に似通った配列があるとする。

今、あるリードがこの二箇所のどちらか片方を読んだものだとすると、マッピングではどちらにもマッピングされうる。
二箇所の配列が全く同一であれば、優劣は付かない。
二箇所の配列が異なっていれば、優劣が付く。

今、二箇所のレファレンス配列は同一であるとする。
片方にだけバリアントがあるとすると、
バリアントのない配列を読んだリードは、二箇所に同確率でマッピングされる。
バリアントを含む配列を読んだリードの、二箇所へのマッピングスコアは、バリアントを含まない場合よりも低くなるが、二箇所で同一の値となる。

二箇所のレファレンス配列が1塩基だけ違っている場合には、ぞれぞれの座位から読まれたリードは、それぞれのレファレンス配列に最も高スコアとなるが、もう片方の座位にも、一塩基バリアントがある場合と同じ程度にはマッピングスコアが出る。

バリアント探索をしているときには、このようなマッピングの良さは捨てたくないレベルであるから、気がもめる。

わずかな差であっても一致の良い方に迷わずマッピングをするのであれば、紛れることはないが、リードのクオリティが下がると、二箇所へのマッピングの良さの差が縮まるから、問題はそれほど簡単ではない。

さらには、二箇所のレファレンス配列に一塩基の違いがあり、その塩基箇所について片方の座位にバリアントがあり、そのバリアントのアレルが、もう片方の座位のレファレンス配列と一致している場合には、バリアント配列からのリードは、本来の座位とは別の座位にマッピングされることになる。

また、コーディング遺伝子とそれに似たpsuedogeneとがあった場合で、両方が増幅されリードを生んでいるにも関わらず、マッピング先をコーディング遺伝子に限定したとすると、pseudogeneに由来するリードはマッピング先がないので、次善のマッピング先として、コーディング遺伝子にマッピングされることになるが、これによって見出されるリード配列とレファレンス配列の不一致はバリアントではない。

したがって、バリアント探索においては

> リード配列とレファレンス配列の異同だけではなく

> レファレンス配列上にありえる、複数の似た配列のセットと、リード配列との異同とを、併せて比較し、

> リードがバリアント配列を反映していたら、どのようなマッピングになるかを考慮する

必要がある。

## マッピング後の評価

### バリアント検出

バリアントの検出とは、サンプルのアレル別配列数が(レファレンスアレル):(バリアントアレル)の比として、$2:0$$1:1$$0:2$かのいずれであるかを識別することである。

以下の『確率的データハンドリング』に、諸条件がその検出にどのような影響を与えるかを述べた。

プールサンプルでのバリアント検出では、$M$人の場合にアレル別配列数が$2M:0,(2M-1):1,(2M-2):2,...,M:M,...,0:2M$を識別することになる。

### 定量比較

#### 同一サンプル内での、座位上の比率の定量

同一座位にマッピングされたリードの情報から、何かしらの定量をすることは、バリアントのジェノタイプ推定・プールサンプルのアレル頻度推定と基本的には同じである。
ジェノタイプのバリアントの方が整数比であることを前提に推定しているのに対して、こちらの定量は連続値であると考える点が異なる。

#### 同一サンプル内での座位間の定量比較

座位別の定量比較をする場合には、素朴にはマッピングされたリード数を比較する。

座位ごとに読まれやすさ$\beta$に違いあることを想定した上で、座位間の定量比較をするのであれば、$\beta$を用いたモデルに照らして、サンプル中のテンプレート量の推定をし、その量を比較することもできる。

また、2座位のテンプレート量の比を問題にするのであれば、同様に座位ごとの読まれやすさを加味した上で、比を推定することができる。

#### 異なるサンプル間での定量比較

絶対量の比較か、(標準配列に対する)相対量の比較かで推定するべき値は変わるが、基本的には、実験が異なるので、両者が同じ土俵に乗るように、標準化するか、総量確認をした上での比較となる。

# 確率的データハンドリング
## Depthとバリアントの検出~理想的な条件~

ここで、Depthがどのようにバリアント検出に影響を与えるかに関する確率的評価をしてみる。

今、サンプルに2つのアレルA,Bが$n_A,n_B$本ずつあるとする。
ゲノムの場合には理想的には$n_A=n_B$である。

この$n_A,n_B$本がそれぞれ、平均Depth $\beta_A,\beta_B$になるような読まれ方をしたとする。
そのときに、実験データとしてA,Bのアレルが$N_A,N_B$本のリードとして得られる確率について考える。

ある確率でぽつりぽつりと起きるような現象があって、ある時間内に平均$\beta$回、起きるとき、実際に$0,1,2,...$回起きる確率はポアソン分布に従うということが知られている。

実際Deep sequencingをして、座位ごとのDepthの分布をとると、ポアソン分布とみなしても悪くなさそうな分布が得られる。

```{r}
beta_A <- 3
N_A <- 0:30
x_A <- dpois(N_A,beta_A)
plot(N_A,x_A,type="h")
```
平均して$beta_A=3$回だとすると、1度も起きない確率もそれなりにあることも含めて、図のように示される。
今、ヘテロ接合体であるとして、$beta_A=beta_B$であると仮定する。
この仮定はどういうことかと言うと、個人のゲノムDNAの段階でA,Bの2つのアレルが等量ずつあり、そのDNA抽出・サンプル調整の過程で、2つのアレルがサンプルに等量ずつ入ることに成功し、さらに、そこから等しい効率でリードが読まれるという状況を設定している。

このとき、2アレル併せて$N=N_A+N_B$本のリードが得られたとする。
そのときの$(N_A,N_B)$の内訳は$(N_A,N_B)={(0,N),(1,N-1),...,(N,0)}$$N+1$通りある。
このそれぞれがどのくらいの確率で起きるかを計算してみる。
そうすると、本当はヘテロ接合体であるのに、観察リードからはホモ接合体のように見える確率が計算できる。

```{r}
N <- 4
# アレルBはアレルAと同条件にする
beta_B <- beta_A
N_B <- N_A
x_B <- dpois(N_B,beta_B)
# A,Bが0:N,N:0で読まれる確率を積で計算
prob.N <- x_A[1:(N+1)]*x_B[(N+1):1]
# これが全体なので、和が1となるように補正
prob.N <- prob.N/sum(prob.N)
plot(0:N,prob.N,type="h",ylim=c(0,max(prob.N)))
# A,Bどちらかのリードしかない確率
prob.N[1]+prob.N[N+1]
```

では、この理想的な条件で総本数(Depth)が変わると誤ってホモ接合体とみなしてしまう確率が計算できる。

```{r}
Ns <- 1:50
prob.Ns <- rep(0,length(Ns))
N_A <- N_B <- 0:max(Ns)
x_A <- dpois(N_A,beta_A)
x_B <- dpois(N_B,beta_B)
for(i in 1:length(Ns)){
  N <- Ns[i]
  prob.N <- x_A[1:(N+1)]*x_B[(N+1):1]
  prob.N <- prob.N/sum(prob.N)
  prob.Ns[i] <- prob.N[1]+prob.N[N+1]
}
plot(Ns,prob.Ns,type="l")
# 確率を提示
round(prob.Ns,5)
```
$N=5$程度で、誤ってホモ接合体とみなす確率は0.05を切っていますから、Depth=5程度で信頼ができるとも考えられます。
しかしながら、ゲノム全体であれば30億塩基対($3\times 10^9$)あり、その100〜数100塩基に1箇所のバリアントがあるとすれば、$10^7$オーダーのバリアント箇所があり、それをエクソン部分(ゲノム全体の1%程度)に限定しても、まだ$10^5$箇所のバリアントが検出されるべきです。5%を見つけそこなうとすると、見つけそこなうバリアント数は$10^3$オーダーになります。
この見つけそこね数を3桁減らすには、「誤ってホモ接合体とみなす確率」を3桁小さくする必要が出ます。そうすると
```{r}
# 10^t で表示
round(log10(prob.Ns),2)
```
これで$-5$程度の値がほしいことになり、それは$N=19$程度であることがわかります。

実際、「このくらいのDepthがほしい」と思って実験計画をしても、実際にはDepthの小さい座位、大きい座位が生じる。
したがって、得られたDepthに応じて、
バリアントを見つけそこねている確率(False Negative)の率の概数を織り込みながら、データを解釈することが必要となる。

## Depthとバリアントの検出~理想的でない条件~

上のシミュレーションでは、ヘテロ接合体ならば、サンプル中の2アレルの割合が等しく、そのリード産生効率も等しいと仮定した。
また、実験エラーはないものとした。

サンプル中のアレル比、そのリード産生効率の不均等は、ヘテロ接合体であるのに、ホモ接合体であるように見えるデータを生む、と言う意味で、感度を下げる働きをする。

逆に、実験エラーがゼロではないことを想定すると(実際、ゼロではない)、バリアントはないのに、バリアントであるかのような実験結果となる(False positive)。

この2つのタイプの「非理想的条件」について、実習する。

なお、ここでは、細かいことよりも、条件によって結果の解釈がどのように変わるかを実感することを主目的とする。

### 不均等

サンプル調整段階でアレル別のDNA分子数の不均等とDeep sequencingにおけるリード産生の不均等に分けられる。

2種類のものから構成されている集まりを考える。
それがたくさんちょうど同数あるところとする。
そこから、適当に抜き出すとき、抜き出したものはおよそ同数であり、抜き出す数が多ければ多いほどその比率は半々に近くなる。
逆に言えば、少ないサンプルでは比率がばらつく。

大容量の実験系・十分な濃度の標本を用いた実験系・よく混ぜ合わせてある実験系では、偏りが生じにくく、比率が半々になりやすいが、その反対の条件では半々からずれる。

Deep sequeincingはさまざまなところでミクロ化しているので、偏りを生じる要素は少なくない。

さらに、それを増幅する過程があることで偏りは大きくなる。

ここでは、どのような過程がどの程度の偏りを生じるかについては触れず、偏りが生じた条件でDepthとFalse negative(ヘテロ接合体なのにホモ接合体のようなデータが得られること)とがどのような関係にあるかを確認する。

### 不均等状態でのDepthとFalse negativeの関係

不均等状態でA,Bのアレルから平均$\beta_A \ne \beta_B$本が観察されるような状況になったとする。
このとき、ヘテロ接合体であるのに観察リードのすべてが同一アレルになる確率(False Negative率)を計算してみる。
上の例で$\beta_A=\beta_B$の場合をやってあるから、$\beta_B = \beta_A \times (1+\alpha); \alpha =0,0.5,1,1.5,2,...,10$として計算する。

```{r}
Ns <- 1:50
beta_A <- 3
N_A <- N_B <- 0:max(Ns)
x_A <- dpois(N_A,beta_A)
alpha <- seq(from=0,to=10,by=0.5)
prob.Ns <- matrix(0,length(alpha),length(Ns))
for(j in 1:length(alpha)){
  beta_B <- beta_A*(1+alpha[j])
  x_B <- dpois(N_B,beta_B)
  for(i in 1:length(Ns)){
   N <- Ns[i]
   prob.N <- x_A[1:(N+1)]*x_B[(N+1):1]
   prob.N <- prob.N/sum(prob.N)
   prob.Ns[j,i] <- prob.N[1]+prob.N[N+1]
  }
}

matplot(Ns,t(prob.Ns),type="l",xlab="Depth")
```

両者の比が大きくなるほど、Depthが大きくないとFalse Positiveになることがわかる。
ここで、両者の比を最大で10としたが、これは、指数関数的増加プロセスが何かの拍子で作用すれば、生じるのにそれほど難しい数字ではない。

また、プールサンプルにおけるアレル頻度推定をする場合には、A,Bアレルが10% vs. 90 % などは、見出したいと考えているだろう。このときのA,Bアレルの比は10である。
これは5人をプールして、10本の染色体のうち1本がバリアントである場合である。
このことに鑑み、より大きなA,Bアレル比、$\beta_B = \beta_A*2^{\alpha};\alpha = 0,1,2,...,6$も実施してみる。これは
```{r}
print(paste("アレル頻度比",2^(0:6),"倍"))
      
print(paste("アレル頻度",1/(1+2^(0:6))))
```
のことである。
```{r}
Ns <- 2^(0:10)
beta_A <- 3
N_A <- N_B <- 0:max(Ns)
x_A <- dpois(N_A,beta_A)
alpha <- 0:6
prob.Ns <- matrix(0,length(alpha),length(Ns))
for(j in 1:length(alpha)){
  beta_B <- beta_A*2^alpha[j]
  x_B <- dpois(N_B,beta_B)
  for(i in 1:length(Ns)){
   N <- Ns[i]
   prob.N <- x_A[1:(N+1)]*x_B[(N+1):1]
   prob.N <- prob.N/sum(prob.N)
   prob.Ns[j,i] <- prob.N[1]+prob.N[N+1]
  }
}

matplot(Ns,t(prob.Ns),type="l",xlab="Depth")
```

### ポアソン分布仮定よりも負の二項分布仮定
以下は、番外編。
部位によらず一定のパラメタで表されるポアソン分布に従うと考えてきたが、部位によって$\beta$の値がばらつくのが実情である。

場所によって$\beta$がばらつくと、どこでも$\beta$が一定であると仮定したときと比べて、Depthの分布のばらつきが大きくなる。

実際、部位ごとに$\beta$の値をいくつにするのかについて、$\beta$値に影響を与える因子で補正してやると、
部位ごとに$\beta$の値のより真実に近いと思われる値が推定できることが確かめられている。

その$\beta$を基にすると、実際のDeep sequencingデータでの、部位別のDepthをポアソン分布とみなすモデルのあてはまりは、$\beta$値をゲノム全体で一様にした場合よりも、実データのあてはまりが良くなる。
しかしながら、まだあてはまりは十分ではない。それをさらに改善するために、負の二項分布を用いることもできる。

これはどういうことかというと、平均してDepthが$\beta$になるような座位で、実際に何本が読まれるかが、ポアソン分布よりもばらつくことを意味している。読まれる本数がばらつくということは、たまたまあるアレルだけ読まれない、という確率が高くなる、ということである。

従って、現実には、ポアソン分布を仮定するよりも、False negativeの起きる確率は高く、その高いFalse negative率を推定するには、負の二項分布モデルを使った計算がより正確だ、ということである。

以下は、負の二項分布をモデルとした際の計算という細部である。

その場合には、さらに、False neagtive率が上がることをシミュレーションしたのが下記である。

ここでのポイントは、False negativeをもたらしている要因を取り込み切らなければ、汲み取っていないFalose positive率があるから、宝探しをしているときには、その見落としの可能性に注意を払い、どうしても確実に見つけたい部分があるならば、そこについては、単純なモデルのFalse negative率で安心せずにDepthを上げる必要がある、ということである。
```{r}
my.dnegbinom <- function(x,k,r){
  choose(k+r-1,k)*(1-x)^r*x^k
}
r <- 20
p <- r*beta_A/(1+r*beta_A)
x_A.NB <- my.dnegbinom(p,N_A,r)
x_B.NB <- my.dnegbinom(p,N_B,r)
prob.N.NB <- x_A.NB[1:(N+1)]*x_B.NB[(N+1):1]
# これが全体なので、和が1となるように補正
prob.N.NB <- prob.N.NB/sum(prob.N.NB)
plot(0:N,prob.N.NB,type="l",ylim=c(0,max(prob.N.NB)))
# A,Bどちらかのリードしかない確率
prob.N.NB[1]+prob.N.NB[N+1] 
```

### コールエラーによるFalse positive

実験データであるから、不確かなデータ、誤った塩基コールは必ずある。
クオリティ値は、塩基ごとに、どのくらいのエラーの可能性があるかについて、数値として提供している。

今、ある座位についてホモ接合体のサンプルを読み$N$リードが得られたとすれば、コールエラー率が$e=0,10^{\alpha};\alpha=-5,-4,-3,-2,-1$のとき、$N$本のすべてにエラーが入らない確率は$1-(1-e)^N$である。


```{r}
# Depth
Ns <- 1:100
# エラー率
es <- c(0,10^((-5):(-1)))
X <- matrix(0,length(Ns),length(es))
for(i in 1:length(Ns)){
  X[i, ] <- 1-(1-es)^Ns[i]
}
matplot(X,type="l")
```

では、エラーコール率が$e=0.01$で一定のとき、読まれる総本数が$N=10$のとき、エラーリードの数が$Ne=0,1,2,...$の確率はいくつかを確認する。

```{r}
e <- 0.01
N <- 10
X <- dbinom(0:N,N,e)
plot(0:N,X,type="h")
```

今、Ne=0の確率は
```{r}
X[1]
```
と、確実と言うには頼りないが、
今、Ne=0,1の場合の確率を合算すると
```{r}
sum(X[1:2])
```
とかなり高いので、$Ne=0,1$の場合ならば、ホモ接合とみなしてもよいと考えるかもしれない。

では$e$が増えるとどうなるだろうか。

```{r}
e <- 0.05
N <- 10
X <- dbinom(0:N,N,e)
plot(0:N,X,type="h")
cumsum(X)
```

示した数字の4番目か5番目当たり($Ne=3,4$)までは、ホモ接合体でも起きうると考えないといけないかもしれない。

特に、ゲノム上の大多数の塩基はホモ接合なので、Deep sequencingで大量の塩基対($10^6$$10^9$)を読むきには、もっと厳しい値$0.99999999$くらいまではホモ接合とみなしておかないとFalse positiveの確率は低くても、False positiveの件数は馬鹿にならなくなる。

$e=0.01$に戻して、$N$を増やしてみよう

```{r}
e <- 0.01
N <- 50
X <- dbinom(0:N,N,e)
plot(0:N,X,type="h")
cumsum(X)
```

$N=50$リードのうち、5,6リードはエラーとして許容しないといけないことがわかる。

### False negative とFalse positiveとのバランスを取る

バリアント検出にあたっては、若干のエラーを許容し、それに伴って、ある割合のリード数がレファレンスと異なっていても、ホモ接合とみなさなくていけない、という要請と、ヘテロ接合であっても、リードの2アレル比は半々からずれるが、どれくらいまでの偏りはヘテロとみなすか、という要請とのせめぎ合いとなる。

その様子を確認する。

```{r}
e <- 0.01
N <- 10
X.homo <- dbinom(0:N,N,e)
beta_A <- beta_B <- N/2
N_A <- N_B <- 0:N
# ヘテロについては理想的なポアソン分布モデルを使用
x_A <- dpois(N_A,beta_A)
x_B <- dpois(N_B,beta_B)
X.hetero <- x_A*x_B[(N+1):1]
X.hetero <- X.hetero/sum(X.hetero)
matplot(0:N,cbind(X.homo,-X.hetero),type="h",main="e=0.01,N=10")
```

上向きの黒いバーはホモ接合のときにバリアントに見えるリード本数が0,1,2,...,N本、観察される確率であり、
下向きの赤い破線バーはヘテロ接合のときにあるアレルが観察される本数が0,1,2,...,N本である確率である。

$N=10$本のリードが観察されるときに、少ないアレルのリード本数が1本のとき、2本のときには、ホモ接合・ヘテロ接合のどちらにしても、目に見える程度のFalse positive, False negativeがあることがわかる。

Nの値が大きくなれば、2つの山の分離は良くなる。
その様子は以下の通り。

```{r}
e <- 0.01
N <- 30
X.homo <- dbinom(0:N,N,e)
beta_A <- beta_B <- N/2
N_A <- N_B <- 0:N
# ヘテロについては理想的なポアソン分布モデルを使用
x_A <- dpois(N_A,beta_A)
x_B <- dpois(N_B,beta_B)
X.hetero <- x_A*x_B[(N+1):1]
X.hetero <- X.hetero/sum(X.hetero)
matplot(0:N,cbind(X.homo,-X.hetero),type="h",main="e=0.01,N=30")
```

また、eの値を小さくすれば、黒の山はごく左端に片寄る。

しかしながら、Nの値が大きい場合に限定すれば、判断できない座位が増えるし、eの値を小さくすれば、活用できるリード数・塩基数が減るので、そのバランスを取ることとなる。
もしくは、バランスを取らずに、各座位のDepth・塩基のクオリティごとに、
False positive,False negativeの確率を算出しながら、情報をすべて活用する、ということになる。

なお、False positive, False negativeの値は、エラー率、ポアソン分布モデルに依存していることも、再度、指摘する。
その想定を変更すれば、おのずと、判断に用いるべき、黒・赤ピークの分離閾値や、False positive, False netative率も変化する。

SNP分割表〜オミックス統計学入門2014

  • こちらの続き
  • SNPのケースコントロール分割表のための配布資料
  • これで90分1回分か2回分と予想
  • 以下はRmdフォーマット。knitrするとhtml化、epub化できます(もちろんフリー)→こちら)
  • その手間を惜しむ人向けにKindleサイトを利用中(こちらにほどなく現れる予定。1US$(処理日の為替レートで円価決定))

# SNPの関連検定

## 3行2列の分割表

一塩基多型(SNP; Single Nucleotide Polymorphism、もしくはSNV; Single Nucleotide Variation)は、DNA配列上のある1塩基が染色体によって異なる箇所のことを言う。

ヒトは常染色体を1対もつので、DNA配列のある位置には1(2本)の染色体の塩基がある。
この個々の染色体が持つ塩基をSNPのアレルと言う。
塩基はA,T,G,Cの4塩基のいずれかである。
二種類の塩基をM,mで表せば、1対の染色体が作るアレルの対はMM,Mm,mmのいずれかとなる。
AとGとをアレルとするSNPではAA,AG,GGの3タイプがあり、この3タイプをMM,Mm,mmと書くことにする、ということである。
アレルの対をディプロタイプと呼んだり、個人が持つ遺伝子多型の型という意味でジェノタイプと呼んだりする。

同じアレルのペアとなっているディプロタイプはホモ接合型、異なるアレルのペアとなっているディプロタイプはヘテロ接合型と言う。

MMとmmはホモ接合型でありMmはヘテロ接合型である。

| |A | B | 行の和 |
|:---:|:---:|:---:|:----:|
|**MM**|  $n_{MM,A}$ | $n_{MM,B}$ |$n_{MM}$ |
|**Mm**| $n_{Mm,A}$ | $n_{Mm,B}$ | $n_{Mm}$ |
|**mm**| $n_{mm,A}$ | $n_{mm,B}$ | $n_{mm}$ |
|**列の和** | $n_{A}$  | $n_{B}$ | $n$ |

# 原因と因果
## ジェノタイプとフェノタイプ
遺伝因子の影響を考えるとき、遺伝子型をジェノタイプ、個体に現れる特徴をフェノタイプと言う。
上の32列の分割表では、ジェノタイプが行、フェノタイプが列となっている。
遺伝疫学研究では、フェノタイプ(A,B)を「病気である(A)」「病気でない(B)」とすることが多い。

ジェノタイプが行、フェノタイプが列である。
22列の分割表のときに、因果関係や条件とその結果を問題にするときには、原因・条件を行とし、結果を列にすることにする、と書いた。
ここでもそのルールを守っている。

ジェノタイプは個体が受精卵という単一の細胞であるときにすでに決まっているタイプであり、すべてのフェノタイプ(個体の特徴)は受精卵が発生して死亡するまでに現れるものであるから、原因と結果の関係がもしあるとすれば、ジェノタイプが原因でフェノタイプが結果である。

従って、ジェノタイプ・フェノタイプの分割表ではジェノタイプが原因・条件を表す行に、フェノタイプが結果を表す列に来る。

『関連因子解析とその臨床応用の基本』にて22列の表を扱ったときには、原因・条件として疾患(X,Y)を行に、検査結果(A,B)を列に取った。
この表では、ジェノタイプを原因・条件として行に、疾患フェノタイプを列に取りたいので、列の疾患フェノタイプを(X,Y)ではなく(A,B)で表していることに注意する。

## ゲノムとそれ以外の"オーム"
DNA配列を全体としてとらえてゲノムと呼ぶ。同様にDNA・染色体が修飾された状態の総体をエピゲノムと言い、転写物の総体をトランスクリプトーム、タンパク質の発現状態全体はプロテオーム、代謝系全体を指してメタボローム、表現型全体をすべて併せてフェノームと言う。

このうちDNA塩基配列以外は、原則として受精卵が発生して個体が死亡するまでの間に変化する。DNA配列も体細胞変異と呼ばれる変化がある。体細胞変異のうち、ある特定のものは癌化と関係している。
このように個体の生涯を通じて変わらないゲノムと、変わるそれ以外とを分けて、前者の変化しないゲノムを狭い意味でのゲノムと言い、その型を(狭い意味での)ジェノタイプと言うことにする。
変化するDNA配列・ゲノムは、癌の場合には「癌ゲノム」と呼ぶが、癌細胞以外にも体細胞変異はあるから、「体細胞変異ゲノム」とでも呼ぶことにする。同様の論理で「体細胞変異ジェノタイプ」と呼ぶことができる。

このように考えると、時間変化しない、「狭義のゲノム」と「時間変化のあるその他」に2分される。同様に「(狭義の)ジェノタイプ」と「フェノタイプ」と分けることができる。

「ジェノタイプ」と「フェノタイプ」とを比較するときには、時間的前後関係がわかるから、ジェノタイプが原因・条件であり、フェノタイプが結果となる。

「フェノタイプ」と「フェノタイプ」とを比較する場合には、時間的前後関係が必ずしも明らかではないから、関連を調べることはできても、因果関係と判断することは難しくなる。

## 3行2列の分割表
32列の分割表に話を戻す。



22列の分割表の色々な計算値を表にすると以下の通りであった。
32列の分割表でも基本的には同じことができるはずであるが、表が異なるので、そのままと言うわけにはいかない。

そのまま適用できるのか、修正・拡張が必要なのかを、個別に検討する。

## 行と列とが独立であるかどうかを評価する場合
### 期待値表

行と列とが独立であると仮定したときの期待値の計算式は

$ne_{i,j} = \frac{n_i n_j}{n}$22列の表の場合と変わらず、表にすれば

|期待値表 |A | B | 行の和 |
|:---:|:---:|:---:|:----:|
|**MM**|  $ne_{MM,A} = \frac{n_{MM}n_A}{n}$ | $ne_{MM,Y}=\frac{n_{MM}n_B}{n}$ |$n_{MM}$ |
|**Mm**|   $ne_{Mm,A} = \frac{n_{Mm}n_A}{n}$ | $ne_{Mm,Y}=\frac{n_{Mm}n_B}{n}$ | $n_{Mm}$ |
|**mm**|   $ne_{mm,A} = \frac{n_{mm}n_A}{n}$ | $ne_{mm,Y}=\frac{n_{mm}n_B}{n}$ | $n_{mm}$ |
|**列の和** | $n_{A}$  | $n_{B}$ | $n$ |

となる。
### $\chi^2$値

この観測表を期待値表に引き比べて、$\chi^2$値を計算することができる。
その式は22列の場合と同じで

$\chi^2 = \sum_{i=MM,Mm,mm;j=A,B}\frac{(n_{i,j}-ne_{i,j})^2}{ne_{i,j}}$

である。

### p値、自由度は2

$\chi^2$値があれば、p値を出すこともできる。
$\chi^2$値からp値を得るには、自由度というものを指定する必要がある。

分割表がa行b列であるときには、自由度が$(a-1)(b-1)$であるので、ここで得られた$\chi^2$値は自由度を$(3-1)(2-1)=2$としてp値に換算することができる。

### 2行2列の式が使えるかどうか

『関連因子解析とその臨床応用の基本』にて22列の表で計算された諸指標は以下の通り。
今32列の表にしたので(X,Y)(MM,Mm,mm)となる。

各指標の式にXとYとが入り乱れて現れている場合には、X,YにMM,Mm,mmのどれを対応付けたらよいのかをすぐには決め難い。
決め難いということは、そのまま使うことができないことを意味する。
逆に、X,Yの片方しか現れない式を持つ指標はそのまま使える。


|そのまま使える|名称 || 意味 | 関係 |
|:---|:---:|:---:|:---:|:----:|
||$\chi^2$|$\sum_{i=MM,Mm,mm;j=A,B}\frac{(n_{i,j}-ne_{i,j})^2}{ne_{i,j}}$ | 帰無仮説からの距離 |p値に換算できる |
||p値| 省略 | 帰無仮説を信じない程度 | $\chi^2$から計算できる |
|×|オッズ比 | $\frac{n_{X,A}n_{Y,B}}{n_{X,B}n_{Y,A}}$  | 行と列との関連の強さ |$\frac{LL(X|A)}{LL(X|B)}$  |
||感度| $\frac{n_{X,A}}{n_{X}}$ | Xという条件でのAの条件付き確率|A観察のときのXの尤度と同じ|
||特異度|$\frac{n_{Y,B}}{n_{Y}}$| Yという条件でのBの条件付き確率|B観察のときのYの尤度|
||A観察のときのXの尤度|$\frac{n_{X,A}}{n_{X}}$|Xという条件でのAの条件付き確率|感度と同じ|
||B観察のときのXの尤度|$\frac{n_{X,B}}{n_{X}}$|Xという条件でのBの条件付き確率||
||A観察のときのYの尤度|$\frac{n_{Y,A}}{n_{Y}}$|Yという条件でのAの条件付き確率||
||B観察のときのYの尤度|$\frac{n_{Y,B}}{n_{Y}}$|Yという条件でのBの条件付き確率|特異度と同じ|
|×|$LL(X/Y|A)$|$\frac{n_{X,A}}{n_{Y,A}}\times \frac{n_{Y}}{n_X}$|A観察のときの尤度比(X/Y)||
|×|$LL(X/Y|B)$|$\frac{n_{X,B}}{n_{Y,B}}\times \frac{n_{Y}}{n_X}$|B観察のときの尤度比(X/Y)||
||Xの事前確率|$\frac{n_X}{n}$|||
||A観察のときのXの事後確率|$\frac{n_{X,A}}{n_A}$||PPVと同じ|
||A観察のときのYの事後確率|$\frac{n_{Y,A}}{n_A}$||1-PPV|
||B観察のときのXの事後確率|$\frac{n_{X,B}}{n_B}$||1-NPVと同じ|
||B観察のときのYの事後確率|$\frac{n_{Y,B}}{n_B}$||NPV|
||PPV|$\frac{n_{X,A}}{n_A}$||A観察のときのXの事後確率と同じ|
||NPV|$\frac{n_{Y,B}}{n_B}$||B観察のときのYの事後確率と同じ|

### 3行2列の表に特徴的な指標

#### 比は2つを比べるもの

##### 尤度比

22列の式をそのまま使えない指標について考えてみる。
それは22列と32列との違いを知る手掛かりになる。

尤度比$LL(X/Y|A) = \frac{n_{X,A}}{n_{Y,A}}\times \frac{n_{Y}}{n_X}$はX,Yの両方を式に含んでいるので、MM,Mm,mmにそのままは使えない。

なぜならA(病気である)という観察がなされたときに、それを引き起こす条件がX,Yの2種類であるならば、その2条件の条件付き確率(尤度)の比をとることはできるが、MM,Mm,mmという3条件を考えているときには、単純に2ジェノタイプを取って比を取ることにわかりやすい意味がないからである。

尤度比が2つの条件の比であるのに、3条件あるという事情がふさわしくないから、と言い換えることもできる。

##### オッズ比

オッズ比$\frac{n_{X,A}n_{Y,B}}{n_{X,B}n_{Y,A}}$も「比」である。
従って、3条件あるときに利用するには、工夫が必要になる。

ジェノタイプが疾患フェノタイプに影響を与えるかどうか、ということに興味があるときには、「疾患を最も起こしにくいジェノタイプ」を基準にして、「あるジェノタイプはどれくらい疾患を起こしやすいか」を知りたい。

今、mmジェノタイプが「疾患を最も起こしにくいジェノタイプ」であるとすれば、MM,Mm,mmの3ジェノタイプのオッズ比は

|オッズ比|||
|:---:|:---:|:---:|
|MM vs. mm | $\frac{n_{MM,A}n_{mm,B}}{n_{MM,B}n_{mm,A}}$  | |
|Mm vs. mm | $\frac{n_{Mm,A}n_{mm,B}}{n_{Mm,B}n_{mm,A}}$  | |
|MM vs. mm | $\frac{n_{mm,A}n_{mm,B}}{n_{mm,B}n_{mm,A}}$  | これは1|

これは(もっともリスクの低いジェノタイプを基準にしたときの)「レラティブ・リスク(相対危険度)」と呼ばれ、分割表のオッズ比がレラティブ・リスクのよい推定値になっていることが多い。

> リスク因子を持っているために、リスク因子が無かった場合に比べて、何倍病気になりやすいか

という値である。したがって、2なら2倍なりやすい、0.5なら1/2倍なりやすい。

## 遺伝的モデルを考慮する〜ジェノタイプに順序がある〜

次の2つの表の$\chi^2$値とオッズ比を計算してみる。
オッズ比の計算に当たっては、mmジェノタイプを基準にすることにする。

| |A | B | 行の和 |
|:---:|:---:|:---:|:----:|
|**MM**|  15 | 10 |25 |
|**Mm**| 20 | 10 | 30 |
|**mm**| 10 | 10 | 20 |
|**列の和** | 45  | 30 | 75 |

| |A | B | 行の和 |
|:---:|:---:|:---:|:----:|
|**MM**|  20 | 10 |30 |
|**Mm**| 15 | 10 | 25 |
|**mm**| 10 | 10 | 20 |
|**列の和** | 45  | 30 | 75 |

この2つの表はMMの行とmmの行の値が入れ替わっているけれど、それ以外はすべて同じである。

$\chi^2$の計算式$\sum_{i=MM,Mm,mm;j=A,B}\frac{(n_{i,j}-ne_{i,j})^2}{ne_{i,j}}$は、3つのジェノタイプについてそれぞれ別々に計算してそれを足し合わせているだけなので、2つの表の$\chi^2$値は等しい。

オッズ比はどうだろうか。

> 上の表では、$OR(MM/mm)=\frac{15\times 10}{10\times 10} =1.5$,$OR(Mm/mm)=\frac{20\times 10}{10\times 10}=2$$OR(mm/mm)=1$> 下の表では、$OR(MM/mm)=\frac{20\times 10}{10\times 10} =2$,$OR(Mm/mm)=\frac{15\times 10}{10\times 10}=1.5$$OR(mm/mm)=1$。

オッズ比だけを抜き出すと、1.5, 2, 1 と、2, 1.5, 1 とである。それぞれは次のように言い換えられる。

> ジェノタイプmmを基準にして、2つのアレルの片方をMに変えたときのレラティブ・リスクが2倍となり、2つのアレルの両方をMに変えるとレラティブ・リスクは少し小さくなり1.5倍となる。


> ジェノタイプmmを基準にして、2つのアレルの片方をMに変えたときのレラティブ・リスクが1.5倍となり、2つのアレルの両方をMに変えるとレラティブ・リスクはさらに大きくなり2倍となる。

どちらの文も真実かもしれないけれど、どちらかと言えば、下の文「2,1.5,1」というオッズ比の並びの方が納得が行くように感じられる。
実際、3ジェノタイプMM,Mm,mmのレラティブ・リスクに順序があることを前提にして解析することの方が普通である。


### ジェノタイプの順序を考慮したモデル

順序があると仮定して解析をすることを、順序があるモデルをあてはめて解析するともいう。
ジェノタイプの順序のモデルにはいくつか代表的なものがあり、モデルにはそれに対応する解析手法がある。


|モデル名|レラティブ・リスクの設定($r>1$)|解析手法名|
|:---:|:---:|:---:|
|優性モデル|RR(MM)=RR(Mm)=$r$,RR(mm)=$1$|2x2表化して$\chi^2$検定|
|劣性モデル|RR(MM)=$r$,RR(Mm)=RR(mm)=$1$|2x2表化して$\chi^2$検定|
|相加モデル|RR(MM)=$2r+1$,RR(Mm)=$r+1$,RR(mm)=$1$|トレンド検定・コクラン・アーミテージ検定|
|相乗モデル|RR(MM)=$r^2$,RR(Mm)=$r$,RR(mm)=$1$|ロジスティック回帰で検定|

このうち、優性モデルと劣性モデルは、いわゆるメンデル型遺伝病に見られるように、責任アレルを1本でも持っていれば発病し(優性モデル)、責任アレルを1本しかもっていない場合には保因者ではあるが疾患フェノタイプは示さず、2本持っている場合には発病する(劣性モデル)というものである。

このようにはっきりと区別できる場合とは異なり、なんらかの影響力を持つアレルがあり、1コピーあるよりも2コピーある方がその影響力が大きいと考えれば、優性モデルと劣性モデルの中間のモデルが必要になる。大多数の遺伝的影響がこのモデルに合致すると思われる。
相加モデルも相乗モデルもどちらもジェノタイプMmに中間的なレラティブ・リスクを設定している。実際、リスクアレルがフェノタイプの発現にどのように影響するかは、個々の分子の生物学的機能の大きさ、発現量、複雑で精密な分子ネットワークによる発現制御、時と場所による発現パターンなどの影響など、複数の要素の複合として現れるので、相加モデルにせよ相乗モデルにせよ、リスクの大きさがモデル式通りになる保証はなく、そうなっていないと考える方が自然である。

### ジェノタイプの順序を考慮するのはなぜか

上述のように、ジェノタイプの順序を考慮したモデル(相加モデル、相乗モデル)はあるが、そのモデルが必ずしも正しいわけではない。
では、なぜそのようなモデルを設定して解析するのだろうか。

端的に言えば、モデルに近いデータを優先して取りあげ、モデルから遠いデータに重きを置かないためである。

例を示す。
次のような3つのデータが得られたとする。

|データ1 |A | B | 行の和 |
|:---:|:---:|:---:|:----:|
|**MM**|  15 | 10 |25 |
|**Mm**| 20 | 10 | 30 |
|**mm**| 10 | 10 | 20 |
|**列の和** | 45  | 30 | 75 |

|データ2 |A | B | 行の和 |
|:---:|:---:|:---:|:----:|
|**MM**|  20 | 10 |30 |
|**Mm**| 15 | 10 | 25 |
|**mm**| 10 | 10 | 20 |
|**列の和** | 45  | 30 | 75 |

|データ3 |A | B | 行の和 |
|:---:|:---:|:---:|:----:|
|**MM**|  16 | 9 |25 |
|**Mm**| 19 | 11 | 30 |
|**mm**| 10 | 10 | 20 |
|**列の和** | 45  | 30 | 75 |

それぞれのオッズ比はmmを基準として、

|データ番号 |OR(MM/mm)|OR(Mm/mm)  |
|:--:|:--:|:--:|
|データ1|1.5|2|
|データ2|2|1.5|
|データ3|1.78|1.73|

この3つの表に、3ジェノタイプの順序を考慮しない検定($\chi^2$検定(自由度2))と、3ジェノタイプの順序を考慮した検定(トレンド検定)を実行してみます

```{r}
d1 <- matrix(c(15,10,20,10,10,10),byrow=TRUE,ncol=2)
d2 <- d1[c(2,1,3),] # 行を入れ替えているだけ
d3 <- matrix(c(16,9,19,11,10,10),byrow=TRUE,ncol=2)
d1
d2
d3
# 順序を考慮しないchi-squared test
chisq.test(d1,correct=FALSE)
chisq.test(d2,correct=FALSE)
chisq.test(d3,correct=FALSE)
# トレンド検定。引数は第1列と、列の和
prop.trend.test(d1[,1],apply(d1,1,sum))
prop.trend.test(d2[,1],apply(d2,1,sum))
prop.trend.test(d3[,1],apply(d3,1,sum))
```

ジェノタイプの順序を考慮しない検定では、データ1とデータ2との$\chi^2$が等しく、データ3のそれはより小さくなっている(X-squared = 1.389,1.389,1.139)。
これは、ジェノタイプ順序を考慮しないのであれば、データ1とデータ2は同程度にジェノタイプがフェノタイプに影響していると考えるのがよく、データ3はそれらに比べて影響があると考えにくい、と読み取れる。
一方、順序を考慮したトレンド検定の方では、統計量が0.3731,1.37,0.8396とあるように、もっとも関連がありそうなのはデータ2、次がデータ3、最後にデータ1が来ることがわかる。

このようにジェノタイプに順序を入れたモデルを適用することで、モデルに合致しないデータには小さい関連指標値が出るようになる。
このような指標値を用いることで、モデルに合致したデータが優先的にスクリーニングすることができる。

## 活用する
以上が、ディプロタイプ型ジェノタイプと2値型フェノタイプに観察される32列分割表の基本事項である。
以下は、その知識を活用するための練習である。

### データを作ってみる
#### Hardy-Weinberg 平衡とランダム・メイティング

今、ある集団にアレルMとアレルmとが割合$P(M)=0.7$,$P(m)=0.3$で存在しているとする。
これはアレルMの染色体が全体の7割、アレルmの染色体が3割である、という意味である。
この集団がランダム・メイティングしており、このアレルが出生・成長・生殖に影響がないとすると、この集団の個体ジェノタイプMM,Mm,mmの頻度は$P(MM)=P(M)^2$,$P(Mm)=2P(M)P(m)$,$P(mm)=P(m)^2$になると考えられる。
アレル頻度とジェノタイプ頻度との関係がこのような関係にあることをHardy-Weinberg平衡(HWE)にある、と言う。

リスクアレルmの頻度が$P(m)=0.2$であり、ジェノタイプ頻度がHWEになっているような集団から、$N_s=100000$人をランダムにサンプリングし、そのジェノタイプを長さ$N_s$のベクトルにし、それを集計せよ。ただし、ジェノタイプMM,Mm,mmにはリスクアレルmの本数0,1,2を対応づけること。

```{r}
P.m <- 0.2
N.s <- 100000
P.M <- 1-P.m
P.MM <- P.M^2; P.Mm <- 2*P.M*P.m; P.mm <- P.m^2
P.g <- c(P.MM,P.Mm,P.mm)
G <- sample(0:2,N.s,replace=TRUE,prob=P.g)
table(G)
```

#### ジェノタイプ別に有病率を設定する

リスクアレルを持たない場合、この病気の有病率を$Q(MM)=0.1$とし、レラティブ・リスクには相乗モデルを仮定し、ジェノタイプMmのレラティブ・リスクを$r=1.5$とし、$N_s$人に確率的にフェノタイプXを割り当てよ。
```{r}
Q.MM <- 0.1
r <- 1.5
Q.Mm <- Q.MM * r
Q.mm <- Q.MM * r^2
Q.G <- c(Q.MM,Q.Mm,Q.mm)
Q.G
Q.s <- Q.G[G+1] # これで各個人のRRが得られる
plot(G,Q.s) # 念のため確認する
plot(jitter(G),jitter(Q.s)) # 点が重なって不安なので少しぶれを入れると安心するかもしれない
R <- runif(N.s) # 各個人用の乱数
X <- R < Q.s # 病気がTRUE,非病気がFALSE
X <- as.numeric(X) # TRUEを1、FALSEを0に換える
table(data.frame(G=G,X=X))
```

#### データを作る処理を関数にする

人数、アレル頻度、レラティブ・リスクを与えて、個人別のフェノタイプを発生させ、さらにそれを32列の表にした。この処理を関数にする。
```{r}
# 入力変数
N.s <- 100000
P.m <- 0.2
Q.MM <- 0.1
r <- 1.5

# 出力
# ジェノタイプG,フェノタイプX,そのテーブルTa
# まずはその枠だけ作ります
my.make.SNPdata <- function(N.s,P.m,Q.MM,r){
  
  return(list(G=G,X=X,Ta=Ta)) # GをGとして、XをXとして、TaをTaとして返す
}
```

関数の入力と出力ができたら、あとは入力の変数と出力の変数をつなぐ処理を関数の{ }の内側に書き込む。
```{r}
# 入力変数
N.s <- 100000
P.m <- 0.2
Q.MM <- 0.1
r <- 1.5

# 出力
# ジェノタイプG,フェノタイプP,そのテーブルTa
# まずはその枠だけ作る
my.make.SNPdata <- function(N.s,P.m,Q.MM,r){
  P.M <- 1-P.m
  P.MM <- P.M^2; P.Mm <- 2*P.M*P.m; P.mm <- P.m^2
  P.g <- c(P.MM,P.Mm,P.mm)
  G <- sample(0:2,N.s,replace=TRUE,prob=P.g)
  Q.Mm <- Q.MM * r
  Q.mm <- Q.MM * r^2
  Q.G <- c(Q.MM,Q.Mm,Q.mm)
  Q.s <- Q.G[G+1] # これで各個人のRRが得られる
  R <- runif(N.s) # 各個人用の乱数
  X <- R < Q.s # 病気がTRUE,非病気がFALSE
  X <- as.numeric(X) # TRUEを1、FALSEを0に換える
  Ta <- table(data.frame(G=G,X=X))
  return(list(G=G,X=X,Ta=Ta))
}
# 使ってみます
test.out <- my.make.SNPdata(N.s,P.m,Q.MM,r)
test.out$Ta
```

### 関連を見出す

データをシミュレーション作成できるようになったので、これを用いて解析することにする。
大きく分けると次の2点をになる。
> 関連があることを、「関連がないという帰無仮説を棄却」するに足るp値が得られたことで示し、

> 関連の強さをジェノタイプ別のレラテイィブ・リスクの推定値として示す。その推定値は分割表から計算するオッズ比であったり、ロジスティック回帰の係数から算出したものであったりする。推定値は点推定値と信頼区間とで示す。

では、上で作成した集団においてケース・コントロール関連解析を実施したものとして、解析をしてみる。
ケース・コントロール関連解析では、フェノタイプごとにサンプリングする。
今、ケースとコントロールをそれぞれ$N_{ca}=500,N_{co}=500$人サンプリングして相加モデルで検定してみる。

```{r}
N.ca <- 500
N.co <- 500
# フェノタイプtest.out$Xの値が1,0の番号を取りだす
cases <- which(test.out$X==1)
conts <- which(test.out$X==0)
# N.ca,N.co人をサンプリングする
cases.sample <- sample(cases,N.ca)
conts.sample <- sample(conts,N.co)
samples <- c(cases.sample,conts.sample)
# サンプルのジェノタイプとフェノタイプ
G.sample <- test.out$G[samples]
X.sample <- test.out$X[samples]
# サンプルの分割表
Ta.sample <- table(data.frame(G.sample,X.sample))
Ta.sample
tr.out <- prop.trend.test(Ta.sample[,1],apply(Ta.sample,1,sum))
tr.out
```

大きな$\chi^2$値と小さなp値が得られる。

> ジェノタイプとフェノタイプには関連がないと考えにくい。
> その判断は相加モデルによって検定した。
> その意味は、3ジェノタイプのリスクが$2r+1,r+1,1$となる場合を最優先としてジェノタイプとフェノタイプの関連を評価・検定した、という意味である。

次にMMを基準として、オッズ比を計算する。
オッズ比は2つのジェノタイプが作る22列の表を用いて計算する。
そのためにパッケージをインストールしてそれを読み込んだ上で、オッズ比とその95%信頼区間を計算する関数を使う。
その関数は、1行目が因子あり、1列目が病気ありを前提にしているので、関数への22列データの受け渡しに際して、それに合うようにしている。
```{r}
# "epiR"パッケージのインストールが済んでいなければ、以下のコメントアウトした行を実行します
#install.packages("epiR")
library(epiR)
OR.Mm <- epi.2by2(Ta.sample[c(2:1),c(2:1)])
OR.mm <- epi.2by2(Ta.sample[c(3,1),c(2:1)])
OR.Mm
OR.mm
```

複数の指標が示されるが、ここでは、Odds ratioに着目する。
点推移定値と信頼区間とが x (y,z) という形で示されている。

データをシミュレーションした際に設定したレラティブ・リスクは、Mm,mmそれぞれ、$1.5,1.5^2=2.25$であった。
サンプルから計算したオッズ比はその値から遠くはないが、等しいわけではない。

また、検定モデルとしては相加モデルを用いたが、それはMmのリスクがmmのそれよりも小さいことを仮定したかったからであって、MM,Mm,mmのリスクが$1,1+r,1+2r$でなければならないと仮定したわけではないので、Mm,mmのレラティブ・リスクは、上記のような計算をしてx (y,z)と信頼区間付きで考えておけばよい。

さて、場合によってはロジスティック回帰を手法として採用するかもしれない。
以下のように実施することができる
```{r}
fit <- glm(X.sample~G.sample,family=binomial())
summary(fit) # display results
```

ここに現れるのは、相乗モデルが最も適切だと考えた上で、ジェノタイプとフェノタイプとに関連があるかどうかの評価結果である。
いくつもの指標と数値が現れるが、prop.trend.test()で算出した$\chi^2$値とp値とに相当するのはCoefficients:に引き続くG.sample のz value とPr(>|z|)の値である。

G.sampleとして与えた、ジェノタイプの値がフェノタイプに与える影響をzという値で表し、それに対応するp値をPr(|z|)と書いている。
$z^2$の値が$\chi^2$値と対応するので、prop.trend.test()$\chi^2$値とこのz値とを、両手法のp値とともに比べてみる。
```{r}
print(c(tr.out$statistic,(summary(fit))$coefficients[2,3]^2))
print(c(tr.out$p.value,(summary(fit))$coefficients[2,4]))
```
だいたい同じであることがわかる。
関連がないという帰無仮説を棄却する上でだいたい同じである、と言うことである。

ジェノタイプMMを基準にしたオッズ比とその信頼区間はモデルによらず算出することができるので、ロジスティック・回帰モデルで解析した後にも、その値を計算し、それを利用することはもちろん可能である。
ロジスティック回帰を用いた場合には、リスク・アレル1本の寄与の程度が計算されているので、それに基づいて、Mm,mmジェノタイプのレラティブ・リスクを推定すれば以下のようになる。

Mmは回帰係数$\beta$を用いて$e^{\beta}$で計算されるから、以下のようになる。
初めに点推定値、次に信頼区間を示す。
```{r}
print(c(exp(coef(fit))[2],exp(confint(fit))[2,])) # 95% CI 信頼区間
```

mmは2本分であるから、1本分の2乗であり、次のようになる。

```{r}
print(c(exp(coef(fit))[2]^2,exp(confint(fit))[2,]^2)) # 95% CI 信頼区間
```

この値も、シミュレーションで用いた1.5,2.25と完全には一致しないが、推定値としては妥当であることが確かめられる。

#### 結果を応用する

みずからサンプルを収集し関連解析を行う場合には、上記のような手順を踏む。
では、そのような研究成果が出ているとき、それはどのように臨床応用すればよいかを考える。

32列の分割表とprop.trend.test()のp値とMM基準のジェノタイプ別オッズ比が報告されたとする。
```{r}
Ta.sample
tr.out$p.value
OR.Mm$rval$OR[c(1,4,5)]
OR.mm$rval$OR[c(1,4,5)]
```

ある人のジェノタイプがmmだったとき、この報告に基づいて

> あなたのジェノタイプはmmなので、あなたが病気になる確率は、MMの人に比べて

```{r}
paste(OR.mm$rval$OR[1],"倍")
```

高いです、と説明することは可能である。
しかしながら、その説明を受けても、結局、その人は自分がどれくらいの確率で病気になるのかがわからない。

これは、ケース・コントロールサンプリングをしたときに、ケースの人数とコントロールの人数の比が、
集団におけるケースとコントロールの人数の比と異なっているために生じる問題である。

今、ある人が病気かそうでないかの情報が全くないものとする。
この人のジェノタイプが判明したとき、この人が病気なのかそうでないのかを予想することを考える。

これは、疾患と臨床検査の分割表
| |A(検査が陽性である) | B(検査が陰性である) | 行の和 |
|:---:|:---:|:---:|:----:|
|**X(病気であるならば)**|  $n_{X,A}$ | $n_{X,B}$ |$n_{X}$ |
|**Y(病気でないならば)**| $n_{Y,A}$ | $n_{Y,B}$ | $n_{Y}$ |
|**列の和** | $n_{A}$  | $n_{B}$ | $n$ |
で言うところの、A検査が陽性であるをある特定のジェノタイプである、と読み換えたときの
PPVと同じである。

この計算をするためには、X,Yの人数比がわからなければ、分割表自体がPPVの計算にそぐわない。

今、ケース・コントロールスタディを実施したので、$n_{ca}=n_{X},n_{co}=n_Y$が一般集団のそれとは異なっている。
一般集団での$n_{ca}=n_X,n_{co}=n_Y$の割合は有病率Qによって$Q,1-Q$となるから、有病率の情報が必要であることがわかる。

今、$Q=0.12$だと考えられていて、ある人のジェノタイプがmmであると判明したとする。
Qを考慮して、mmかそうでないかで22列の表の値を再計算することにする。
ここではジェノタイプを列に取ることにする。ジェノタイプとの関係では、病気か否かはA,Bで表して来たことに注意する。

| |mm | mm以外 | 行の和 |
|:---:|:---:|:---:|:----:|
|**A(病気である)**| $n_1=$? | $n_2=$? |$n\times Q$ |
|**B(病気でない)**| $n_3=$? | $n_4=$? | $n \times (1-Q)$ |
|**列の和** | $n_1+n_3=$?  | $n_2+n_4=$? | $n$ |

で?の値をケース・コントロールの観察表から算出したい。

$n_1 : n_2 = n_{mm,A} : n_{MM,A}+n_{Mm,A}$ であり、$n_3 : n_4 = n_{mm,B} : n_{MM,B}+n_{Mm,B}$であることを考慮して、

| |mm | mm以外 | 行の和 |
|:---:|:---:|:---:|:----:|
|**A(病気である)**| $nQ\frac{n_{mm,A}}{n_A}$ | $nQ\frac{n_{MM,A}+n_{Mm,A}}{n_{A}}$ |$n\times Q$ |
|**B(病気でない)**| $n(1-Q)\frac{n_{mm,B}}{n_{B}}$ | $n(1-Q)\frac{n_{MM,B}+n_{Mm,B}}{n_{B}}$ | $n \times (1-Q)$ |
|**列の和** | $n_1+n_3$ | $n_2+n_4$| $n$ |

という表が得られます。

この表に基づいて、いわゆる$PPV = \frac{n_1}{n_1+n_3}$を計算すれば、ジェノタイプがmmとわかったときの病気である確率が得られます。
式で書くと

$\frac{Q\frac{n_{mm,A}}{n_A}}{Q\frac{n_{mm,A}}{n_A}+(1-Q)\frac{n_{mm,B}}{n_{B}}}$

実際にQを与えてこの値を計算してみると
```{r}
Q <- 0.125
n.mm.A <- Ta.sample[3,2]
n.A <- sum(Ta.sample[,2])
n.mm.B <- Ta.sample[3,1]
n.B <- sum(Ta.sample[,1])

(Q*n.mm.A/n.A)/(Q*n.mm.A/n.A + (1-Q)*n.mm.B/n.B)
```

となる。
ジェノタイプ情報がなかったときに$Q=0.125$と予想していたものが、かなり大きな値となった。
これが、ジェノタイプが判明したときの病気である事後確率である。

### ジェノタイプ情報を用いた疾病確率の計算に関する注意

上記の計算で明らかにしたのは、実は、「病気かどうかわからない」状態のときに、「ジェノタイプがmmと判明したとして」、その人が「有病率Q」で表される「疾病集団の一員である」(事後)確率である。

今、Qが生まれてから死ぬまでにこの病気にかかる人の割合(生涯有病率)を表し、ジェノタイプを調べ、病気になる確率を計算した対象者が新生時であれば、この計算でよい。

しかしながら、そもそも生涯有病率にしろ、ある一時点での有病率(時点有病率)も正確な数字は不明なので、上の計算は、幅のあるおよその数字である。

また、実際には「ある年齢の人」が「これから病気になる」確率を知りたいとすると、「その年齢までは病気になっていない人のうち、その年齢以降に病気になる割合」をQとしなくてはならない。
この割合は、年齢とともに下がるが、どのように下がるかがわかっている疾患はあまりない。
たとえば、20歳の人が死ぬまでに肺癌になる確率を訪ねる場合と、85歳の人が同じ質問をする場合とでは異なる考え方が必要だと言うことである。

このように、有病率Qに関する曖昧さや難しさの他に、レラティブ・リスクの信頼区間のことも考慮に入れる必要がある。
上の計算では、観測表の度数をそのまま使っているので、オッズ比で言うところの点推定値のみを使い、信頼区間は無視していることとに相当する。
有病率の曖昧さや年齢の要素の曖昧さがあるので、オッズ比に関するところだけを厳密にしても、適当とは考えない。
したがって、今の段階では、

> 「複数の理由で曖昧さのある点推定値」を「ジェノタイプ情報を得た後の発病事後確率」の目安として計算できる

ということに留める。
時間があれば、この点について回を改めて取り扱う。

オミックス統計学入門2014

  • 京都大学大学院医学研究科の修士対象講義は、2012年度に「インフォマティクス」2013年度に「統計学処理の基礎概念」を扱っています。2014年度は、遺伝疫学やオミックスデータ解析の論文の主張を読み取ったり、それを臨床展開することを念頭に置いた「解析初心者」を相手にした一番「簡単な」内容になる予定です(こちら)。
  • そのテキストの1つ目
  • Rmdフォーマットで書いてepub化してみました

(こちら)

  • 以下をRmdファイルとして保存してknitr()するだけ(やり方解説)
  • これで90分1回分かなー
  • こちらにそのうち現れる予定

# 関連因子解析とその臨床応用の基本
ここでは、因子が因子と関連があることと、ある因子についての情報を知りえたときに、その因子と関連のある別の因子について予測するということに関する、基礎事項を扱います。

「関連がある」ことが論文で報告されているときに、そのような関連はどのように示されるのかという点は前者に、
「関連すると報告された因子」を臨床に活用する場合にどう考えるのかという点は後者に相当します。

この基本を22列の分割表を例にして扱います。

## 分割表の基本:2行2列の分割表
### 行と列
| |A | B | 行の和 |
|:---:|:---:|:---:|:----:|
|**X**|  $n_{X,A}$ | $n_{X,B}$ |$n_{X}$ |
|**Y**| $n_{Y,A}$ | $n_{Y,B}$ | $n_{Y}$ |
|**列の和** | $n_{A}$  | $n_{B}$ | $n$ |

分割表の基本は2x2分割表です。行X、行Yと2行あり、A列、B列と列数も2です。
行の和、列の和は表の周辺にある数なので、周辺度数と言います。

## 行と列に関する2つの見方
### 関連
#### 関連と掛け算、期待値

行と列とに関連がないときには、行の和と列の和との情報(周辺度数)だけがあれば、表のセルの値を決めることができると考え、その値を「関連がない」という仮定の下での期待値と言い、そのような値で出来た表を「関連がない」という仮定の下での期待値表と言います。

|期待値表 |A | B | 行の和 |
|:---:|:---:|:---:|:----:|
|**X**|  $ne_{X,A}=\frac{n_X \times n_A}{n}$ | $ne_{X,B}=\frac{n_X \times n_A}{n}$ |$n_{X}$ |
|**Y**| $ne_{Y,A}=\frac{n_X \times n_A}{n}$ | $ne_{Y,B}=\frac{n_X \times n_A}{n}$ | $n_{Y}$ |
|**列の和** | $n_{A}$  | $n_{B}$ | $n$ |

#### 関連の強さ
##### $\chi^2$とp値
関連の強さは「関連がない」という仮定の下での期待値から、観測表がずれている程度で測ります。

$\chi^2 = \frac{(n_{X,A}-ne_{X,A})^2}{ne_{X,A}}+\frac{(n_{X,B}-ne_{X,B})^2}{ne_{X,B}}+\frac{(n_{Y,A}-ne_{Y,A})^2}{ne_{Y,A}}+\frac{(n_{Y,B}-ne_{Y,B})^2}{ne_{Y,B}}$

$\chi^2 = \sum_{i=X,Y,j=A,B}\frac{(n_{i,j}-ne_{i,j})^2}{ne_{i,j}}$

と書くこともできます。
分割表のセルの値がちょうど期待値になっているときに$\chi^2=0$であり、それ以外は必ず$\chi^2>0$となっています。

ずれの程度を計算する方法として、このようにすると好都合であり、この値を、2x2表の行と列とが関係ないかどうかを測るための$\chi^2$(カイ二乗値)と言います。

$\chi^2$値を用いて、「帰無仮説が正しいとすると、観測表が得られるのは珍しいことになるから、帰無仮説は信じないことにしよう〜帰無仮説の棄却」をするときには、$\chi^2$値をp値に変換します。
p値は小さいほど「帰無仮説を信じない方がよい」ことになります。

p値の小ささは「帰無仮説を信じない方がよい」程度の指標になりますが、関連の強さを評価するには、次のオッズ比も必要です。

##### オッズ比

$OR=\frac{n_{X,A}n_{Y,B}}{n_{X,B}n_{Y,A}}$という値をオッズ比と言います。

この値は、分割表のセルの値がちょうど期待値になっているときに$OR=1$になり、それ以外では、$1$より大きくなったり、小さくなったりします。$OR,\frac{1}{OR}$の値が$1$から離れれば離れるほど、帰無仮説から遠いことになります。

関連のありなし、関連の強さを知るためには、$\chi^2$とそれに対応するp値とを使って、帰無仮説が正しくないと信じるための根拠を与え、関連の強さを表す数値としてオッズ比を求めます。

行と列とを入れ替えても値が変わらないことにも注意します。

#### 関連の相対性、行と列の入れ替え
(A,B)の持ち具合と(X,Y)の持ち具合に関連があるかどうかを問題にします。
(A,B)(X,Y)とは対等なので、行と列とを入れ替えても構いません。

| |X | Y | 行の和 |
|:---:|:---:|:---:|:----:|
|**A**|  $n_{X,A}$ | $n_{X,B}$ |$n_{A}$ |
|**B**| $n_{Y,A}$ | $n_{Y,B}$ | $n_{B}$ |
|**列の和** | $n_{X}$  | $n_{Y}$ | $n$ |

関連の強さを表す$\chi^2$の値は、行列を入れ替えても変わりません。

### 因果・条件、感度・特異度、尤度・尤度比、事前確率・事後確率、PPV・NPV
#### 行と列との入れ替えは、不可
関連を考えるときには、行と列とを入れ替えることが可能でしたが、今度は入れ替えてはいけない場合を扱います。
XであるならばA(または)B、YであるためにA(またはB)である、というような場合には、(X,Y)は原因や条件を表し、(A,B)はその結果を表しています。
このような場合は、(X,Y)(A,B)とを入れ替えることはできません。

以降では、原因・条件とその結果とを問題にするときには、原因を行((X,Y)の方)に、結果を列((A,B)の方)に取ることにします。

#### 感度と特異度、条件付き確率
次のような疾患診断のための検査を考える。

| |A(検査が陽性である) | B(検査が陰性である) | 行の和 |
|:---:|:---:|:---:|:----:|
|**X(病気であるならば)**|  $n_{X,A}$ | $n_{X,B}$ |$n_{X}$ |
|**Y(病気でないならば)**| $n_{Y,A}$ | $n_{Y,B}$ | $n_{Y}$ |
|**列の和** | $n_{A}$  | $n_{B}$ | $n$ |

病気であるという条件の下(X)での検査の陽性(A)(条件付き確率)$\frac{n_{X,A}}{n_{X}}$と計算され、それは感度と言われる。
同様に病気でないという条件の下(Y)での検査の陰性(B)(条件付き確率)$\frac{n_{Y,B}}{n_{Y}}$と計算され、それは特異度と言われる。

#### PPV(Positive Predictive Value;陽性的中率)とNPV(Negative Predictive Value;陰性的中率)、尤度・尤度比、事前確率・事後確率
次のような疾患診断のための検査を考える。

| |A(検査が陽性である) | B(検査が陰性である) | 行の和 |
|:---:|:---:|:---:|:----:|
|**X(病気であるならば)**|  $n_{X,A}$ | $n_{X,B}$ |$n_{X}$ |
|**Y(病気でないならば)**| $n_{Y,A}$ | $n_{Y,B}$ | $n_{Y}$ |
|**列の和** | $n_{A}$  | $n_{B}$ | $n$ |

> 「病気であると、検査が陽性になる」

と考えているので、「病気である・ない」を行に、「検査の陽性・陰性」を列にとっている。

##### PPV
検査が陽性(A)であったという事実が与えられた下で、病気である(X)か、病気でない(Y)かは、それぞれ

> 病気である条件の下での検査の陽性率(条件付き確率)

> 病気ではない条件の下での検査の陽性率(条件付き確率)

であり、$\frac{n_{X,A}}{n_{X}}$,$\frac{n_{Y,A}}{n_{Y}}$と計算される。
このとき、この2つの値は、「検査が陽性であるという事実があるとき」の病気である・ないという仮説の尤度と呼ばれる。

2つの尤度の比は検査が陽性のときの尤度比と呼ばれ、$LL(X/Y|A)=\frac{\frac{n_{X,A}}{n_{X}}}{\frac{n_{Y,A}}{n_{Y}}}=\frac{n_{X,A}}{n_{Y,A}}\times \frac{n_{Y}}{n_X}$である。

この検査が陽性のときの尤度比は、検査結果から病気であるかどうかを判断するときに役に立つ値である。

検査結果が無かったときに病気であるかないかを$\frac{Pre(X)}{Pre(Y)}$の比で想定していたとしたら、$\frac{Pre(X)}{Pre(Y)}\times LL(X/Y|A)$に考えを改める、という意味がある。

実際、検査をする前に、$\frac{n_{X}}{n_{Y}}$の比であると思っていたとすると

> $\frac{n_X}{n_Y}\times LL(X/Y|A) = \frac{n_X}{n_Y}\times \frac{n_{X,A}}{n_{Y,A}}\times \frac{n_{Y}}{n_X} = \frac{n_{X,A}}{n_{Y,A}}$

となる。

検査をする前に思っていた比$\frac{z}{w}$$z,w$を使って$\frac{z}{z+w}$とすると、これは、検査前に病気であると考える「事前確率」。

また、検査後の比$\frac{n_{X,A}}{n_{Y,A}}$を使った、$\frac{n_{X,A}}{n_{X,A}+n_{Y,A}} = \frac{n_{X,A}}{n_A}$は「事後確率」。

この検査が優性のときの尤度比と同様に、$\frac{n_{X,A}}{n_A}$も病気であるかないかの判断に役立つ。これをPPVと言う。

##### NPV

同様に、検査が陰性(B)であったという事実が与えられた下で、病気である(X)か、病気でない(Y)かは、それぞれ
> 病気である条件の下での検査の陰性率(条件付き確率)

> 病気ではない条件の下での検査の陰性率(条件付き確率)

であり、$\frac{n_{X,B}}{n_{X}}$,$\frac{n_{Y,B}}{n_{Y}}$と計算される。
このとき、この2つの値は、「検査が陽性であるという事実があるとき」の病気である・ないという仮説の尤度と呼ばれる。

2つの尤度の比は検査が陰性のときの尤度比と呼ばれ、$LL(X/Y|B)=\frac{\frac{n_{X,B}}{n_{X}}}{\frac{n_{Y,B}}{n_{Y}}}=\frac{n_{X,B}}{n_{Y,B}}\times \frac{n_{Y}}{n_X}$である。

この検査が陰性のときの尤度比は、検査結果から病気であるかないかを判断するときに役に立つ値である。

検査結果が無かったときに病気であるかないかを$\frac{Pre(X)}{Pre(Y)}$の比で想定していたとしたら、$\frac{Pre(X)}{Pre(Y)}\times LL(X/Y|B)$に考えを改める、という意味がある。

実際、検査をする前に、$\frac{n_{X}}{n_{Y}}$の比であると思っていたとすると

> $\frac{n_X}{n_Y}\times LL(X/Y|B) = \frac{n_X}{n_Y}\times \frac{n_{X,B}}{n_{Y,B}}\times \frac{n_{Y}}{n_X} = \frac{n_{X,B}}{n_{Y,B}}$

となる。

検査をする前に思っていた比$\frac{z}{w}$$z,w$を使って$\frac{z}{z+w}$とすると、これは、検査前に病気であると考える「事前確率」であるのは、PPVのときと同じ。

また、検査後の比$\frac{n_{X,B}}{n_{Y,B}}$を使った、$\frac{n_{X,B}}{n_{X,B}+n_{Y,B}} = \frac{n_{X,B}}{n_B}$は「事後確率」。

この検査が優性のときの尤度比と同様に、$\frac{n_{X,B}}{n_B}$$1-\frac{n_{X,B}}{n_B}=\frac{n_{Y,B}}{n_B}$病気であるかないかの判断に役立つ。$\frac{n_{Y,B}}{n_B}$をNPVと言う。

## 2行2列分割表の色々な計算値
|名称 || 意味 | 関係 |
|:---:|:---:|:---:|:----:|
|$\chi^2$|$\sum_{i=X,Y,j=A,B}\frac{(n_{i,j}-ne_{i,j})^2}{ne_{i,j}}$ | 帰無仮説からの距離 |p値に換算できる |
|p値| 省略 | 帰無仮説を信じない程度 | $\chi^2$から計算できる |
|オッズ比 | $\frac{n_{X,A}n_{Y,B}}{n_{X,B}n_{Y,A}}$  | 行と列との関連の強さ |$\frac{LL(X|A)}{LL(X|B)}$  |
|感度| $\frac{n_{X,A}}{n_{X}}$ | Xという条件でのAの条件付き確率|A観察のときのXの尤度と同じ|
|特異度|$\frac{n_{Y,B}}{n_{Y}}$| Yという条件でのBの条件付き確率|B観察のときのYの尤度|
|A観察のときのXの尤度|$\frac{n_{X,A}}{n_{X}}$|Xという条件でのAの条件付き確率|感度と同じ|
|B観察のときのXの尤度|$\frac{n_{X,B}}{n_{X}}$|Xという条件でのBの条件付き確率||
|A観察のときのYの尤度|$\frac{n_{Y,A}}{n_{Y}}$|Yという条件でのAの条件付き確率||
|B観察のときのYの尤度|$\frac{n_{Y,B}}{n_{Y}}$|Yという条件でのBの条件付き確率|特異度と同じ|
|$LL(X/Y|A)$|$\frac{n_{X,A}}{n_{Y,A}}\times \frac{n_{Y}}{n_X}$|A観察のときの尤度比(X/Y)||
|$LL(X/Y|B)$|$\frac{n_{X,B}}{n_{Y,B}}\times \frac{n_{Y}}{n_X}$|B観察のときの尤度比(X/Y)||
|Xの事前確率|$\frac{n_X}{n}$|||
|A観察のときのXの事後確率|$\frac{n_{X,A}}{n_A}$||PPVと同じ|
|A観察のときのYの事後確率|$\frac{n_{Y,A}}{n_A}$||1-PPV|
|B観察のときのXの事後確率|$\frac{n_{X,B}}{n_B}$||1-NPVと同じ|
|B観察のときのYの事後確率|$\frac{n_{Y,B}}{n_B}$||NPV|
|PPV|$\frac{n_{X,A}}{n_A}$||A観察のときのXの事後確率と同じ|
|NPV|$\frac{n_{Y,B}}{n_B}$||B観察のときのYの事後確率と同じ|

## Rの練習
22列の行列を分割表に見立て、それぞれの値の計算をしてみること。
```{r}
# 行列を作る
t <- matrix(sample(10:20,4),2,2)
t
# 周辺度数を出す
n.row <- apply(t,1,sum)
n.col <- apply(t,2,sum)
n <- sum(t)
n.row
n.col
n
# 期待値表を作る
t.e <- n.row %*% t(n.col) / n
t.e
# 期待値表の周辺度数を確認する
apply(t.e,1,sum)
apply(t.e,2,sum)
# chi^2
chi2 <- sum((t-t.e)^2/t.e)
chi2
# p値化する
pchisq(chi2,df=1,lower.tail=TRUE)
# chisq.test()を使う。ただし、Yates補正はしない
chisq.test(t,correct=FALSE)
# OR
t[1,1]*t[2,2]/(t[1,2]*t[2,1])
# 感度
t[1,1]/n.row[1]
# 特異度
t[2,2]/n.row[2]
# LL(X/Y|A)
t[1,1]/t[2,1] * (n.row[2]/n.row[1])
# Xの事前確率
n.row[1]/n
# PPV
t[1,1]/n.col[1]
# NPV
t[2,2]/n.col[2]
```