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\')
参照