独立したN人の患者のアウトカム\(y\)を考える。
i番目の患者において得られた観測値を\(y_i(i=1–N)\)とすると、
\(y_i\)は確率密度関数\(p(y_i;Θ_i)\)で表される確率分布に従う。
仮に確率分布が正規分布なら、\(y_i\sim N(μ_i,\sigma^2_i)\)ということになるが、ここではあえて確率分布を明示しない。
確率密度関数\(p(y_i;Θ_i)\)において、Θは確率密度関数を規定するパラメータである。
\(y_i\)が全て観測されているとき、このデータと最も当てはまりが良い確率密度関数のパラメータΘを推定したい。
尤度原理
データから計算された尤度関数が高くなるようなパラメータの値がよい推定値である
に基づいて、パラメータΘを推定する方法を最尤推定法(MLE)と呼ぶ。
確率密度関数\(p(y_i;Θ_i)\)は、患者iにおける確率密度関数であり、N人の観測値はすべて独立していることから、N人のアウトカムがそれぞれ観測値\(y_1,y_2,...,y_i\)になる同時確率分布は
\[ p(\vec y;\vecΘ )=p(y_1;Θ_1)*p(y_2;Θ_2)*...*p(y_N;Θ_N)=\prod_{i=1}^Np(y_i;Θ_i) \]
と表せる。
尤度関数は\(y_i\)のもとで\(Θ_i\)となる尤もらしさを表す関数であるから、\(p(Θ_i;y_i)\)と表すことができる。
尤度関数 \(L(\vecΘ)=p(\vecΘ;\vec y)=\prod_{i=1}^Np(Θ_i;y_i)\) とすると
対数尤度関数は
\[ l(\vecΘ)=\log[L(\vecΘ)]=\log(\prod_{i=1}^Np(Θ_i;y_i))\\ =\log[p(Θ_1;y_1)*p(Θ_2;y_2)*...*p(Θ_N;y_N)]\\ =\log p(Θ_1;y_1)+\log p(Θ_2;y_2)+...\log p(Θ_N;y_N)\\ =\sum_{i=1}^N\log[p(Θ_i;y_i)] \]
と表せる。
最尤推定量(MLE)を推定するときは、総積\(\prod\)を総和\(\sum\)に変換できて計算上便利なので\(L(\vecΘ)\)の代わりに\(l(\vecΘ)\)を用いる。
尤度原理より、MLEは尤度関数(対数尤度関数)の最大値と一致するので、\(l(\vecΘ)\)の最大値を求めれば良い。
多くの確率分布において、\(l(\vecΘ)\)は上に凸の関数形をしているので、\(l(\vecΘ)\)の最大値は、\(l(\vecΘ)\)の頂点を求めることと同じである。つまり、\(l(\vecΘ)\)の1階微分=0の解を求めれば良い。
すなわち、\({l(\vecΘ)}'=\frac{\partial{l(Θ)}}{\partial{Θ}}=0\)の解\(\hatΘ\) がMLEとなる。
ここで、\({l(\vecΘ)}'=\frac{\partial{l(Θ)}}{\partial{Θ}}=0\)をスコア方程式、\(U(Θ)={l(\vecΘ)}'=\frac{\partial{l(Θ)}}{\partial{Θ}}\)をスコア関数と呼ぶ
対数尤度関数は多くの確率分布において上に凸の関数形をしているので、二次関数で近似することができる。\(\hatΘ\)をMLE、SEを標準誤差とすると、
\[ l(Θ)\fallingdotseq\frac{-1}{2}(\frac{\hatΘ-Θ}{SE})^2 \]
と近似することができる。
Wald信頼区間は、対数尤度関数を二次関数で近似できることを利用した信頼区間の構成方法である。
\[ \frac{-1}{2}(\frac{\hatΘ-Θ}{SE})^2=-1.92 \]
という方程式をΘについて解くことで、以下のように信頼区間を求めることができる。
\[ 95\%\ Wald\ CI=\hatΘ±1.96\ SE \]
MLEは、対数尤度関数\(l(Θ)\)から導かれるスコア関数\(U(Θ)={l(\vecΘ)}'=\frac{\partial{l(Θ)}}{\partial{Θ}}\)の最大値で推定できた。
MLEを良く推定できるような対数尤度関数とはどんな形になるのか考えてみる。
対数尤度関数は上に凸の形をしており、そのピークの点がMLEになる。そして、そこから-1.92だけ下がった二点が95%Wald信頼区間であった。
ということは、上に凸のピークが鋭く尖っているほど、MLEとその信頼区間は狭く確度高く推定できるというふうに考えられる。
ここで、簡単のために対数尤度関数の近似式を以下(A),(B)のように考えてみる
\[ l(Θ)=-Θ^2 ...(A)\\ l(Θ)=-5Θ^2...(B) \]
ピークが鋭くなっているのは(B)の方であり、MLEに対する分散(95%CI)は(A)に比べて小さくて推定精度が高いことがわかる。
このピークの鋭さは、 “傾きの変化率”が大きい(急な)ほうが鋭く、“傾きの変化率が”小さい(緩やか)なほうが鈍いピークになる。これを定量的に表すと、対数尤度関数の傾きの傾き、つまり対数尤度関数の二階微分になる。
すなわち、
$$ (A)について一階微分するとl’(Θ)=-2Θ\二階微分すると l’‘(Θ)=-2\ (B)について一階微分するとl’(Θ)=-10Θ\ 二階微分するとl’’(Θ)=-10
$$
となり、(B)のほうが傾きの変化率が負に大きいことがわかる。
定式化すると
\[ I(Θ)=E[-{l''(\vecΘ)}]\\ =E[-\frac{\partial^2{l(Θ)}}{\partial^2{Θ}}] \]
となる。ここで、傾きの変化率は負の値なので、負の符号をつけて正の値にすることで大きさの比較をしやすくし、さらに期待値を取ることで、関数全体でのピークの鋭さを表している(と理解したが、間違っていたら教えてください)。
\(**I(Θ)\)をFisher情報量(行列)**と呼び、この値が大きいほど、対数尤度関数のピークが鋭く、MLEの推定精度が高いことを意味する。
つまりFisher情報量が多いほうが、推定精度が高い。情報量が少ないほうが、推定精度が低いことになる。これは直感的にも理解できる。
スコア関数\(\vec U(Θ)={l(Θ)}'=\frac{\partial{l(Θ)}}{\partial{Θ}}\) とするとき、
Fisher情報量行列は\(\vec I(Θ)=E[-\frac{\partial^2{l(Θ)}}{\partial^2{Θ}}]=E[\vec U(Θ)\vec U^T(Θ)]=Var[U(Θ)]\)と表せる
行列ではなく、Fisher情報量の場合、\(I(Θ)=E[-l''(Θ)]=E[U^2(Θ)]=Var[U(Θ)]\) と表せる
対数尤度関数のスコア関数を\(U(Θ)=I'(Θ)=\frac{\partial{l(Θ)}}{\partial{Θ}}\)
Fisher情報量を\(I(Θ)=E[-l''(Θ)]=E[U^2(Θ)]=Var[U(Θ)]\)とする
任意のパラメータΘの不偏推定量の分散とFisher情報量の関係性は
\[ Var(\hatΘ)≧ \frac{1}{I(Θ)} \]
という不等式が成り立つ。
Fisher情報量が大きいほど、分散の下限は小さく、Fisher情報量が小さいほど、分散の下限は大きくなる関係性を示している。
\[ MLE(\hatΘ)\xrightarrow[n\to\infty]{\mathrm{d}}N(Θ,\vec I(Θ)^{-1}) \]
MLEのサンプルサイズが大きいとき、その推定値の従う確率分布は平均Θ、分散\(\vec I(Θ)^{-1}\)(つまりFisher情報量行列の逆数)の多変量正規分布に従う。
ここで、MLEの従う分布の分散\(\vec I(Θ)^{-1}\)は、Cramer–Raoの下限そのものである。
MLEの重要な性質として