Rでロジスティック回帰分析の結果を見る方法
実習
内蔵データturnoutデータを読み込んでみよう。
turnoutは1992年のアメリカ総選挙に関するデータであり、race（人種）、age（年齢）、educate（教育水準）、income（収入）によるvote（投票の有無）を把握できる。このデータは投票をしたかしなかったかという従属変数に関心があるので、ロジスティック回帰分析を使用できる。
ロジスティック回帰分析は一般的な回帰分析とは異なり、glm()関数を通じてモデルを立てる。このときオプションとしてfamily=binomial()を入れてやれば、自動的にロジスティック回帰分析をしてくれる。突然「バイノミアル」が出てくる理由は、従属変数が0か1かの二項変数だからである。
結果の解釈
結果を見ると、一般的な多重回帰分析と似ているようで異なる点が多い。各変数に対する回帰係数の見方は同じだが、F検定がなくなり、逸脱度devianceというものが登場した。
逸脱度とは、ロジスティック回帰モデルがどれほどデータを説明できていないかの尺度と見なして差し支えない。零逸脱度null devianceは何の変数もなく定数項だけがあるときの逸脱度であり、データが全くない最悪の状況だと言える。こうして得られた残差逸脱度residual devianceは小さければ小さいほど良く、カイ二乗分布に従うため、カイ二乗適合度検定を通じてモデルが適合しているか確認できる。
glm()関数の名前は「一般化された回帰分析」を意味してはいるが、基本的に回帰分析であるため、多重共線性を考慮しないわけにはいかない。幸いにも、VIFはvif()関数を使用することで簡単に求めることができる。例では多重共線性の問題はないように見える。
モデル診断
一方、残差プロットが非常に衝撃的だが、幸いにもロジスティック回帰分析において残差プロットは何の意味も持たない。したがって、わざわざ確認する必要もなく、確認してみて異常でも気にする必要はない。先に述べた適合度検定があれば十分である。
さて、実際にいくつかの方法を通じてモデルが適合しているか確認してみよう。
カイ二乗適合度検定を行う方法は非常に簡単だ。有意水準$\alpha$に対して検定したければ、qchisq()関数に$(1 - \alpha)$と残差逸脱度の自由度を入れて臨界値を計算すればよい。上のように臨界値より残差逸脱度が小さければ、モデルは適合していると見なす。
ホスマー・レメショウ適合度検定は、ロジスティック回帰分析で使用される代表的な適合度検定であり、同様にカイ二乗統計量を通じてモデルが適合しているか検定してくれる。ResourceSelectionパッケージのhoslem.test()関数にモデルの実際の従属変数と適合値を入れることで検定できる。帰無仮説は「モデルが適合している」であるため、上のように有意確率が高く帰無仮説を棄却できなければ、モデルは適合していると見なす。しかし、フランク・ハレルによれば、いくつかの弱点があり、今ではもはや推奨されないという1。
rmsパッケージのlrm()関数を使用すると、モデルの尤度比検定をしてくれる。上の二つの方法が適合度検定であり帰無仮説が「モデルが適合している」であるのとは異なり、尤度比検定では帰無仮説を棄却してこそモデルが適合していると見なす。どうせrmsパッケージを使っているなら、わざわざglm()関数でロジスティック回帰モデルを作ってそのアウトプットを入れるのではなく、そのままlrm()関数にモデルを直接入れても構わない。例で言えば、lrm(vote~.,data=turnout)と入力することで、すべての過程を一度に実行してくれる。
予測
最後に、ロジスティック回帰モデルを利用して直接確率を計算してみよう。確率はpredict()関数にtype='response'オプションを入れることで簡単に計算される。
もちろん、我々はすでにデータを持っていて真の値を知っているため、このような予測は大きな意味を持たない。したがって、newdataオプションを通じて新しいデータを与え、そのときもうまく当てられるかを確認する過程が必要である。
コード
以下は例のコードだ。
install.packages("Zelig")
install.packages("car")
install.packages("ResourceSelection")
install.packages('rms')
library(Zelig)
data(turnout); head(turnout)
out0<-glm(vote~.,family=binomial(),data=turnout); summary(out0)
library(car)
vif(out0)
win.graph(4,4); plot(out0$residuals,main="잔차그림")
qchisq(0.95,df=1995)
library(ResourceSelection)
hoslem.test(out0$y,fitted(out0))
library(rms)
lrm(out0)
lrm(vote~.,data=turnout)
predict(out0,type='response\')