この記事は、遥か昔のこちらの記事の続きのようなものです。また何度も何度も恐縮ですが、今回の記事内容も付け焼き刃で書いているので色々間違っている可能性があります。お気付きの方は是非ご指摘くださいm(_ _)m
各方面のエコノメトリシャンの方々と上記記事を書いた際に議論*1したことがあるのですが、その時は基本的に統計モデリングを行う際は以下のような判別表に従ってモデルを使い分けるべきだという話になったのでした。
| 確率分布 | 特徴 |
|---|---|
| ポアソン分布 | データが正の離散値、平均値30ぐらいまで、標本平均＝標本分散 |
| 負の二項分布 | データが正の離散値、平均値30ぐらいまで、標本平均＜標本分散 |
| 二項分布 | データが離散値、ゼロ以上でしかも有限 (0, 1, 2, ... N) |
| 正規分布 | データが連続値もしくは離散値でも平均値が十分大*2 (-∞～∞) |
| 対数正規分布 | 同上、ただし正の値、範囲 (0～∞) |
| ガンマ分布 | データが連続値、範囲 (0～∞) |
ところが、現実にはこの判別表に従うとかえってモデリングの精度が悪くなるケースというのが実データを相手にしていると割と頻繁にあって、意外と判断に迷うことが多いんですね。端的に言うと「○○『率』のように明らかにロジスティック回帰の方が当てはまりが良さそうなものであっても、データの取得状況によっては目的変数が正規分布しているように見えて、実際に線形回帰するとR^2でも交差検証(CV)でも精度が同等もしくは優れている場合はどうするべきか」というような話です。
特に「単純にreasoning / interpretationが目的なので推定パラメータ（偏回帰係数）の大小や符号が分かればOK」という場面でどうかとなると、結構悩ましいところです。その場合は確かに線形回帰してもそれなりの結果になることが経験上見込めるからです。加えて線形回帰の方が当然ながら計算負荷も軽いので、それなりの結果になるなら線形回帰すればええやんという話になるのは想像に難くありません。
モデル精度原理主義的には「そんなのCV精度の高い方を選べばええやん」でおしまいだという気もするんですが、それでもデータの生成過程が明確に見えている状況であえて線形回帰にしてもいいんだっけ？というのが個人的にあるので、とりあえずシミュレーションデータでチェックしてみようと思います。
シミュレーションデータで試してみる
ベタッとRで書いてみます。ちなみに、多分これでロジスティック回帰のシミュレーションデータの作り方は合ってると思うんですが。。。間違っていたら悲しいので、おかしなところを見つけた方是非容赦無くご指摘くださいm(_ _)m
# 説明変数 > set.seed(71) > x1 <- runif(5000, 0, 1) > set.seed(72) > x2 <- runif(5000, 0, 1) > set.seed(73) > x3 <- runif(5000, 0, 1) > set.seed(74) > x4 <- runif(5000, 0, 1) > set.seed(75) > err <- rnorm(5000, 0, 1) # 真の偏回帰係数 > b1 <- 6 > b2 <- 3 > b3 <- -3 > b4 <- -6 # 回帰子 > p <- b1*x1 + b2*x2 + b3*x3 + b4*x4 + err # ロジット変換 > y <- plogis(p) > plot(p, y) # データフレームにする > d <- data.frame(y=y, x1=x1, x2=x2, x3=x3, x4=x4)
シミュレーションデータはこれで生成できました。試しに普通にロジスティック回帰してみます。
> d.glm <- glm(y~., d, family=binomial) 警告メッセージ: eval(family$initialize) で: 二項 glm で整数でない成功数がありました! > summary(d.glm) Call: glm(formula = y ~ ., family = binomial, data = d) Deviance Residuals: Min 1Q Median 3Q Max -1.5287 -0.1925 0.0006 0.1918 1.4651 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.09489 0.13048 -0.727 0.467 x1 5.19189 0.17355 29.916 <2e-16 *** x2 2.63917 0.14396 18.333 <2e-16 *** x3 -2.47161 0.14389 -17.177 <2e-16 *** x4 -5.15022 0.16935 -30.412 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 3283.97 on 4999 degrees of freedom Residual deviance: 522.16 on 4995 degrees of freedom AIC: 2914.6 Number of Fisher Scoring iterations: 6
各パラメータを見てみると普通に生成した時の係数に近い結果になりました。次に、区間[0.4, 0.6]でデータを切り取って、これをそれぞれロジスティック回帰と線形回帰してみます。
# 上記区間だけ抜き出す > d2 <- d[d$y>=0.4&d$y<=0.6,] # ロジスティック回帰 > d2.glm <- glm(y~., d2, family=binomial) 警告メッセージ: eval(family$initialize) で: 二項 glm で整数でない成功数がありました! > summary(d2.glm) Call: glm(formula = y ~ ., family = binomial, data = d2) Deviance Residuals: Min 1Q Median 3Q Max -0.258996 -0.096203 0.002195 0.094381 0.219714 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.03648 0.29660 -0.123 0.902 x1 0.35723 0.58223 0.614 0.540 x2 0.20728 0.38355 0.540 0.589 x3 -0.15676 0.38059 -0.412 0.680 x4 -0.32663 0.58364 -0.560 0.576 (Dispersion parameter for binomial family taken to be 1) Null deviance: 7.7568 on 552 degrees of freedom Residual deviance: 7.3232 on 548 degrees of freedom AIC: 770.57 Number of Fisher Scoring iterations: 3 # 線形回帰 > d2.lm <- lm(y~., d2) > summary(d2.lm) Call: lm(formula = y ~ ., data = d2) Residuals: Min 1Q Median 3Q Max -0.129281 -0.048069 0.001069 0.047073 0.109597 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.49089 0.00855 57.414 < 2e-16 *** x1 0.08924 0.01677 5.320 1.51e-07 *** x2 0.05178 0.01105 4.686 3.52e-06 *** x3 -0.03916 0.01097 -3.570 0.000388 *** x4 -0.08160 0.01682 -4.852 1.59e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.05768 on 548 degrees of freedom Multiple R-squared: 0.05611, Adjusted R-squared: 0.04922 F-statistic: 8.145 on 4 and 548 DF, p-value: 2.199e-06
パラメータの大きさ自体は合っていませんが、符号と比だけ見ると線形回帰でもそれなりに元のパラメータ同士の関係性を表していることが分かります。そこで、ここだけ拡大して元データ・ロジスティック回帰での再当てはめ結果・線形回帰での再当てはめ結果の3つを比べてみます。
# 上記区間のx軸の値を取得する > p2 <- p[as.integer(row.names(d2))] # 重ね合わせてプロットしてみる > matplot(p2, cbind(d2$y, predict(d2.glm, newdata=d2, type='response'), predict(d2.lm, newdata=d2)), pch=c(19,19,1), cex=c(1,1,3), xlab='', ylab='') # 相関係数だけざっくり見る > cor(predict(d2.glm, newdata=d2, type='response'), predict(d2.lm, newdata=d2)) [1] 0.9999998
この区間だけで見ると、嫌な言い方ですが「ロジスティック回帰でも線形回帰でも等しく当てはまりが悪い」ことが分かります。一方で、ロジスティック回帰と線形回帰双方の再当てはめ結果はほぼ一致することが分かります。この領域ではどちらで計算しても同じようなモデルになりそうです。
一方、線形回帰を全領域でやるとこうなります。
> d.lm <- lm(y~., d) > summary(d.lm) Call: lm(formula = y ~ ., data = d) Residuals: Min 1Q Median 3Q Max -0.5482 -0.1174 -0.0018 0.1195 0.5076 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.489555 0.008282 59.11 <2e-16 *** x1 0.727627 0.008112 89.70 <2e-16 *** x2 0.358074 0.008007 44.72 <2e-16 *** x3 -0.333494 0.008069 -41.33 <2e-16 *** x4 -0.731140 0.007930 -92.20 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.164 on 4995 degrees of freedom Multiple R-squared: 0.7991, Adjusted R-squared: 0.7989 F-statistic: 4967 on 4 and 4995 DF, p-value: < 2.2e-16
本来の偏回帰係数とは似ても似つかないパラメータの推定結果になりました。ただし「符号と比」についてはやはり狭い区間の時と同様に真の偏回帰係数の関係が保たれているようです。そこで、これをやはり同様に全区間で再当てはめして比べてみましょう。
# 重ね合わせてプロットしてみる > matplot(p, cbind(d$y, predict(d.glm, newdata=d, type='response'), predict(d.lm, newdata=d)), pch=c(19,19,1), cex=c(2,1,1), xlab='', ylab='') # ざっくり相関係数を見てみる > cor(d$y, predict(d.glm, newdata=d, type='response')) [1] 0.9268167 > cor(d$y, predict(d.lm, newdata=d)) [1] 0.8939221 > cor(predict(d.glm, newdata=d, type='response'), predict(d.lm, newdata=d)) [1] 0.9627183
当たり前ですが、本来の定義域が[0, 1]なのでロジスティック回帰で再当てはめすると綺麗に元データに沿うように収まるものの、線形回帰で再当てはめすると物の見事に上下に飛び出してしまいます。相関係数のようなざっくりとした指標だと分かりづらいですが、プロットで見比べると一目瞭然です。
ということである意味予想された結果ながら、「厳密に偏回帰係数を求めようと思ったらロジスティック回帰でなければいけない」ものの「単純にreasoning / interpretationのためにパラメータの大小や比だけ分かれば良い状況なら線形回帰で代用しても（場合によっては）差し支えない」ということが分かりました。
おそらく最も妥当と思われる説明
以前よくお世話になっていたアンコレこと@uncorrelatedさんから、ズバリこのようなコメントをいただきました。これ考えてみたら当たり前で、GLMであればロジスティック回帰に限らず、テイラー展開した際に一次の項までで十分に近似できてしまうようなサンプルに対しては、事実上線形回帰で表せるわけです*3。そりゃそういうケースでは線形回帰（≒一次の項）であってもR^2見ようがCV誤差見ようが精度はそれなりになるよな、と。。。
感想
そのうち別の記事でだらだら書こうかと思っていますが、実務の現場だと原理主義的にベストのモデルを選択するよりも「とりあえず目の前の課題に適したモデルでクイックに結果を出すべき」というシチュエーションが少なくありません。今回は「GLMであるべきものを線形回帰で代用しても問題ないかどうか」というポイントに絞ってみましたが、その他のケースでも同じような結論に至るものが実は結構あるのかな？と思った次第です。
追記
上記のシミュレーションですが、明らかにロジスティック分布の「ノイズなし」バージョンになっていてこれだと幾ら何でもひどくないか的な雰囲気があり。。。どなたかベストのシミュレーション生成方法をご存知でしたら是非ご教示くださいm(_ _)m