カテゴリカル変数に対するロジスティック回帰分析
カテゴリカル変数を目的変数にとったロジスティック回帰分析についてまとめてみます．
主に『カテゴリカルデータ解析 (Rで学ぶデータサイエンス 1)』 藤井良宜 2010年 共立出版 の第6章を参考にしました．
ロジスティック回帰分析とは
カテゴリカル変数として表された目的変数を，その他の変数で説明するモデルを立てて分析する手法のひとつです． ロジスティック回帰分析では，説明変数 x と目的変数 y の関係が直線的ではなく，シグモイド型に変化していくモデル (たとえば，x が大きくなれば y も大きくなるが，ある程度のところを超えると増加の度合いが鈍っていくというようなかたち) を仮定します．目的変数には，2値のカテゴリー，3つ以上のカテゴリー，順序のあるカテゴリーをとることができます．目的変数が2値の場合
まず，目的変数が2値で，説明変数が連続変数の場合を説明します．書籍にならい，パッケージ vcd にあるスペースシャトルのデータ SpaceShuttle を例にとります．スペースシャトルの打ち上げテスト時に，O-リングがひとつでも故障したかどうかを調べたデータです．
library(vcd) head(SpaceShuttle) FlightNumber Temperature Pressure Fail nFailures Damage 1 1 66 50 no 0 0 2 2 70 50 yes 1 4 3 3 69 50 no 0 0 4 4 80 50 <NA> NA NA 5 5 68 50 no 0 0 6 6 67 50 no 0 0
ここでは，Temperature と Pressure (どちらも連続的な説明変数) によって，Fail (故障ありかなしという2値の目的変数) を説明するモデルを立ててみます．一般的なかたちは
となります (a, b は定数)．二項分布 (または多項分布) にしたがうデータに対して，連結関数にロジット変換を用いた一般化線形モデル (GLM) ※1 をあてはめて解析する手法と言うこともできます．そうすると，[tex:frac p{1-p} は，ふたつのカテゴリーの確率の比と考えることができます．
解析では，GLMの枠組みでモデルを構築しますので，glm 関数を使用します．
shuttle.tp <- glm(Fail ~ Temperature + Pressure, data = SpaceShuttle, family = binomial)
summary(shuttle.tp)
Call:
glm(formula = Fail ~ Temperature + Pressure, family = binomial,
data = SpaceShuttle)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.2130 -0.6089 -0.3870 0.3472 2.0928
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 14.359775 7.442899 1.929 0.0537 .
Temperature -0.241540 0.109722 -2.201 0.0277 *
Pressure 0.009534 0.008703 1.095 0.2733
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 28.267 on 22 degrees of freedom
Residual deviance: 18.972 on 20 degrees of freedom
(1 observation deleted due to missingness)
AIC: 24.972
Number of Fisher Scoring iterations: 5結果を見ると，Pressure は Fail の説明にそれほど効いていないようです．
ですので，Temperature だけで Fail を説明するモデルを立ててみます．
shuttle.t <- glm(Fail ~ Temperature, data = SpaceShuttle, family = binomial)
summary(shuttle.t)
Call:
glm(formula = Fail ~ Temperature, family = binomial, data = SpaceShuttle)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.0611 -0.7613 -0.3783 0.4524 2.2175
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 15.0429 7.3786 2.039 0.0415 *
Temperature -0.2322 0.1082 -2.145 0.0320 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 28.267 on 22 degrees of freedom
Residual deviance: 20.315 on 21 degrees of freedom
(1 observation deleted due to missingness)
AIC: 24.315
Number of Fisher Scoring iterations: 5AICを比較してみると，前者 (24.972) に比べて後者 (24.315) のモデルのほうが若干小さい値となっているので，Temperature のみで Fail を説明するモデルのほうがリーズナブルなようです．
モデルで予測される値を出力してみるとこんな感じ．
library(ggplot2)
fail.t.pred <- predict(shuttle.t, type = "response")
fail.obs <- as.integer(SpaceShuttle$Fail[!is.na(SpaceShuttle$Fail)]) - 1
shuttle.df <- data.frame(
Temperature = rep(SpaceShuttle$Temperature[!is.na(SpaceShuttle$Fail)], 2),
Fail = c(fail.obs, fail.t.pred),
Class = c(rep("Observed", length(fail.obs)), rep("Predicted", length(fail.t.pred))))
gg <- ggplot(shuttle.df, aes(Temperature, Fail))
gg + geom_point(aes(colour = Class))
でも，あまりフィットしてるように見えませんね…
目的変数が3つ以上でカテゴリー間に順序がないとき
ここでも書籍にならって，iris のデータを使います．フィッシャーが計測したアヤメのデータですね．head(iris) Sepal.Length Sepal.Width Petal.Length Petal.Width Species 1 5.1 3.5 1.4 0.2 setosa 2 4.9 3.0 1.4 0.2 setosa 3 4.7 3.2 1.3 0.2 setosa 4 4.6 3.1 1.5 0.2 setosa 5 5.0 3.6 1.4 0.2 setosa 6 5.4 3.9 1.7 0.4 setosa
図示するとこうなります．
library(ggplot2) gg <- ggplot(iris, aes(Petal.Length, Petal.Width)) gg + geom_point(aes(color = Species))
setas, versicolor, virginica の3種を目的変数に，雄しべの計測値を説明変数にとり，アヤメの種を推定するモデルを立ててみます．
今度は，目的変数のカテゴリー数が2より多いので，多項分布モデルを適用して，各カテゴリーの確率を考える必要があります．ひとつのカテゴリーを基準として，それと他のカテゴリーの確率の比を考えます．
以下が，このときのモデルになります．
解析には，パッケージ nnet の multinom 関数を使用します．これは，多項分布モデル専用の glm 関数のような感じです (使い方は glm とほぼ同じだが，多項分布専用なので，family の指定はない)．
install.packages("nnet")
library(net)
iris.multi <- multinom(Species ~ Petal.Length + Petal.Width, data = iris, maxit = 1e4)
# weights: 12 (6 variable)
initial value 164.791843
iter 10 value 12.657828
iter 20 value 10.374056
iter 30 value 10.330881
iter 40 value 10.306926
iter 50 value 10.300057
iter 60 value 10.296452
iter 70 value 10.294046
iter 80 value 10.292029
iter 90 value 10.291154
iter 100 value 10.289505
iter 110 value 10.289258
iter 120 value 10.287417
iter 130 value 10.286160
final value 10.286136
convergedsummary(iris.multi)
Call:
multinom(formula = Species ~ Petal.Length + Petal.Width, data = iris,
maxit = 10000)
Coefficients:
(Intercept) Petal.Length Petal.Width
versicolor -24.48672 7.26604 8.626438
virginica -69.57214 12.99170 19.045286
Std. Errors:
(Intercept) Petal.Length Petal.Width
versicolor 59.26616 50.07007 110.0120
virginica 60.79077 50.12256 110.0756
Residual Deviance: 20.57227
AIC: 32.57227 モデルが複雑なので，内部でニューラルネットワークを利用した計算を動かしているそうですが，勉強不足なのでここでは触れません (汗)．
しかし何にせよ，setas に比べて versicolor と virsinica は Length も Width も大きく (Coefficients が正のため)，それらが大きくなるほど virsinica である確率が高くなる (virsinica の Coefficients のほうが大きいため)，という解釈ができるようです ※2．
目的変数が3つ以上でカテゴリー間に順序があるとき
パッケージ MASS に含まれる polr 関数で比例オッズモデルを利用するようですが，理解しきれていないのでまた今度．注
※1 一般化線形モデルについては，生態学のデータ解析をご参照ください．※2しかしこの場合は，判別分析を採用したほうが良いような気もしますが…
