この文書は京都大学医学部法医学講座にて2010年から始まった法数学勉強会で提供した話題をまとめたものです。
同勉強会はDNA多型を用いた個人識別を取り巻く話題の数学的側面を勉強しようと言う会で、
私はDNA多型に関する統計学を遺伝疫学分野にて専門にしていることから、法医学講座の玉木敬二教授にお誘いいただいて参加してきたものです。
多型についてはわかっていても法医学分野での問題設定・情報の解釈と活用は、私にとってまったく初めてのテーマであり、参加のたびに多くのことを勉強させていただいております。
このように門外漢として出発した私の文書ですので、内容は的外れであったり、基本的な誤解があったりするかもしれません。
その点はあらかじめご容赦いただきたいと思います。
それとは逆に門外漢であるが故に、法医学分野の常識にとらわれない考えもあるかもしれず、それについては、ポジティブな側面があるかもしれないという温かい眼で見ていただき、多少なりとも参考になることがあれば幸いです。
本文書は章立てをしてはいますが、全体の統一をとるための作業をほとんどしておりません。
複数の章で同様の記述が繰り返されることがあるかもしれませんが、ご容赦ください。
勉強会で使用したスライドは[こちら](http://www.genome.med.kyoto-u.ac.jp/StatGenet/lectures/2010/HouSuugaku/DNAevidence2010Aug.ppt 勉強会で使用したスライド)です。
以下の文章と勉強会での話の内容には違う点もありますが、参考までにリンクを張っておきます。
DNA鑑定で一番多い状況は、「このDNA試料の型が仮説Xと仮説Yとのどちらによりよく合致するか」を尤度比で数値化するものだと思われます。
この章では、特定２仮説の比較ではなく、
比較するべき仮説が２つとは限らない場合とはどういう場合なのかを問題にします。
ある仮説の下で、ある事象が起きる確率を
\(
Pr(\mathbf{x}|\mathbf{\theta})
\)
と書きます。
\(\mathbf{x}\)はいろいろな事象を、\(\mathbf{\theta}\)は仮説の条件を表します。
\(\mathbf{x}\)と\(\mathbf{\theta}\)との間の縦棒"|"は\(A|B\)で
>『Bが成り立つという条件でのA』
と読みます。
従ってこの式は
>『\(\mathbf{\theta}\)が成り立っているときに\(\mathbf{x}\)が起きる確率』
と読みます。
条件が付いた確率なので、条件付き確率と呼びます。
\(\mathbf{x}\)は太字で書いてあります。これは\(\mathbf{x}\)は一つだけでなくたくさん(場合によっては\(0 \le x \le 1\)のようにある範囲について数限りない値のことを想定する場合もあります。
DNA鑑定では、個人が持つジェノタイプ(複数マーカーの組合せジェノタイプ)の一つ一つが
\(\mathbf{x}\)の一つ一つの要素に対応します。
だれしも、何かしらのジェノタイプを持ちますし、あるジェノタイプを持ったら、他のジェノタイプは持てませんから、\(\mathbf{x}\)のすべてについて確率を計算して
その値を足しあわせれば、１になります。
\(
\sum_{x \in \mathbf{x}} Pr(\mathbf{x}|\mathbf{\theta}) =1
\)
この式は、
>\(\mathbf{x}\)が離散的な場合に全事象の条件付き確率を足し合わせると\(1\)になる
と言っています。
\(
\int_{x \in \mathbf{x}} Pr(\mathbf{x}|\mathbf{\theta}) dx =1
\)
この式は\(\mathbf{x}\)が連続的な場合に同じことを言っています。
DNA鑑定で、ある試料についてジェノタイプが実験的に確定していれば、その特定のジェノタイプ\(x_0\)が問題になります。
今、ある条件(その試料がある個人由来である、という仮説)を取れば、この確率は
\(
Pr(x_0|\mathbf{\theta})
\)
となります。
これは確率の式そのものですが、ジェノタイプとして観測し終えたものを想定しているので、呼び方を変えることにして、
> \(\mathbf{\theta}\)の下での尤度
と呼びます。
確率と尤度の関係については『確率と尤度は同じものの見方が違うだけ』の章も参考にしてください。
今、試料の由来主は２人のみに限られるとすると、その２人のそれぞれが$x$を観測する尤度としては
\(
Pr(x|\mathbf{\theta_1})
\)
\(
Pr(x|\mathbf{\theta_2})
\)
の２つが考えられます。
尤度が高い方が、「真の由来主」だろうと思えます。
その「真の由来主」らしさを表す方法の一つが尤度比です。
\(
\frac{Pr(x|\mathbf{\theta_1})}{Pr(x|\mathbf{\theta_2})}
\)
この式は、２人目に比べて１人目が真の由来主らしさは何倍かを数字にしたものです。
別の方法で「真の由来主」らしさを数字にする方法があります。
\(
\frac{Pr(x|\mathbf{\theta_1})}{Pr(x|\mathbf{\theta_1})+Pr(x|\mathbf{\theta_2})}
\)
とする方法です。
この式は尤度の相対的比率です。
この場合は最小で0、最大で1です。
尤度比は「●倍」で大小を表します。
差がなければ「１倍」です。
尤度の相対的比率は単位はありません。「0.9」なら0.9です。
差が無い時は「0.5」です。
どちらも数字が違うだけで、本来の意味は変わりません。
DNA試料の由来主が３人以上いる場合は、仮説も３つ以上になります。
その場合、尤度比はあくまでも、２つの仮説の尤度の比なので、\(n\)個の仮説に対して
\(\frac{n(n-1)}{2}\)個の尤度比が計算できます。
もし、尤度比を使って、ある特定の仮説の真偽を考えようとすれば、
この\(\frac{n(n-1)}{2}\)個の数字についての判断が必要になります。
簡単な考え方としては、\(\frac{n(n-1)}{2}\)個の尤度比の最小値でさえも
ある値より大きい、というような基準でしょう。
たくさんの尤度比を判断にあたってどのように使うかについてのコンセンサスは存在していないと思います。
一方、尤度の相対的比率であれば、
\(
\frac{Pr(x|\mathbf{\theta1})}{\sum_{i} Pr(x|\mathbf{\theta_i})}
\)
と、２仮説の比較の場合と同じ枠組みです。
DNA鑑定では２つの異なる立場からの主張が対立することがあると思います。
片方の主張(主張X)は、
> 『試料のジェノタイプはAさん以外の人に由来すると仮定すると、あまりにも尤度が低いので、Aさん以外の人に由来するという仮定は間違っていると思います。Aさん由来のはずです』
というものです。
もう片方の主張(主張Y)は、
>『試料がAさん由来である可能性もありますが、誰かひとりでもAさん以外の誰かがいて、その人由来であるという尤度が十分に低くなければ、Aさんであると断定することはできないはずです』
というものです。
主張Xでは、Aさんである、という仮説以外のすべての仮説について確率を出して、その総和確率(場合によっては事前確率で重みづけした総和)が十分に小さいことを示せばよいです。
他方、主張Yでは、たくさんの仮説があるなかで、ただ一つでも、Aさんとの違いが十分大きくない人がいることを示せばよいです。
主張Xでは総和を問題にしていますが、多くの場合は、現実的な計算上の理由から、平均値を使います。
主張Yでは、特定の１つだけについて問題にします。
もし、特定の１つの仮説だけは、「Aさんである」という仮説と十分な差がなく、でもそれは１つの仮説だけで、そのほかの大量の仮説は「Aさんである」という仮説と十分な差がある、という場合には、主張Xと主張Yとはどちらも成り立ってしまいます。
たくさんの仮説(DNA鑑定ではたくさんの候補者)があって、そのうちわけが不均一であるときに、このようなどっちつかずが起きます。
そのような不均一の最たるものは、一卵性双生児同胞の存在です。
したがって、仮説をどのように設定するかは個々の主張にも影響しますし、複数の主張でどっちつかずになるような場合にも大きく影響してきます。
前回は確率と尤度と仮説に関する内容でした。
また、尤度比は２つの仮説のもっともらしさを比較する指標であることも扱いました。
今回は、その尤度比を使った検定の話です。
同じく尤度比を使いますが、使い方が少し違うので、その違いを押さえることを目的とします。
前回書いた通り、ある仮説の下である事象が観察される確率を、観察事象が得られたあとで問題にして、仮説のもっともらしさを表す数値とみなしたとき、それを尤度と言います。
仮説のことをモデルと言うこともあります。
少し言葉のニュアンスが違うとすると、「モデル」という場合には、
変数を導入してその変数を使って確率・尤度が計算できる、という印象が強いかもしれません。
たとえば、「サイコロを３回振ったら、4,2,4が出た」という事象に対して、「振ったサイコロは理想的なサイコロである」というのは、「仮説」ですし、
「振ったサイコロは６つの目のどれも等確率で出る」と言うのも仮説です。
「振ったサイコロは1,2,3,4,5,6の目を持ち、それぞれの目が出る確率は\(p_i=\frac{1}{6};i=1,2,...,6\)である」と変数を持ちだして、確率的に起きる現象であることを説明してあるとモデルらしさが強いという感じでしょうか。
「サイコロを３回振ったら、4,2,4が出た」と言う事象は「サイコロを３回振ったら、1,2,3,4,5,6の目がそれぞれ、\(\mathbf{x}=(0,1,0,2,0,0)\)回ずつ出た」ともいいかえられます。
「６つのサイコロの目の出る確率が\( \mathbf{p}=(p_i);\sum_{i=1}^6 p_i=1;i=1,2,...,6\)」であるというモデルの尤度は
\(
L(\mathbf{ p }|\mathbf{ x })=\begin{pmatrix}\sum_{i=1}^6 x _ i\\x_1,x_2,...,x_6\end{pmatrix} \prod_{i=1}^6 p_i^{x_i}
\)
ただし、
\(
\begin{pmatrix}\sum_{i=1}^6 x _ i\\x_1,x_2,...,x_6\end{pmatrix} = \frac{(\sum_{i=1}^6 x _ i)!}{x_1!\times x_2! \times ...\times x_6!}
\)
となります。
この関数を尤度関数と言います。
\(p_i^{x_i}\)のようなべき乗があると計算が面倒臭いので、対数を取ってやります。
対数を取ったものを対数尤度関数と言います。
\(
LL( \mathbf{ p }|\mathbf{ x }) = \ln{ L(\mathbf{ p }|\mathbf{ x }) } = \ln{ \begin{pmatrix}\sum_{ i = 1 } ^ 6 x _ i\\
x _ 1,x _ 2 ,..., x _ 6\end{pmatrix}} + \sum_{ i = 1 } ^ 6 x _ i \ln{ p _ i }
\)
今、\(\mathbf{x}\)という観察結果があったときにどんな\(p_i\)が尤もらしいかを考えるとき、
尤度関数・対数尤度関数の
\(
\begin{pmatrix}\sum_{i=1}^6 x _ i\\x_1,x_2,...,x_6\end{pmatrix}
\)
の部分は\(p_i\)によらないので、どの目が出やすいかについての検討(\(p_i\)の検討)のときには、この項は無視できることになります。
すると、\(\mathbf{p_1}=(p_{1,1},...,p_{1,6})\)と\(\mathbf{p_2}=(p_{2,1},...,p_{2,6})\)とを比べるときには
、その対数尤度の差(尤度の比の対数を取ったものなので、対数尤度比と言います)は
\(
\sum_{ i = 1 } ^ 6 x _ i \ln{ p _ {1,i }} - \sum_{ i = 1 } ^ 6 x _ i \ln{ p _ {2,i }} = \sum_{i=1}^6 x_i(\ln{\frac{p_{1,i}}{p_{2,i}}})
\)
のように単純な式になります。
ここで、サイコロの目の出る確率\(\mathbf{p}\)について、２つのモデルを考えます。
一つ目のモデルは、「理想的なサイコロであるモデル」で、\(\mathbf{p^1}=(\frac{1}{6},\frac{1}{6},...,\frac{1}{6})\)です。
二つ目のモデルは、なんでもいいから適当に決めてよいモデル\(\mathbf{p^2} = (p_i)\)です。
一つ目のモデルは6つの確率の値\(p_i\)のどれも、自由に決めることができませんでした。
二つ目のモデルは、\(\mathbf{p^2}\)の6つの確率の値のうち\(p_1,...p_5\)までを自由に(ただし\(p_1,...,p_5 \ge 0\)で、\(p_6 =1-\sum_{i=1}^5 p_i \ge 0\) という制約はありますが、その範囲では自由に)決めることができました。
かたや、まったく自由がなく、かたや、５つが自由でした。
この自由に値を選べる変数の数の差を自由度と呼びます。
二つ目のモデルは一つ目のモデルに比べて、自由度が５だけ大きいのです。
この自由度は後に出てくる尤度比検定のところで使います。
たとえば\(\mathbf{p^2}=(p_1,p_2,...,p_6) = (0,\frac{1}{3},0,\frac{2}{3},0,0)\)とするのも自由です。
これは、二つ目のモデルの尤度が最も大きくなるように選びました。
このように尤度を最大にするモデル変数の値を最尤推定値と呼びます。
最尤推定値を尤度関数の微分で出すことについては『最尤推定値 対数尤度関数の微分』の章を参照してください。
２つのモデルがあって、両者の違いは自由な変数の個数であるとき、変数が自由なモデルの方が
尤度が大きく、従って対数尤度も大きくなります。
２つのモデルの対数尤度の差(２つのモデルの尤度の比の対数を取ったものなので、対数尤度比と言います)の２倍の値が有用になります。
\(
2\ln{\frac{L(\mathbf{p^2}|\mathbf{x})}{L(\mathbf{p^1}|\mathbf{x})}}=2(LL(\mathbf{p^2}|\mathbf{x})-LL(\mathbf{p^1}|\mathbf{x}))
\)
この統計量は帰無仮説(サイコロが理想的であるという仮説)が成り立っているときには、自由な変数の数を自由度としたカイ二乗分布に従うことが知られています。
従って、この統計量を用いて帰無仮説を棄却すべきかどうかの検定(棄却検定)が可能となります。
Rで検定してみます。
---
```{r fig.width=7, fig.height=6}
p1 <- rep(1/6,6)
p2 <- c(0,1/3,0,2/3,0,0)
x <- c(0,1,0,2,0,0)
LL1 <- lgamma(sum(x)+1)-sum(lgamma(x+1)) + sum(x*log(p1))
LL2 <- lgamma(sum(x)+1)-sum(lgamma(x+1)) + sum(x[which(x!=0)]*log(p2[which(x!=0)]))
S <- 2 * (LL2-LL1)
df <- 5
pchisq(S,df,lower.tail=FALSE)
```
---
帰無仮説が成り立っているときに尤度比検定を繰り返すと、p値が一様分布する様子は、『尤度比検定〜サイコロの例〜』章を参照してください。
尤度比検定は「尤度比」という言葉が入っているけれども、多くのDNA鑑定の場合のように、２つの仮説の尤度の比を比べているときには用いません。
なぜなら、DNA鑑定での２つの尤度は、それぞれ「決め打ちにした仮説」であるのに対して、
尤度比検定の方の尤度は、かたや、何かしらモデル変数の自由を奪って制約したもの、
かたや、モデル変数を自由に動かしたものです。
従って、対立仮説の方の尤度は、自由に変数の値を動かせる世界の中で、最大の尤度になっています。
DNA鑑定の尤度比検定では「もう１つの仮説」の尤度が、何かしらの自由な世界の頂点であるというわけではない点が違います。
実際、尤度比検定で統計量をカイ二乗統計量に照らして評価するのは、変数を自由に動かす世界に連続的充満している仮説群があることを前提にしたいるから可能であるのです。
では、この尤度比検定はDNA鑑定では用いないのか、というと、そうではありません。
実験データの評価をするときに変数を用いたモデルを立てたり、ベイジアンネットワークでモデル変数を用いたりするときには、立てたモデルが意味のあるものであるのかどうかの評価に尤度比検定の考え方を使います。
サイコロの例でもみたとおり、変数を自由にすると尤度は大きくなりますが、「変数を自由にしたのなら、それに見合うだけ十分に大きくならないのなら、その尤度の増加は意味のないものかもしれない。意味のありなしを尤度比検定の考え方で評価しよう」という使い方です。
2011年3月11日、東日本大震災とそれに引き続く大津波があり、多数の方が亡くなられました。
私ごとながら、出生地である宮城県多賀城市も大きな被害を受けました。
さて、この大震災により、多人数の行方不明者と多数のご遺体が県境を越える形で身元鑑定を必要とする事態が発生しました。
DNA鑑定はその強力なツールの一つとなりました。
実務上、いくつかの課題があると法医学講座より情報提供を受け、それにそって取り組みました。
まずは、血縁者からのDNA提供が得られている状況での行方不明者のジェノタイプ推定と
それに基づく尤度計算です。
また、複数の行方不明者と複数のご遺体との多対多マッチング問題がありました。
また、ジェノタイプ情報はそれ以外の情報(ご遺体の発見場所・身体的特徴)などと併せて解釈されるべきであることも問題となりました。
これらについて、何回かに分けて法数学勉強会でも取り上げました。
本章はそのうち、血縁情報を用いたDNA鑑定用の尤度計算法に関してになります。
勉強会で使用したスライドは[こちら](http://www.genome.med.kyoto-u.ac.jp/StatGenet/RY%e8%ac%9b%e7%be%a9%e9%96%a2%e4%bf%82/2011/%e5%ae%b6%e7%b3%bb%e3%82%b8%e3%82%a7%e3%83%8e%e3%82%bf%e3%82%a4%e3%83%97%e6%83%85%e5%a0%b1%e3%81%ae%e7%a2%ba%e7%8e%87%e5%b0%a4%e5%ba%a6%e8%a8%88%e7%ae%9720110520.pptx こちら)です。
ある家系図があって、その家系がある集団に属しているとき、家系図上の
メンバーのジェノタイプを活用して、ジェノタイプ不明なメンバーの
ジェノタイプの確率を計算するとは、どのような手順で行われているのでしょうか？
一つのやり方を示します。
- 家系図を描く
- ジェノタイプを提供しているすべての個人をプロットする
- 親子関係をつなぐ
- 片親しかいない場合には、もう一人の片親を仮に描いてつなぐ
- 親子以外の血縁関係も途中にはいるべき血縁者を仮に描いておく
- 親子関係をすべて列挙する
- 親子関係は両親と子のトリオ
- ジェノタイプのない個人もすべて「仮」に描いてあるので、すべての人は次の二通りに分かれる
- 親が２人とも描かれている
- 親がどちらも描かれていない
- 親がどちらも描かれてない人は、「集団から生まれた」とみなす
- 個人のジェノタイプを確率的に決める
- ジェノタイプが提供されている人の場合は、１つのジェノタイプの確率が１
- 「集団から生まれた人」のジェノタイプは集団のジェノタイプの頻度に応じて確率的に決まる
- それ以外の人(両親がともに家系図に描かれているが、本人のジェノタイプは提供されていない人)のジェノタイプは、確率的に計算する
問題は「それ以外の人(両親がともに家系図に描かれているが、本人のジェノタイプは提供されていない人)のジェノタイプは、確率的に計算する」を実際にどうするか。
それについては細かい話になるので、『近親婚のない家系図でのジェノタイプ確率計算』の章を参照してください。
ただし、ジェノタイプの確率的推定方法に関して、家系図は２種類に分かれることは有用な知識と思いますので、コメントしておきます。
家系内に近親婚がある場合、これは、家系図をグラフで描いたときに閉じた輪がある場合です。
- - -
```{r fig.width=7, fig.height=6}
library(kinship2)
test1 <- data.frame(id =c(1,2,3,4,5,6,7,8,9,10,11), mom =c(0, 0, 2, 2,2, 2, 4, 10,8,0,0), dad =c(0, 0, 1, 1, 1, 1, 11,6,7,0,0),sex =c(0, 1, 0, 1, 1, 0, 0, 1, 0, 1,0))
affected<-rep(0,11)
affected[9] <- 1
ptemp<-pedigree(id=test1$id,dadid=test1$dad,momid=test1$mom,sex=test1$sex,affected=affected)
plot(ptemp)
```
- - -
図では７と８の間の関係がいとこ婚です。
『９→７→４→(１，２）→６→８→９』
という輪ができていて、これが家系内の近親婚に対する輪です。
ジェノタイプの確率計算においては、このような『輪』がない場合には、簡単ですが、
『輪』があると面倒臭くなります。
これはグラフ(家系図は個人(や個人が持つ染色体)を点とし、伝達関係を辺としたグラフ)を扱う理論であるグラフ理論において、サイクル(ぐるりとめぐる『輪』)があるかないかで、アルゴリズムが大きく変わることに対応しています。
グラフ理論はDNA鑑定領域であれば、ベイジアン・ネットワークもグラフです。
ベイジアン・ネットワークでもグラフのサイクルの有無は問題にされます。
本章は、前章に引き続き大震災関連課題を取り扱います。
多人数の行方不明者と多人数のご遺体とのマッチングに関する話です。
勉強会で使用したスライドは[こちら](href{http://www.genome.med.kyoto-u.ac.jp/StatGenet/lectures/2010/HouSuugaku/taninzuuprofile.pptx こちら)です。
ある行方不明者\(m\)のジェノタイプがわかっており、ある遺体\(b\)のジェノタイプがわかったとき、本人確認は、それが一致するかどうかを基にします。
通常の刑事上のDNA鑑定とだいたい同じ枠組みです。
\(m\)のジェノタイプがわかっておらず、その血縁者のジェノタイプが得られたとき、
\(m\)のジェノタイプは一意に決まらず、確率的に決まりますが、
それは計算できます。
もし、\(m\)が１人で、\(b\)が１体であったなら、
- \(m\)は\(b\)である
- \(m\)以外の誰かが\(b\)であって、\(m\)は\(b\)以外の誰かである
という２つを問題にすることになります。
複数の人\(M=\{m_1,m_2,...,m_{Nm}\}\)が行方不明や犠牲になり、
複数のご遺体\(B=\{b_1,b_2,...,b_{Nb}\}\)が見つかった場合というのもある。
今、行方不明者数とご遺体数が同数であるときを考えます。
考えるべき仮説
- \(m_1 = b_1,m_2 = b_2,m_3 = b_3,...,m_{Nm} = b_{Nm}\)
- \(m_1 = b_2,m_2 = b_1,m_3 =b_3,...,m_{Nm} = b_{Nm}\)
- ...
となりますが、この仮説の数は
\(
N_m ! = N_m \times (N_m-1) \times (N_m-2) \times ...\times 2 \times 1
\)
となります。
それぞれの仮説の下で観察データ(この場合は、\(M\)のジェノタイプ(血縁者試料からの推定ジェノタイプを含む)と\(B\)のジェノタイプのこと)を観察する尤度が計算できます。
原理としては、これを比較すればよいわけです。
複数人について$N!$通りの仮説を考えてそれぞれの尤度を計算して、
尤度が最大のものを選び出す、というのは、言うのは簡単ですが、
\(N\)の数が少し多くなると、\(N!\)がものすごく大きくなるので、非現実的です。
- - -
```{r}
N <- 1:20
factorial(N)
```
- - -
\(N=5\)で120通り、\(N=10\)で300万通りを超え、\(N=12\)で約５億通りです。
しかしながら、何人かが集まって、その人たちをうまく割り当てたい、という希望は、
DNA鑑定に限らず、日常茶飯なことなので、
この手の課題は「割り付け問題」と呼ばれ、比較的うまく対処する方法があります。
たとえば、\(N!\)通りのすべてについてしらみつぶしに調べることは無理だけれど、
最もよい割り付け方を見つけ出す、というアルゴリズムはあって、
コンピュータを使えばわかります。
『ハンガリアン法による最適割り付け』の章に簡単に計算できる様子を
示しました。
最もよい割り付け(最適割り付け)を\(N!\)通りから見つけ出せたら、それが「正解」なのでしょうか？
最適割り付けは「尤度」が
>「最も高い割り付け方である」
ということしか言っておらず、
>「他の割り付け方はあり得ない」
と言っているわけではありません。
もしも、２番目に高い尤度が、最大尤度よりも明らかに小さいことがわかれば、安心して
最大尤度をもたらす割り付けを採用することができます。
最大尤度をもたらす割り付けを見つけるアルゴリズムはあると書きましたが、
上から\(k=1,2,...\)番目までの尤度をもたらす
割り付けまでを見出すアルゴリズムもあります。
もちろん\(k\)がそれほど大きくない場合です。
これができるのであれば、\(k=1,2\)を調べ、その２つの尤度に大差があれば、
選ぶべき割り付けは決まります。
もし\(k=1,2\)の割り付けを見つけ、\(k=1\)を選ぶことが適切であると決まらなかったら
どうしたらよいでしょうか？
決まらない、と言うのは、
>\(k=1,2\)の尤度がそれほど違わない場合で、
\(k=1\)の方が尤もらしいけれど、\(k=2\)である可能性を捨てきれない
という場合です。
すべての行方不明者とすべてのご遺体との対応が一発でわかれば、とても素晴らしいですが、ある特定の行方不明者を探している人の立場からすれば、
- 私の探している人はこのご遺体なのか、そうでないのか
- 私の探している人はこのご遺体であるかもしれない、と言うことだが、それ以外のご遺体である可能性はないのか
というように、もっぱら、自分の探している人について興味があるだけで、
それ以外の行方不明者とご遺体とのマッチング関係については
興味がほとんどないでしょう。
また、この質問に答えることが、ご遺体を遺族が引き取ることができるかの
決断のためになすべきことであることも多いでしょう。
では、この個別の希望に対して、割り付け問題はどのようにアプローチすることができるでしょうか？
今\(N\)人の行方不明者とご遺体があるときに\(m_i\)いう行方不明者についての
質問に答えたいとします。
\(m_i=b_1,m_i=b_2,...,m_i=b_N\)という\(N\)通りの仮説があります。
割り付け問題では\(N!\)通りの仮説があったのに比べると、
ずいぶん少ない数です。
ただし、\(m_i=b_j\)という仮説は、実際には
\(m_i\)以外の\(N-1\)人の行方不明者と\(b_j\)以外の\(N_1\)体のご遺体とが
どのように対応づけられるかのすべてが合算されています。
\(N-1\)人と\(N-1\)体の総当たり割り付け問題が封じ込められているわけです。
ここには\((N-1)!\)通りの割り付け方法がありますから、この確率・尤度をすべて計算する
ことにしたのでは、場合の数が大きすぎる問題の
解決にはなりません。
そうではなくて\((N-1)!\)通りの確率・尤度を個別に計算して
足し合わせる代わりに、すべての確率・尤度の和をいっぺんに計算することが
できれば、場合の数問題に対応したことになります。
そのような方法は規模がそれほど大きくなければないわけではないようです。
そんなアプローチが適切であるかもしれません。
棄却検定では正確確率検定と呼ばれる方法があります。
本章では、得られた情報から分割表を作成し、正確確率を計算することで、ある容疑者が犯人である確率を正確に計算することについて考えます。
勉強会で使用したスライドは[こちら](href{http://www.genome.med.kyoto-u.ac.jp/StatGenet/lectures/2011/Housuugakuwadai1-2012-0324.pptx こちら)です。
>ある小学校では、すべての男子生徒３００人に「消しゴム屋」さんから１個ずつ消しゴムが配られた。
>「ナルト」柄が５０個、「ワンピース」柄が１００個、「ドラえもん」柄が１５０個であった。
>ある日、音楽室で、１個の消しゴムの落し物が発見された。
>「ナルト」柄だった（５０個/３００個）
>その消しゴムには
>>「ほなみたまちゃんLOVE」
>と書かれていた！
>さあ、この消しゴムは誰のものか？
>女子児童たちの捜査が始まった。
>捜査員の一人であるみぎわさんが消しゴム屋に迫った。
>>「消しゴム屋さん、『ナルト』柄の消しゴムを配った男子のリストを出しなさいよー」と。
>リストが得られれば、『ナルト』所有者50名に絞られるからだ。
>そのリストに載っている50名は、1/50の確率で「たまちゃんLOVE」であるとわかるはずだ。
>しかし、その情報は得られない。
>>「あいにく、個人情報保護法の規定により、リストをお出しすることはできません」
>と、消しゴム屋さん。
>ただし、
>>「特定の児童さんを指名していただけば、その児童さんにどのタイプを渡したかはお答えできます」と。
> ...
>>「たまちゃんLOVE」の『犯人』って誰だと思う？
>>「やっぱり、同じクラスの男子？」
>>「クラブが一緒？」
>>「委員会つながりかも？」
>>「以前に同じクラスだった・・・とか」
>>「帰り道が一緒だからだったりして」
>と、女子全体の詮索が続く。
>それを聞いた、まるちゃん母、
>>「事前確率なんて考えるのはやめなさい、下手な考え、休むに似たりって言うじゃない」
>そうとばかりは言えないのだが…
>そんなところへ、音楽の先生から耳より情報が得られた。
>>「音楽室はいつも鍵がかかっていて、鍵は私しか持っていないから、その日、音楽室に出入りできた児童は、その日に音楽の授業があったクラスの子、だけよ」
>と。
>>「その日、音楽の授業があったのは、３クラス、１２０人のうち、男子は、６０人」
>>「３００人のうち、６０人が候補者だ」
>話を整理しよう。
- 「その日、音楽の授業があったのは、３クラス、１２０人のうち、男子は、６０人」
- 「音楽の授業は、３００人のうち、６０人」
- 「『ナルト』柄は、３００人のうち、５０人」
-「じゃあ、音楽のあったクラスで、『ナルト』消しゴムを持っていたのは何人？」
>ここから得られるのは
>>「６０ｘ５０/３００ ＝１０人 くらいね！」
>さて、消しゴム屋さんからの情報も使うとしよう。
>>「特定の児童さんを指名していただけば、その児童さんにどのタイプを渡したかはお答えできます」と。
>みぎわさんは
>>『絶対、花輪君のを教えてもらうわ！誰も文句ないわよね！』
>と方針を決定。
>それを受けて、消しゴム屋さんいわく
>>「花輪君ですね、花輪君は、『ナルト』柄でした」
>>みぎわ：「!弩?恐狂!!逆散魔*警鱗!鬱+閉静無..止....」
>「花輪君以外の男子(５９人)が、『ナルト』柄である確率は 50/300 = 1/6→49/299」
>「そのうちの誰かが、落とし主だとすると、それが『ナルト』柄である確率は49/299」
>「花輪君が落としたのなら、それが『ナルト』柄である確率は 1」
>「誰かが落としたのなら、それが『ナルト』柄である確率は 49/299」
>「尤度比は1/(49/299) = 299/49 = 6.10」
>「花輪君が落としたと考える方がよさそうじゃない？」
>>みぎわ：「叫怒天髪射」
>「候補の男子が59人もいるのに、いきなり、『 6.10倍』っていうのは、おかしいんじゃねーか？」
>「花輪君が落としたのなら、それが『ナルト』柄である確率は 1」
>「誰か１が落としたのなら、それが『ナルト』柄である確率は 49/299」
>「誰か２が落としたのなら、それが『ナルト』柄である確率は 49/299 」
>…
>「誰か59が落としたのなら、それが『ナルト』柄である確率は 49/299 」
>「花輪君 vs. その他59人 は 1: 59x49/299 = 0.103 」
>>みぎわ：「莞喜爾天静昇」
> こうだったじゃろう
>>「音楽の授業は、３００人のうち、６０人」
>>「『ナルト』柄は、３００人のうち、５０人」
>>「じゃあ、音楽のあったクラスで、『ナルト』消しゴムを持っていたのは何人？」
>>「６０ｘ５０/３００ ＝１０人 くらいね！」
>「本当のところは、音楽授業の６０人中、何人が『ナルト』柄じゃったんかいのー」
>１人、２人、…、５０人の場合が全部ありえるはずじゃなぁ
>１人だけが『ナルト』柄で、それが花輪君なら、「たまちゃんLOVE」は花輪君で確定じゃ
>２人が『ナルト』柄なら、花輪君かもしれないし、もう１人かもしれなくて、確率は1/2
…
>k 人が『ナルト』柄で、花輪君がそのうちの１人なら、確率は1/k
>これを計算するとどうなるかの？
『ナルト』柄の確率は
\(
p=\frac{50}{300}=\frac{1}{6}
\)
音楽室に出入りした50人中、$k$人が『ナルト』柄である確率は
\(
pr(k) = \frac{50!250!60!240!}{300!k!(50-k)!(60-k)!(190+k)!}
\)
こんな$2\times 2$表で考えれば計算式はわかりやすい。
|条件 |ナルト | ナルト以外 | 行和 |
|:-----------|------------:|:------------:|:------------:|
| 音楽室に居た | k | 50-k |50 |
|音楽室に居なかった| 60-k | 190+k | 250 |
|列和 | 60 | 240 | 300 |
$k$の大きさに比例して『ナルト』柄が落ちやすくなるから
\(
Pr(k) = k \times pr(k) = k\times \frac{50!250!60!240!}{300!k!(50-k)!(60-k)!(190+k)!}
\)
$k$人の誰もが「たまちゃんLOVE」と書く確率は同じであれば
\(
Hanawa(k) = \frac{1}{k} \times k \times pr(k) = \frac{50!250!60!240!}{300!k!(50-k)!(60-k)!(190+k)!}
\)
$k=0$の場合は『ナルト』柄消しゴムが見つかっているから、
もうありえないので、
すべての場合の和は
\(
\sum_{i=1}^{50} Pr(k) = \sum_{i=1}^{50}k\times \frac{50!250!60!240!}{300!k!(50-k)!(60-k)!(190+k)!}
\)
花輪君が「たまちゃんLOVE」であるのは
\(
\sum_{i=1}^{50} Hanawa(k) = \sum_{i=1}^{50}\frac{50!250!60!240!}{300!k!(50-k)!(60-k)!(190+k)!}
\)
この方法では
>花輪君が「たまちゃんLOVE」である尤度は，0.100
>花輪君が「たまちゃんLOVE」である尤度と、他の誰かが「たまちゃんLOVE」である尤度の比は、0.111
場合によっては作成される分割表が複雑になります。
そのように複雑になった分割表では、周辺度数による分割表制約がベイズ流の事後確率としての性質を強くすることもわかりました。
これまではジェノタイプ情報やその他の情報を用いてある仮説がどれくらい尤もらしいかを考えてきました。
実際のところ、この尤もらしさの評価が提供しているのは尤度比であって、それは、事前確率を事後確率に変化させる要素であって、事前確率が不明なときには、利用の意味が大きく損なわれます。
従って、DNA鑑定においてジェノタイプ情報を的確に活用するためには
ジェノタイプがもたらす尤度比のみではなく、その一歩手前の事前確率に関して、
どのようなことを統計学・数学が提供できるかを考えることは重要です。
本章はこの点に関する話です。
事前確率は情報がない状態もしくは、情報が限られている状態で思い描く、仮説の確率のことです。
例でそのことを確認します。
ある夏の日の最高気温が何度になるかを予想するとします。
最高気温ですから、20-40度くらいの何かしらの値になるでしょう。
今、これ以外に情報がないとき、
>「20-40度」くらい
という気持ちが事前確率です。
ただし、20度の確率はどれくらいで、30度の確率はどれくらいで35度の確率はどれくらいで40度の確率はどれくらいなのか、
32.3度の確率はどうなのか、と連続した実数値に対して、「思い入れ」の強弱のカーブが描けるはずです。
可能性の話をしますから、絶対零度から無限大までという取りうる温度の数値のすべてに対して、「思い入れ」の強弱をつけ、それを足し合わせ(積分)したときに１になるように、スケールを調整したもの、それが事前分布です。
朝、６時半、外から元気に体操する声が聞こえたとします。
その情報から、雨が降っていないようだし、『朝、晴れていた夏の日の最高気温』になりそうだと考えて、少し情報修正するかもしれません。
このように情報修正したとすれば、「体操」という情報を使って、
事前確率の分布を事後確率の分布に変えています。
このようにして得た事後確率の分布も、次に「クマゼミがいっせいに『シャワシャワ』と鳴く声が聞こえた」という情報を受けて、さらに上方修正されるかもしれません。
ですから、事後確率はさらなる情報に関して事前確率として働きます。
このように情報に応じて変化する事前確率の分布ですが、どのような分布を描くかは、予想する人の経験等にも左右されますから、個人個人で異なります。
心の中にある事前確率分布、と書いてきましたが、予想するべき場所と日付がわかれば、過去の気象記録を活用することは、いい手でしょう。
天気などについてまったくわからなければ、天気などによらず、過去の最高気温分布をそのまま使うのが良さそうです。
もし、予想したい日が朝から晩まで晴れと思われているなら、過去にあったそんな日の最高気温だけを使えば良いでしょう。
孤島にサルの群れが暮らしていて、あるメス猿から子ザルが生まれたとします。
この子ザルの父親はどのオス猿なのか、ということを知りたいとします。
サル山にボス猿がいれば、それが父親である確率は相当高いことでしょう。
どれくらい高いかは、「心の中」にあるかもしれないですし、
何かしらの調査の結果に基づいて、ある程度の絞り込みができるかもしれません。
もし、この子ザルがちょうど「ボス猿不在の時代」の生まれだとします。
どのように考えればよいでしょうか。
相変わらず「心の中」に分布はありますが、もし、母猿とオス猿すべての接触時間に関する情報が得られるのであれば、それを客観的に活用して、各オス猿について「父親である確率」の高低を調整することができるでしょう。
シンデレラが落として行ったガラスの靴を頼りに、王子がシンデレラを探すのも、靴を用いて、女性の事前分布を変えて、それらしい女性に絞り込むためです。
同じようなことはそこここでやられています。
ただし、それを行う時に、「シンデレラではありえない」「シンデレラであるかもしれない」の２つにはっきりと分けることができるか、最高気温予想のように、どの値がどれくらいありそうかを相対的に重みづけするかでは、少し考え方が違うでしょう。
DNA鑑定の場合には、多くの場合、はっきり２分するというよりは、色々な仮説の事前確率が高かったり低かったり色々だけれど、0であることはない、としておくことが重要です。
しかしながら、DNA鑑定を用いて、判決という決断を下さなければならないとき、「決断」は0か1かの選択をすることですから、どんなに頑張っても、ある程度の曖昧さを切り捨てる「痛み」が生じることはやむをえないと考えることが適当です。
事前分布のことを扱ったので、共役事前分布にも触れておきます。
共役事前分布はベイジアン・ネットワークなどでよく使うものです。
例を出します。
画鋲を投げたときに、針が上を向いて止まる場合と針が下を向いて止まる場合とがあります。
今、画鋲が針を上に向けて止まる確率$p$、下に向けて止まる確率$1-p$について知りたいとします。
画鋲を投げる実験を開始する前に、$p$の予測を「心」に聞いてみます。
それを実験前の予測分布として$p_0$とします。
$p_0$は$0$から$1$の値を取る分布であるはずです。
どんな分布を思い描いてもよいのですが、「ベータ分布」というのを思い描くと便利ですよ、というのが、「共役事前分布」とは何か、という話とつながります。
画鋲の裏(針が下を向いている状態)が出る確率が0である確率はないだろうし、1である確率もないだろう。
0.5ということもなさそうだ。
0.5よりは大きいところにピークを持たせて、$p=0,1$で0になるような分布が描きたい。
ベータ分布と呼ばれる分布は次の図のような形をしており、
$0$から$1$の間で積分すると$1$になるという性質がある。
形と言い、積分すると$1$になる性質といい、事前分布としてとてもよい。
```{r fig.width=7, fig.height=6}
p <- seq(from=0,to=1,length=100)
b <- dbeta(p,4,3)
plot(p,b,type="l")
abline(h=0)
```
実験を開始しました。
$n$回投げたら、$k$回の裏、$n-k$回の表が出るでしょう。
裏の出る確率を$p$としたらこの確率は
\(
\begin{pmatrix}n\\k\end{pmatrix} p^k (1-p)^{n-k}
\)
となります。
これが二項分布です。
$n=5,k=0,1,2,3,4,5,p=0.6$の場合のそれを見てみます。
```{r fig.width=7, fig.height=6}
n <- 5
k <- 0:5
p <- 0.6
a <- dbinom(k,n,p)
plot(k,a,type="h")
```
$p$の値を決めると二項分布が描けました。
$p$を決めるというのは、仮説を決めたことになり、二項分布を描いたのは、
その仮説の下での生起確率を計算したことになります。
生起確率があれば、尤度があることは冒頭の章でも出てきた内容ですから、それについて
考えてみます。
尤度とは、ある仮説が、ある観察結果をもたらす確率のことですが、その確率を仮説の側からみたものです。
ある特定の仮説ではなく、すべての仮説について、ある特定の観察結果をもたらす確率がわかるとよいです。
画鋲の表裏の場合の仮説とは$p$が$0$から$1$のどれかです。
$p$は連続値ですから、$0-1$の範囲で連続関数のはずです。
実は、このような関数がベータ分布です。
実際、その計算式は
\(
\begin{pmatrix}n\\k\end{pmatrix} p^k (1-p)^{n-k}
\)
と同様です。
ただし、$p$について$0-1$で積分をして1にならないといけませんから、
\(
\frac{p^k (1-p)^{n-k}}{\int_0^1 t^k (1-t)^{n-k} dt}
\)
となります。
今、画鋲を実際に$n$回投げて$k$回の裏$n-k$回の表が出た、という情報があれば、$p$の尤度はベータ関数で計算できるから、それをそっくりそのまま使えば良いです。
しかしながら、昔、画鋲を投げたとき、裏が$K$回、表が$N-K$回くらいだったな、と言うような、おぼろげな記憶があったとしたら、どうしたらよいでしょうか？
ベータ分布で実際に投げた$n,k$を使って
\(
\frac{p^k (1-p)^{n-k}}{\int_0^1 t^k (1-t)^{n-k} dt}
\)
とする代わりに、実験前の事前分布を
\(
\frac{p^{K} (1-p)^{N-K}}{\int_0^1 t^K (1-t)^{N-K} dt}
\)
と想定してやれば、
その後、実際に実験をした後に信じるべきなのは
\(
\frac{p^{K+k} (1-p)^{(N+n)-(K+k)}}{\int_0^1 t^{K+k} (1-t)^{(N+n)-(K+k)} dt}
\)
と簡単に指数を変えるだけで済みます。
また、実験結果だけを用いて、$p$の値を推定するよりも、
実験する回数を減らしてもよくなるかもしれません。
逆に言うと、実験結果だけを用いて$p$の分布を推定する、というのは、事前確率として、おぼろげな記憶がない、つまり$N=0,K=0$という状態からスタートして予想する、ということでもあります。
たしかに$N=0,K=0$のベータ分布は次に示すように$0-1$のどれも等確率になった分布のことです。
```{r fig.width=7, fig.height=6}
p <- seq(from=0,to=1,length=100)
b <- dbeta(p,1,1)
plot(p,b,type="l")
```
ベータ分布は二項分布の共役事前分布(二項分布を観測の分布としたとき、仮説の側の分布としてとるとうまく行く分布)はベータ分布であるという話でしたが、
多項分布のそれはディリクレ分布、正規分布のそれはある場合は正規分布、ある場合はガンマ分布、というようにいくつかのものが知られています。
DNA鑑定において、ジェノタイプを調べて尤度比を求めたりするというのは、ある人が仮説を持っていて、その事前確率に差が小さいために決断できないときに、尤度比を提供することで、
事後確率に変化させ、決断しやすくするためにやっていると言えます。
この事前確率→情報提供→事後確率→事後確率に基づく決断という枠組みにこの章は基づいています。
勉強会で使用したスライドは[こちら](href{http://www.genome.med.kyoto-u.ac.jp/StatGenet/lectures/2013/Houigakkai2013June.pptx こちら)です。
複数のオプションがあって、どちらかを選ばなくてはならないことはたくさんあります。
服を選ぶ、レストランのメニューから料理を選ぶ…。
ものすごく悩ましくて決められないこともあれば、即決即断できることもあります。
また、同じ条件が与えられても、すぐに決める人もいれば、なかなか決めない人もいます。
ここで｢同じ条件」なのに個人差がある、ということが重要です。
DNA鑑定では、決断に資する情報を提供します。この情報は決断のための条件です。
情報が同じ・条件が同じでも、個人によって選択行動は違います。
では、情報以外に決断に影響している要因は何なのでしょうか？
このような問に答える、もしくは答えようとする理論が決断理論です。
決断理論には、哲学、心理学、経済学、数学などが絡みます。
DNA鑑定とそこから出てくる尤度や確率の情報は決断理論の中の数学(の中の統計学)に含まれます。
仮説間の尤度比が何倍だったら、決断できるか、というのは、心理学や哲学的な色彩を帯びてきます。
特にその決断が「重大」なものであるとき、心理学・哲学は大きな意味を持つかもしれません。そもそも「決断が重大であるとは何か」というような哲学的な疑問も生じてきます。
生物は太古の昔から、餌を探しにあっちに行こうかこっちに行こうかというような決断に迫られたり、繁殖戦略として有性生殖をしようか無性生殖にしようか、卵生にしようか胎生にしようか、というような決断・選択を繰り返してきました。
そのような過程をへて、現在、生き延びている生物種は、「長い目で見たときに生き残る」戦略をとってきたものたちであろう、と考えられています。
この「長い目でみて、サバイバルレースに勝ち残る戦略」として、情報科学が教える最適戦略は、「情報がないならないなりに、あるならあるなりに選択・決断すべし」「選択・決断にあたって、よさそうなものを選ぶのがよいが、少しは、よくないかもしれないものも選んでおくがよい」というものだと言うことです。
このように選択・決断にばらつきを持たせるのがよいという戦略家の末裔である人間が、同じ条件を与えられたときに、決断にばらつきが出るのは、妥当かもしれません。
そんな視点から考えると、DNA鑑定で提供するべきは、数学的に正確な情報であり、決断するべき人が使いやすい形で提供することが大事なのであって、その結果、各個人がどのような決断をするかには、『完璧な』介入はできないと考えることも適当かもしれません。
人が決断するとき、
- 事前分布
- 情報
- 事後分布
の３要素があります。
法数学はこの３つのポイントのそれぞれで、役割を持つことができると思います。
以下、ポイントごとに考えてみます。
決断するべき人の決断支援をするとして、仮説の尤度の違いとはどういうことかなどを理解可能な形式で説明することに尽きるように思います。
判断については無色透明を目指すことが数学的には妥当かと思います。
事前分布は「心の中」にあるという話が前の章で出てきました。
心の中の事前分布に客観性を持たせるために情報を的確に活用する、という話もその際にしました。
情報を的確に利用する、というのは言うのは簡単ですが、実行するのは非常に難しいです。
ですから、難しいから、「心の中」の事前分布が個人個人で違っていてもよいことを認めた上で、その個人差を小さくする情報活用の仕方についての検討は重要であると思います。
具体的には、考慮に値する仮説の広さについて適切な推定をし、それを提示することは重要だと思われます。
また、ジェノタイプに基づく尤度・尤度比等が提供されたとき、事前分布が変われば、同じジェノタイプ情報でも事後分布は変わります。
個人というのは、電卓ではないので、次のようにはできません。
>先日はこういう事前分布だった、そのときにジェノタイプ情報を使うと、これこれの事後分布になりますよ、と教えました。
>事前分布が変わったって？じゃあ、先日提供した尤度比を新しい事前分布に掛けてください、新しい事後分布が出ますから
この計算ができるわけではないので、要請があるたびに、事前分布に違いが生じるたびに、ジェノタイプ情報を適用しなおして事後分布を解釈しやすい無色透明な説明で提示することが適当だろうと思います。
DNA鑑定に関する多くの数学的・統計学的テーマがこの部分に集中して当てはまります。
特に気になることを挙げてみます。
取り扱いたい仮説が不特定多数を含む場合に、不特定多数の平均値や代表値を用いることが適切な場合があります。
特に、尤度の期待値を計算するときなどは、それが適切であることも多いです。
しかしながら、尤度の信頼区間をもとめたり、最大値・最小値を求めたりする場合には、不特定多数が作っている分布自体を考慮しないとわからなくなります。
極端な例としては、疑われている集団が赤の他人の場合とその中に血縁者が含まれている場合とがあると、集団の平均を用いることで失われる情報は後者でより大きくなる、というような例です。
そのほか、この範疇の課題として、勉強会にて教えていただいたものを以下に列挙します。
- 試料が複数人の混合である場合
- 試料が希少で曖昧な実験結果しか得られない場合
- 実験結果のクオリティを考慮しなくてはならない場合
- 複数の試料がえられるも、整合性に欠ける場合
- 集団に関する情報が限局的である問題
ひとまず、本文書は、ここまでで一区切りとします！
成功率が$0 \le p \le 1$であるときに、20回試行して、$k=0,1,2,...,20$回成功する確率は
\(
\begin{pmatrix}20\\k\end{pmatrix} p^k(1-p)^{20-k}
\)
です。
これを$p=0,0.01,0.02,...,0.99,1$について計算してプロットしてみます。
```{r setup}
library(rgl)
knit_hooks$set(rgl = hook_rgl)
```
```{r}
p <- seq(from=0,to=1,by=0.01)
N <- 20
k <- 0:N
pk <- expand.grid(p,k)
pr <- dbinom(pk[,2],N,pk[,1])
```
```{r, rgl=TRUE}
plot3d(pk[,1],pk[,2],pr,xlab="Success rate",ylab="No.successes",zlab="Probability/Likelihood")
```
横軸はSuccess rate(成功率)で0から1までの値を取ります。
No.successes(成功回数)は図の奥行きです。
縦軸は、その確率・尤度を表しています。
0,1,2...,20の21通りなので、左から右に向かう、21本の山を持つ点線が表示されています。
一番手前の点線が左から右に向かって下降しています。
これは20回やって0回成功したときには、成功率は確率0と考えるのが最も尤もらしく、
成功率が高くなるにしたがって尤もらしくないと考える、という意味です。
手前から奥に向かって、20回のうち、0,1,2,...,20回成功したそれぞれの場合に、
どのような成功確率が尤もらしいかを描いてある、と読むことができます。
つまり、この点線の山が20回やって、k回成功した場合の尤度です。
一方、success rateの軸に着目すると、成功率を0.01刻みで101通り計算してプロットしてあります。
ある特定の成功率に着目すると、図の奥に向かって、k=0,1,2,...,20についてそれぞれとびとびに21個の値を追いかけることができます。
これが特定の成功率pにおけるk回成功する生起確率です。
このように確率と尤度は、同じ値を仮説(成功率)の軸から眺めるか、結果(成功回数)の軸から眺めるかの違いです。
最尤推定値は尤度関数・対数尤度関数が最大となるようなパラメタの値です。
関数が最大となるとき、関数は極大値を取っていることを利用します。
極大値では関数の一次導関数が0になっていることを利用します。
\(
LL(\mathbf{p}|\mathbf{x}) = \ln{\begin{pmatrix}+\sum_{i=1}^6 x_i\\x_1,x_2,...,x_6\end{pmatrix}} \sum_{i=1}^6 x_i \ln{p_i}
\)
を微分します。
$p_6 = 1- \sum_{i=1}^5 p_i$
なので
$\frac{d p_6}{d p_i} = -1; i=1,2,...,5$
に気をつけて
\(
\frac{\partial LL(\mathbf{p}|\mathbf{x})}{\partial d p_i} = \ln{\begin{pmatrix}\sum_{i=1}^6 x_i\\x_1,x_2,...,x_6\end{pmatrix}} + \sum_{i=1}^6 x_i \ln{p_i}
\)
式を整理して
\(
\frac{\partial LL(\mathbf{p}|\mathbf{x})}{\partial d p_i} = x_i\frac{1}{p_i}-x_6 \frac{1}{p_6}=0
\)
従って
\(
\frac{x_1}{p_1}=\frac{x_2}{p_2}=...=\frac{x_6}{p_6}
\)
であるから、
\(
p_i = \frac{x_i}{\sum_{i=1}^6 x_i}
\)
理想的なサイコロを振って６つの目の出た数を観察する。
その上で出た目に比例した確率(最尤推定値)を対立仮説の生起確率ベクトルとして
尤度比検定を行う。
この試行を何度もくりかえすと
尤度比検定のp値は一様分布になることを示す。
一様分布に従う値をソートしてプロットすると対角線上に並びます。
```{r fig.width=7, fig.height=6}
# 理想的サイコロの目の出る確率
p1 <- rep(1/6,6)
# n.iter回、試行試行する
n.iter <- 1000
# １試行あたりサイコロをN回振る
N <- 100
# 各試行のp値を格納するベクトル
ps <- rep(0,n.iter)
# 試行の繰り返し
for(i in 1:n.iter){
# N回のサイコロの目
s <- sample(1:6,N,replace=TRUE,prob=p1)
# カウントする
s.table <- tabulate(s)
# 対立仮説の最尤推定値
p2 <- s.table/N
# 帰無・対立仮説の対数尤度
LL1 <- lgamma(sum(s.table)+1)-sum(lgamma(s.table+1)) + sum(s.table*log(p1))
LL2 <- lgamma(sum(s.table)+1)-sum(lgamma(s.table+1)) + sum(s.table[which(s.table!=0)]*log(p2[which(s.table!=0)]))
# 対数尤度比の２倍
S <- 2 * (LL2-LL1)
# 自由度は５
df <- 5
ps[i] <- pchisq(S,df,lower.tail=FALSE)
}
# p値をソートしてプロット
plot(sort(ps))
```
- ある集団のジェノタイプ頻度情報が与えられている
- その集団において、ジェノタイプ情報つきの家系がどのくらいの頻度で存在しているかを計算する
- 独立した座位は座位ごとに計算して、積をとる
- 家系内には「集団から直接生まれた子」がいる
- 家系内の「集団から直接生まれた子」とは、家系図で、親がない個人のこと。家系が「単名」の場合には、両親がいて、その両親が二人とも「集団から直接生まれた子」である場合に相当する
- 「集団から直接生まれた子」は、集団に存在しうるジェノタイプのすべてを取りうる
- そのような「集団から直接生まれた子」からスタートした家系の構成員のジェノタイプは、「集団から直接生まれた子」からアレルの伝達によって決まる
- あるジェノタイプ情報つきの家系が「集団の子」である確率は、「集団から直接生まれた子」から、アレルが伝達されて、観察しているようなジェノタイプを持つ構成員となる確率である
- 地道なアプローチ
- すべての個人は、すべてのジェノタイプを取る可能性があると考え、そのすべての場合について、「起きえるか」「起きえないか」、「起きるとしたら、それはどのくらいの確率で起きるか」を計算する
- 「起きるとしたら、どのくらいの確率」というのは、アレルの伝達が、50:50であることを利用して計算する
- 50:50でないのは「集団から直接生まれた子」のジェノタイプ頻度であって、それのみである
- そのようにして、すべての場合を確認した上で、「観察ジェノタイプ」に合致した場合を足し合わせればよい
- 地道なので、わかりやすい
- 大変だけれど、数え上げもいつかは終わる
- 現実的な時間内には終わらない
- 地道な方法は計算が終わらない
- 地道に計算しないで、同じ値を算出したい
- やり方１
- 地道なやり方では、すべての個人がすべてのジェノタイプを取りうるというところからスタートした
- 家系情報とジェノタイプが与えられると、ジェノタイプが確定している人は、そのジェノタイプしかとりえない
- ジェノタイプが与えられていない個人については、家系図上の位置により、メンデルの法則を満足するジェノタイプしかとることができない
- メンデルの法則を満足することを条件に、「取りうるジェノタイプ」を限定する
- 限定されたジェノタイプの場合分けについて、網羅的に計算する
- 計算するにあたって、核家族ごとに計算して、核家族をまたがった確率は、２つの（複数の）核家族に所属する個人のジェノタイプの場合分けにのみ留意して、確率を掛け合わせる
- この方法が\href{http://d.hatena.ne.jp/ryamada22/20110509}{こちら}の記事である
- 核家族の連結をするというアルゴリズムなので、ループのある家系図に対応できない（ループがあるときは、そのループを含む範囲を一塊で取扱い、そのうえで、連結しながら計算することは可能なはず（ただし、面倒くさそうなので、実装していない）
- やり方２
- 別のやり方で計算を省略する
- やり方１は２倍体ジェノタイプに関してメンデルの法則を適用して、場合の数を減らしたし、核家族の連結という考え方も、２倍体を基本とした取扱いである
- こちらのやり方２は、染色体に着目する
- 染色体に着目することで、要素数は２倍となる
- 要素数が増えれば、場合の数を考慮するべき対象も２倍となる
- その代わり、個々の要素の取りうる場合は減少するし、
- 個々の要素の間の関係が単純となるので、処理が単純になるメリットがある
- また、２倍体の場合分けでは、「AND、OR」といった処理が多いため、「数式」での確率計算に適さないのに対して、１倍体での場合分けは、その点もメリットがある
- さて、その具体的なやり方は次節で。
- 概要
- ハプロタイプグラフ
- 家系図(２倍体の伝達関係）から１倍体の伝達の様子（ハプロタイプグラフ）を作る
- 家系図はすべての個人がつながっているが、ハプロタイプグラフはつながっていない
- 母方と父方のグラフに分かれる
- 「子はかすがい」
- 「子が産まれると母方・父方のハプロタイプグラフが子を介して連結する」
- 場合分け
- 伝達パターンで場合分け
- 親が子にアレルを伝えるときに、そのアレルが、その親の母方のものか父方のものかの場合わけ
- こうすると、「枝分かれ」の無い、木グラフになっている
- 母方・父方アレルに振り分ける場合分け
- ２倍体のジェノタイプ（アレル２本分の情報）を母方アレル・父方アレルに振り分ける
- メンデルチェック
- アレルの伝達パターンとして確定したものが得られており
- 父方・母方の伝達アレルも確定して考えているから
- そのような伝達パターン・父母パターンが起こり得るかどうかの確認はとれる
- 木グラフ上のすべてのアレルが同一であれば、「起こり得る」ことで、それ以外は「起こり得ない」こと
- 木の観測確率
- すべての木は、一つの「祖先染色体」を持つ
- 「祖先染色体」は、「集団からの直接の子」の染色体である
- 「祖先染色体」があるアレルを持つ確率は、（HWEを仮定すれば）集団のアレル頻度に等しいので、それが確率
- この集団のアレル頻度と、「場合分けされた木の確率（これは、伝達分岐・父母パターンの割り付けが作る確率なので、$\frac{1}{2^k}$の形をしている
- したがって、木の観測確率は、アレル頻度と$\frac{1}{2^k}$とでできている
- ハプロタイプグラフ全体の観測確率
- ハプロタイプグラフは、複数の連結グラフでできているから、個々の連結グラフを場合分けをにより、木グラフにしたうえで、確率を計算し、すべての連結グラフ（木）での確率の積をとり
- それをすべての場合について足し合わせる
- ソースは[こちら](http://d.hatena.ne.jp/ryamada22/20110518/1305621596 こちら)
```{r}
# packages
# 依存パッケージ
chooseCRANmirror(ind=45)
packages.need <- c("kinship","MCMCpack","gtools","sets","paramlink","rSymPy")
install.packages(packages.need)
library(kinship)
library(MCMCpack)
library(gtools)
library(sets)
library(paramlink)
library(rSymPy)
MakeHaplotypeGraph2<-function(p){
ns<-length(p[,1])
ret<-matrix(0,ns*2,2)
for(i in 1:ns){
tmpmom<-p[i,2]
tmpdad<-p[i,3]
ret[2*i-1,1]<-2*tmpmom-1
ret[2*i-1,2]<-2*tmpmom
ret[2*i,1]<-2*tmpdad-1
ret[2*i,2]<-2*tmpdad
}
ret[which(ret<0)]<-0
ret
}
SeparateGraphs<-function(hG){
ns<-length(hG[,1])
assigned<-rep(0,ns)
cnt<-1
for(i in 1:length(hG[,1])){
if(hG[i,1]!=0){
tmpself<-assigned[i]
tmpp1<-assigned[hG[i,1]]
tmpp2<-assigned[hG[i,2]]
M<-max(tmpself,tmpp1,tmpp2)
if(M==0){
M<-cnt
cnt<-cnt+1
assigned[c(i,hG[i,1],hG[i,2])]<-M
}else{
trioID<-c(i,hG[i,1],hG[i,2])
trio<-c(tmpself,tmpp1,tmpp2)
for(j in 1:length(trio)){
if(trio[j]!=0){
assigned[which(assigned==trio[j])]<-M
}else{
assigned[trioID[j]]<-M
}
}
}
}
}
G<-list()
cnt<-1
hG2<-cbind(1:ns,hG)
for(i in 1:max(assigned)){
if(length(which(assigned==i))!=0){
G[[cnt]]<-hG2[assigned==i,]
cnt<-cnt+1
}
}
G
}
MakeDecsendPattern<-function(hG){
Alt<-which(hG[,1]!=0)
nAlt<-length(Alt)
s<-expand.grid(rep(list(c(1,2)),nAlt))
for(i in 1:(2^nAlt)){
tmp<-matrix(0,nAlt,2)
for(j in 1:nAlt){
tmp[j,1]<-Alt[j]
tmp[j,2]<-hG[Alt[j],s[i,j]]
}
}
Ns<-length(hG[,1])/2
Ls<-length(s[,1])
s2<-matrix(0,Ls,2*Ns)
cnt<-1
for(i in 1:Ns){
if(p[i,2]==0){
s2[,i*2-1]<-1
s2[,i*2]<-2
}else{
s2[,i*2-1]<-s[,cnt]
cnt<-cnt+1
s2[,i*2]<-s[,cnt]
cnt<-cnt+1
}
}
s2
}
MakePatMatPattern<-function(g){
GenotypePlus<-which(g[,1]!=0)
heteros<-rep(2,length(GenotypePlus))
heteros[which(g[GenotypePlus,1]==g[GenotypePlus,2])]<-1
heterolist<-list()
for(i in 1:length(heteros)){
heterolist[[i]]<-1:heteros[i]
}
s3<-expand.grid(heterolist)
s3
}
AssignHaplotype<-function(g,s){
Ns<-length(g[,1])
hapg<-rep(0,Ns*2)
GenotypePlus<-which(g[,1]!=0)
for(j in 1:length(GenotypePlus)){
hapg[GenotypePlus[j]*2-1]<-g[GenotypePlus[j],unlist(s[j])]
hapg[GenotypePlus[j]*2]<-g[GenotypePlus[j],(unlist(s[j])-1.5)*(-1)+1.5]
}
hapg
}
SelectTrees<-function(sepG,v){
ret<-list()
for(i in 1:length(sepG)){
tmp<-sepG[[i]][,1:2]
for(j in 1:length(tmp[,1])){
tmp[j,2]<-sepG[[i]][j,v[tmp[j,1]]+1]
}
tmp2<-rep(0,max(tmp))
cnt<-1
for(j in 1:length(tmp[,1])){
if(tmp[j,1]*tmp[j,2]!=0){
M<-max(tmp2[tmp[j,1]],tmp2[tmp[j,2]])
if(M==0){
M<-cnt
cnt<-cnt+1
tmp2[tmp[j,1]]<-tmp2[tmp[j,2]]<-M
}else{
tmp2[tmp[j,1]]<-tmp2[tmp[j,2]]<-M
}
}
}
ret[[i]]<-tmp2[tmp[,1]]
}
ret
}
MendelCheck<-function(sTrees,sepGraphs,hapg){
ret<-TRUE
ret2<-c()
ret3<-c()
for(i in 1:length(sTrees)){
for(j in 1:max(sTrees[[i]])){
nodes<-sepGraphs[[i]][,1][which(sTrees[[i]]==j)]
numanc<-length(which(sepGraphs[[i]][,2][which(sTrees[[i]]==j)]==0))
tmp<-hapg[nodes][which(hapg[nodes]!="0")]
if(length(as.set(tmp))!=1){
if(length(as.set(tmp))>1){
ret<-FALSE
}
ret2<-c(ret2,"0")
ret3<-c(ret3,1)
}else{
ret2<-c(ret2,hapg[nodes][which(hapg[nodes]!="0")][1])
ret3<-c(ret3,numanc)
}
}
}
return(list(mendel=ret,alleles=ret2,nAncestor=ret3))
}
ProbPerMarker<-function(p,g,As,Ps,hG,sepGraphs,s2,Expression=TRUE){
cnt<-0
s3<-MakePatMatPattern(g)
if(Expression){
Vars<-list()
for(i in 1:length(As)){
Vars[[i]]<-Var(letters[i])
}
retExpression<-Var("0")
}
retProb<-0
for(i in 1:length(s3[,1])){
hapg<-AssignHaplotype(g,s3[i,])
for(j in 1:length(s2[,1])){
sTrees<-SelectTrees(sepGraphs,s2[j,])
MendelOK<-MendelCheck(sTrees,sepGraphs,hapg)
if(MendelOK$mendel){
cnt<-cnt+1
tmpexp<-Var(1/length(s2[,1]))
tmpProb<-rep(0,length(MendelOK$alleles))
for(k in 1:length(tmpProb)){
if(!MendelOK$alleles[k]=="0"){
aid<-which(As==MendelOK$alleles[k])
if(Expression){
tmptmpexp<-Vars[[aid]]
if(MendelOK$nAncestor[k]>1){
for(l in 2:MendelOK$nAncestor[k]){
tmptmpexp<-tmptmpexp*Vars[[aid]]
}
}
tmpexp<-tmpexp*tmptmpexp
}
tmpProb[k]<-Ps[which(As==MendelOK$alleles[k])]
}else{
tmpProb[k]<-1
}
}
if(Expression){
retExpression<-retExpression+tmpexp
}
retProb<-retProb+prod(tmpProb^MendelOK$nAncestor)/length(s2[,1])
}
if(!MendelOK$mendel){
}
}
}
if(!Expression) retExpression=Var("0")
list(prob=retProb,express=retExpression)
}
CalcLikeRatioPerMarker<-function(p,g1,g2,searched,As,Ps,hG,sepGraphs,s2,Expression){
retProb1<-ProbPerMarker(p,g1,As,Ps,hG,sepGraphs,s2,Expression)
retProb2<-ProbPerMarker(p,g2,As,Ps,hG,sepGraphs,s2,Expression)
expres1<-retProb1$express
expres2<-retProb2$express
ProbGenPop<-1
genExp<-Var(1)
for(i in 1:length(searched)){
tmpallele1<-g2[searched[i],1]
tmpallele2<-g2[searched[i],2]
P1<-Ps[which(As==tmpallele1)]
P2<-Ps[which(As==tmpallele2)]
tmp<-P1*P2
tmpVar1<-Var(letters[which(As==tmpallele1)])
tmpVar2<-Var(letters[which(As==tmpallele2)])
tmpgenExp<-tmpVar1*tmpVar2
if(P1!=P2){
tmp<-2*tmp
tmpgenExp<-2*tmpgenExp
}
ProbGenPop<-ProbGenPop*tmp
genExp<-genExp*tmpgenExp
}
retProb1x<-retProb1$prob*ProbGenPop
list(P1=retProb1x,P2=retProb2$prob,Pratio=retProb2$prob/retProb1x,expression=expres2/(expres1*genExp))
}
CalcLikeRatioAllMarker<-function(p,gs,searched,poolid,Gpool,Alleles,Probs,Express=TRUE){
CLRAM<-list()
ped<-MakePedigreeFromFamilyInfo(p)
plot(ped)
hG<-MakeHaplotypeGraph2(p)
sepGraphs<-SeparateGraphs(hG)
s2<-MakeDecsendPattern(hG)
for(na in 1:length(Alleles)){
As<-Alleles[[na]]
Ps<-Probs[[na]]
g<-gs[,,na]
g2<-g
g2[searched,]<-Gpool[poolid,,na]
CLRAM[[na]]<-CalcLikeRatioPerMarker(p,g,g2,searched,As,Ps,hG,sepGraphs,s2,Express)
}
CLRAM
}
PrintStringCLRPM<-function(CLR){
LocusName<-c("D8S1179","D21S11","D7S820","CSF1PO","D3S1358","TH01",
"D13S317","D16S539","D2S1338","D19S433","VWA","TPOX","D18S51","D5S818","FGA")
st<-paste("Locus","LR","Expression",sep="\t")
st2<-paste("Locus","LR","Expression",sep=" ")
print(st2)
cprod<-1
for(i in 1:length(CLR)){
tmp<-paste(LocusName[i],CLR[[i]]$Pratio,sympy(CLR[[i]]$expression),sep="\t")
tmp2<-paste(LocusName[i],CLR[[i]]$Pratio,sympy(CLR[[i]]$expression),sep=" ")
st<-paste(st,tmp,sep="\n")
print(tmp2)
cprod<-cprod*CLR[[i]]$Pratio
}
st2<-paste("cumulative LR",cprod,sep="\t")
st<-paste(st,st2,sep="\n")
print(st2)
st
}
```
```{r}
# 実行コマンド
## 式情報とかをため込むので、たくさんの家系で回すと重い
## 適宜、回し方は調整を要す
searched<-3
st<-""
for(fi in 1:6){
p<-pedigrees[[fi]]
gs<-genotypesFamily[[fi]]
CLout<-CalcLikeRatioAllMarker(p,gs,searched,fi,Gpool[,,],Alleles,Probs,Express=TRUE)
st<-paste(st,PrintStringCLRPM(CLout),sep="\n")
}
st
```
行に行方不明者、列にご遺体をとり、
行列Mの第i行j列には、第i行方不明者が第j遺体である確率が入っているとする。
最適割り付けは、以下のようにすぐ出せます。
ハンガリアン法と呼ばれるアルゴリズムです。
[こちら](http://www.tcp-ip.or.jp/~n01/math/cost.pdf こちら)を参考に。
```{r}
k<-100
library(clue)
# 適当に行列を作ってやる
M<-matrix(abs(rnorm(k^2)),k,k)
# "max":最大になるように割り付け
out<-as.vector(solve_LSAP(M,maximum=TRUE))
# 結果を出力
print(out)
# ここから、「最大値」を計算してみる
VAL<-sum(diag(M[1:k,out]))
# print(out)と同じ値
print(VAL)
```
京都大学医学部法医学講座の玉木敬二先生、ならびに、同講座の皆さま、特に、真鍋翔さん、川合千尋さんには、多くのことを教えていただき、非常に貴重なディスカッションの機会をいただきました。
ここに感謝いたします。
また、法数学勉強会にご参加いただきました皆様には、退屈だったりわかりにくかったり、誤解したまま話したり、ということが多々ありましたこと、この場を借りてお詫びしますとともに、それにもかかわらず、さまざまなトピックについて一緒に勉強する時間を共有させていただきましたことを感謝いたします。
