2018年の統計応用・医薬生物学にロジスティック回帰の式とモデル選択について問題が出ていたので基本的な概観を書いてみます。
基本さえ押さえていれば計算が煩雑でないので、知っていれば結構簡単な問題だったと思うのですが、逆に数理的な背景を知らないとさっぱりです。統計応用はそのパターンが多いですね・・・。
目次：
- 一般線形モデル（general linear model）
- 一般化線形モデル(generalized linear model)とリンク関数
- ロジスティック回帰分析と調整オッズ比
- 赤池情報規準（AIC）
- カルバックライブラー情報量
一般線形モデル（general linear model）
ロジスティック回帰分析は一般化線形モデルと呼ばれるものの一種です。まずそもそも一般化線形モデルってなんやねんと思うわけですが、その前に一般線形モデル(general linear model)を見てみましょう。
知りたい数値である結果変数yとそれに影響を及ぼす説明変数x()として
\(y=\beta_0+\beta_1x_1+\beta_2x_2+…+\beta_px_p+\epsilon\)
のような関係で見たモデルを一般線形モデルと言います。βはパラメータと呼ばれ、それぞれのxの影響の大きさに関連します。\(\epsilon\)は誤差項ですね。また右辺をまとめて線形予測子と呼ばれています。
この式に基づいて分析していくのは、いわゆる重回帰分析と呼ばれる方法になります。最小二乗法と呼ばれる方法でβを求めていきます。この辺は以前にも記事を書きました。
一般化線形モデル(generalized linear model)とリンク関数
これに対してロジスティック回帰分析は一般化線形モデルと呼ばれるモデルの一種です。ロジスティック回分析では結果変数と説明変数が次のような関係になります。
\(\log\frac{y}{1-y}=\beta_0+\beta_1x_1+\beta_2x_2+…+\beta_px_p\)
左辺がロジット関数と呼ばれる関数の形になっています。左辺の関数の形をリンク関数と呼びます。これが先ほどの一般線形モデルのように、そのままyではないものが一般化線形モデルです。
このような結果変数・説明変数の式関係をもとに分析をするのがロジスティック回帰分析です。
結果変数・説明変数の関係性を線形予測子・リンク関数を用いてモデリングしますが、そのやり方によってモデルが変わるわけですね。
ロジスティック回帰分析と調整オッズ比
先ほどの線形予測子とリンク関数の式を色々と変形してみると、とても便利な形になっていることがわかります。まず結果変数yを左辺に持ってくるようにしますと
\(y=\frac{\exp(\beta_0+\beta_1x_1+…+\beta_px_p)}{1+\exp(\beta_0+\beta_1x_1+…+\beta_px_p)}\)
となります。この式から0<y<1であることがわかりますので、結果変数は確率を示すのに適していることがわかります。
よくあるのは2値データの起こる確率をロジスティック回帰で予測するというものですね。
結果変数をθとして治療で何らかの反応を示す確率としてみます。つまり、先ほどのyをθで置き換えて
\(\theta=\frac{\exp(\beta_0+\beta_1x_1+…+\beta_px_p)}{1+\exp(\beta_0+\beta_1x_1+…+\beta_px_p)}\)
とします。
得られたデータに関して反応があったときにy=1、ないときにy=0とすると
\(P(Y=1)=\theta\)
\(P(Y=0)=1-\theta\)
となります。
そうするとYは確率θのベルヌーイ分布に従う確率変数であることがわかります。
このとき、パラメータβの値が各説明変数のオッズ比の対数を取ったもの(対数オッズ比)を示しています。過去問でも出されていたので、それを確認してみます。
まず上述の条件の時、対数オッズは
\(\log\frac{\theta}{1-\theta}=\beta_0+\beta_1x_1+\beta_2x_2+…+\beta_px_p\)
となります。
さて、ここで疾患の反応確率に関わるある因子を持つ患者を\(x_1=1\)、因子を持たない患者を\(x_1=0\)としてみると対数オッズ比は
\(\log\frac{\frac{P(Y=1|X_1=1)}{P(Y=0|X_1=1)}}{\frac{P(Y=1|X_1=0)}{P(Y=0|X_1=0)}}\\=\log\frac{\beta_0+\beta_1+\beta_2x_2+…+\beta_px_p}{\beta_0+\beta_2x_2+…+\beta_px_p}\\=\beta_1\)
となります。
よって\(\beta_1\)が対数オッズ比となることが分かりました。
ここから説明変数\(x_1\)が1増えるとオッズ比がどう変動するかがわかります。この時のオッズ比は他の説明変数を固定したときの変化を示しているため、調整オッズ比(adjusted odds ratio)と言われます。
赤池情報規準（AIC）
これも問題に出ていたので触れておきます。
説明変数を選択する際に「データへの当てはまりの良さ」を優先すると当然ながら説明変数を増やせば増やすほど精度は良くなります。
ただ、説明変数があまりに増えると煩雑であったり、また次に得られたデータを予測する際には役に立たない可能性も十分あります。
そこで、「予測の良さ」に焦点を当てたときに、どの説明変数のモデルであれば良いのかを評価する規準の一つが赤池情報規準です。
以前その意味と導出に関しては記事を書きました。
式としては
AIC=-2log対数尤度＋2×最尤推定を行なったパラメータ数
となります。
先ほどの例で考えますと、対数尤度はYがベルヌーイ分布を取ることを利用して導出できます。得られたデータがn個あるとして、それぞれについて
\(\theta_i=\frac{\exp(\beta_0+\beta_1x_{1i}+…+\beta_px_{pi})}{1+\exp(\beta_0+\beta_1x_{1i}+…+\beta_px_{pi})}\)
が成り立ちます。(i=1,2,3,…,n)
まず尤度関数は
\(L=\Pi_{i=1}^n{\theta_i}^{y_i}(1-\theta_i)^{1-y_i}\)
となります。
ここで
\(\theta^{y_i}=\frac{\exp{y_i}(\beta_0+\beta_1x_{1i}+…+\beta_px_{pi})}{\{1+\exp(\beta_0+\beta_1x_{1i}+…+\beta_px_{pi})\}^{y_i}}\)
であり
\((1-\theta)^{1-y_i}=\frac{1}{\{1+\exp(\beta_0+\beta_1x_{1i}+…+\beta_px_{pi})\}^{1-y_i}}\)
なので、対数尤度は
\(\log{L}=\sum y_i(\beta_0+\beta_1x_{1i}+…+\beta_px_{pi})-\sum \log\{1+\exp(\beta_0+\beta_1x_{1i}+…+\beta_px_{pi})\}\)
となります。
よって、AICはβを\(\hat\beta\)で最尤推定したとき
\(AIC=-2\sum y_i(\hat\beta_0+\hat\beta_1x_{1i}+…+\hat\beta_px_{pi})+2\sum \log\{1+\exp(\hat\beta_0+\hat\beta_1x_{1i}+…+\hat\beta_px_{pi})\}+2(p+1)\)
となります。
カルバックライブラー情報量
さて、これもまた問題に出ていたので触れてみますが、AICはカルバックライブラー情報量が小さくなるようにしてモデルの良さを評価する方法です。
カルバックライブラー情報量とは真の分布があると仮定したとき、モデルとなっている分布と真の分布の違いの大きさを数値化したものです。真の分布との“距離“とは厳密には違うようです（真の分布からモデルとなる分布を見るか、モデルとなる分布から真の分布を見るかで数値が変わるため）。
式としては真の分布をq(x), モデルとなる分布をf(x)としたとき、真の分布による期待値をとって
\(E\left[\log\frac{q(x)}{f(x)}\right]\)
と表されます。
具体例として今回のようなベルヌーイ分布の場合を見てみます。
真の分布を確率θのベルヌーイ分布q(x)として、モデルとなる分布を確率πのベルヌーイ分布f(x)としてみると
\(E[log\frac{q(x)}{f(x)}]=E[logq(x)]-E[logf(x)]\\=\sum_{x=0}^1q(x)logq(x)-\sum_{x=0}^1q(x)logf(x)\\=(1-\theta)log(1-\theta)+\theta log\theta-(1-\theta)log(1-\pi)-\theta log\pi\\=\theta log\frac{\theta}{\pi}+(1-\theta)log\frac{1-\theta}{1-\pi}\)
となります。
これだけわかっていれば過去問も簡単に解けるのですが、範囲が広い分対処する自信が無くなりますね、、、。各分野のあんまり細かいところは出ない気がするのでそれなりに幅広く基本を理解しておく必要があるのかなと思ってます。
参考文献：
いつも参考にしてます。
コメントを残す