離散変数の固定／ランダム効果モデル
太郎丸 博
2015年6月18日
1 Dichotomous Variable: 二値変数が従属変数の場合
二値変数が従属変数の場合でも、連続変数とみなして固定効果モデルやランダム効果モデルでパラメータを推定することは可能である。実際、かなりよい近似が得られるという報告も複数ある。しかし、不適切な場合もあるので、ロジスティック回帰分析やプロビット分析のような手法を使うほうが無難であるし、理論上好ましいといえる。
以下では、まず横断的データのロジスティック回帰分析について簡単に紹介し、それをパネルデータに応用する方法についても概観する。
モデルからの予測値が 0 より小さくなったり、1 より大きくなったりする場合、単純な線形モデルをあてはめるのは不適切であろう。横断的データの分析の場合、二値の従属変数には、ロジスティック回帰分析やプロビット分析が用いられるが、これを固定効果モデルやランダム効果モデルに応用するやり方を論じていこう。
1.1 Problem: 線形回帰の問題
例えば、鉄棒で逆上がりができるかどうかを従属変数、逆上がりの練習を始めて何日目かを独立変数として、回帰分析をしたとしよう。仮に図1 のような結果が出たとする。
x <- rep(1:30, each=10)
xc <- cut(x, c(0, 5, 10, 15, 20, 25, 31))
n <- length(x)
odds <- exp(-8.5 + 0.6*x)
y <- rbinom(n, 1, p= odds/(1 + odds))
m.y <- tapply(y, xc, mean)
par(mfrow=c(1,2), mar=c(3,3,2,0.5), mgp=c(2, 0.7, 0))
plot(y ~ jitter(x), ylim=c(-0.1, 1.2), main="習得率と回帰直線つきの散布図")
lm1 <- lm(y ~ x)
abline(lm1, col="blue", lty=2)
lines(2.5 + (0:5) * 5, m.y, col="red", lwd=2)
legend("right", col=c("blue", "red"), lty=c(2,1), legend=c("回帰直線", "習得率"))
plot(jitter(residuals(lm1)) ~ jitter(x), main="残差のプロット")
m.resid <- tapply(residuals(lm1, "response"), xc, mean)
lines(2.5 + (0:5) * 5, m.resid, col="red", lwd=2)
legend("topright", col=c("red"), lty=1, legend=c("残差の平均"))
図1: 二値変数を線形回帰分析した場合の予測値と残差
回帰直線は、逆上がりができるかどうかの二値変数の平均値を予測しているとも解釈できるので、それは習得率を予測しているモデルでもあるといえる。実際の習得率（５日刻みで平均値を計算）と比較すると、そこそこ近似できてはいるが、３０日に近くなると、予測値が１を超えてしまう。このように場合によっては絶対にありえない値を予測してしまうため、二値変数を通常の回帰分析で予測するのは、あまり好ましくない。しかしかなりよく近似できる場合もあるので、実践的には必ずしもダメというわけではない。
もう一つの問題は、残差の期待値と分散である。残差をプロットしてみると、 x が 25 あたりでは習得率の予測値が 1 で、実際にほとんど全員できるので、残差の分散はほとんどゼロになるが、10～15 日目あたりはできる子とできない子の数がほとんど同じなので分散は、\(0.5^2=0.25\) ぐらいになる。期待値も 0 より大きくなったり小さくなったりを繰り返している。これは OLS の前提を満たしていない。
1.2 Logit: ロジスティック回帰分析
そこで、ロジスティック回帰分析やプロビット分析がよく用いられる。ロジスティック回帰分析とは、モデルから予測される従属変数が 1 をとる確率を \(\hat p_i\) とすると、 \[ \log \frac{\hat p_i}{1- \hat p_i} = \beta_0 + \beta_1 X_i \qquad (1) \] というモデルをたて、最尤法で推定する（サンプル・サイズが十分に大きいのが前提）。また、実際の被説明変数の値は 0 か 1 しかとらないので、被説明変数は確率が \(p_i\) の二項分布に従うと仮定する。なお上のモデルは説明変数が一つだけだが、複数あっても考え方は同じである。
確率は 0 から 1 までの間の値しかとらないが、\(\frac{\hat p_i}{1- \hat p_i}\) は \(p_i\) が 1 に近づくにつれていくらでも大きな値をとる。この対数をとると、上限がなくなったことに加えて、さらに \(p_i\) が 0 にちかづくにつれていくらでも小さい値をとるので、モデルからの予測値が絶対にとりえない値をとってしまうという問題は回避される（図2 を参照）。
p <- (1:99)/100
odds <- p/(1-p)
logit <- log(odds)
par(mfrow=c(2,2), mgp=c(2, 0.7, 0), mar=c(3,3,2.5,0.2))
plot(p[1:91], odds[1:91], type="l", xlab="p",
main="p と p/(1-p) の関係", ylab="p/(1-p)")
plot(p, logit, type="l", main="p と log(p/(1-p)) の関係", ylab="p/(1-p)")
plot(logit, p, type="l", xlab="x",
main="ロジスティック回帰分析における \n x と p の関係の例")
図2: 確率とオッズ、ロジット、独立変数との関係
式 (1) のような関係が成り立つときに x によって p がどう変化するかを表したのが、図2 の下のグラフで、p が 0 の近くでは変化は緩やかだが、p が 0.5 のあたりにかけて変化は激しくなり、また p が 1 に近づくにつれて緩やかに変化するようになる。線形回帰モデルでは x と y の関係は単純な線形であるが、ロジスティック回帰モデルでは、x と p の関係は、基本的に非線形である。
1.3 glm: R でロジスティック回帰分析
R でロジスティック回帰分析をする場合、glm() という関数を使い、
glm(モデル式, data=データフレーム名, family=binomial)
とすればよい。ほとんど lm() と同じ書式であるが、 “family=binomial” を指定する点だけが異なる。以下は計算例。
glm1 <- glm(y ~ x, family=binomial)
summary(glm1)##
## Call:
## glm(formula = y ~ x, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.03115 -0.11839 0.01144 0.11838 2.29898
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.99121 1.28526 -6.996 2.64e-12 ***
## x 0.64224 0.08999 7.137 9.55e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 412.88 on 299 degrees of freedom
## Residual deviance: 102.38 on 298 degrees of freedom
## AIC: 106.38
##
## Number of Fisher Scoring iterations: 7pred1 <- predict(glm1, type="response")
plot(y ~ jitter(x), ylim=c(-0.1, 1.2))
abline(lm1, col="blue", lty=2)
lines(x, pred1, col="red", lwd=1.5)
legend("bottomright", col=c("blue", "red"), lty=c(2,1), legend=c("線形回帰分析の予測値", "ロジスティック回帰分析の予測値"))
図3: 散布図に5年刻みの習得率と線形回帰の予測値、ロジスティック回帰の予測値を加えたもの
図3 を見ると、単純な線形モデルで予測するよりも、ロジスティック回帰のほうが実際の習得率に近い予測ができているのがわかるだろう。
summary() の結果も lm() の summary() とほとんど同じであるが、異なる点だけ、解説しておく。Deviance Residuals とは、実際の Y の値（0 か 1）とモデルからの予測確率 \(\hat p_i\)の差を単純に計算しているのではなく、 \[ \sqrt{2 \times (Y_i \log \frac{Y_i}{\hat p_i} + (1-Y_i) \log \frac{1 - Y_i}{1 - \hat p_i})} \] を計算し、\(Y_i - \hat p_i < 0\) のときだけマイナス１を掛け合わせたものである。これは単純な残差と同じように、実際の Y の値と予測された確率の差が小さいほど小さくなるが、モデルの尤度に対する影響の大きさ（逸脱度残差の絶対値が大きいほど尤度を下げる）を示す。これも極端に大きかったり、小さかったりすると外れ値と言える。deviance は逸脱度と訳されることが多いと思うが、回帰分析における残差平方和と同じように、モデルと実際のデータの分布の乖離の度合いを示す指標で、 \[ -2 \times \log \mathrm{尤度} \] で計算できる。尤度とは、このモデルが正しい場合に今のサンプルが得られる確率であり、大きいほどよい推定値・モデルであると考えられる。最尤法はこの尤度が最大になるようにパラメータの値を決める。尤度の計算は割愛。
glm() の summary() には Null deviance と Residual deviance が表示されるが、Null deviance とは切片のみのモデルから計算した逸脱度であり、Residual deviance がそのモデルの逸脱度である。説明変数に効果があれば必ず Residual deviance のほうが小さくなり、その落差が大きいほど誤差が減少したと考えられ、当てはまりのよいモデルであると考えられる。それゆえ、逸脱度を使って決定係数のようなものを計算することにも意味はある。McFadden の疑似決定係数とはまさにそれであり、 \[ R^2_\mathrm{McFadden} = 1 - \frac{\mathrm{Residual \; deviance}}{\mathrm{Null \; deviance}} \] で計算される。
2 Fixed Effects Model: ロジスティック回帰分析の固定効果モデルとランダム効果モデル
パネルデータでも二値変数を扱うことはあり、二値変数を被説明変数として固定効果モデルを推定したい場合はある。その場合、モデルは、 \[ \log \frac{\hat p_{it}}{1- \hat p_{it}} = \alpha_i + \beta_1 X_{it} \qquad (2) \] と書ける。このモデルを推定するのに、単純な線形モデルの時のように個体内偏差や階差をとることにはできない。被説明変数が対数変換されているうえに確率のオッズまでとられているからである。そこで、まず思いつくのは、単純な線形モデルのときにそうしたように、切片をすべて愚直に推定するという方法である。これは計算できるが、推定結果がゆがむことが知られている。
そこで、いくつかの方法が提唱されているが、ここでは条件付きの尤度を最小化する方法と、ハイブリッド・モデルを使った方法を紹介しよう。
2.1 Condtional Likelihood 条件付き最尤法で固定効果モデルの推定
条件付き尤度とは、当該のモデルが正しいというだけでなく、個体 i に関して \(Y_{it}=1\)になった回数は所与として、その条件の下で現在のサンプルが得られる確率を計算する。計算の詳細は省略するが、R では survival パッケージの clogit() 関数で
clogit(モデル式, data= データフレーム名) # モデル式の最後に “+ strata(個体ID)” をつける
とすることで計算できる。モデル式の最後に “+ strata(個体ID)” をつける以外は、lm() と同じように引数を指定すればよい。以下は計算例。
head(d1) # x の真の回帰係数は 0.5 になるように作った人工的データ## id t x y
## 1 1 1 -2.108865 0
## 2 1 2 -2.106576 1
## 3 1 3 -1.511525 1
## 4 1 4 -2.229997 1
## 5 1 5 -2.767079 1
## 6 2 1 -2.452040 0library(survival)
cl1 <- clogit(y ~ x + strata(id), data=d1)
summary(cl1)## Call:
## coxph(formula = Surv(rep(1, 1000L), y) ~ x + strata(id), data = d1,
## method = "exact")
##
## n= 1000, number of events= 742
##
## coef exp(coef) se(coef) z Pr(>|z|)
## x 0.61515 1.84993 0.09961 6.176 6.58e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## x 1.85 0.5406 1.522 2.249
##
## Rsquare= 0.042 (max possible= 0.424 )
## Likelihood ratio test= 43.07 on 1 df, p=5.29e-11
## Wald test = 38.14 on 1 df, p=6.583e-10
## Score (logrank) test = 41.46 on 1 df, p=1.2e-10summary() の結果はやはり lm() とほとんど同じだが、異なる点を解説しよう。 “number of events” は被説明変数が 1 になっているケース数である。 “exp(coef)” は係数の指数をとった値で、ロジスティック回帰分析の場合、被説明変数が 1 増えると、被説明変数が 1 になる確率のオッズ (\(\frac{\hat p_i}{1- \hat p_1}\)）が exp(coef) 倍になる（証明は省略）。この例では、x の値が 1 ふえると y が 1 になるオッズは約 1.7 倍になるということである。“exp(-coef)” は x が 1 減少した時にオッズが何倍になるかを示しており、“lower.95”, “upper.95” は “exp(-coef)” の 95% 信頼区間を示している。“Rsquare” 逸脱度を使った擬似決定係数の一種のようだが、詳細は不明。“Likelihood ratio test” の値は、切片のみのモデルと当該のモデルの逸脱度の差である。これは「当該モデルの説明変数の係数はすべて 0 である」、という帰無仮説のもとで、自由度が説明変数の数のカイ二乗分布をすることが知られている。その検定結果が示されているのだが、この例の場合は、自由度が 1 で p が 0.001 未満なので、すべての説明変数の係数が 0 という帰無仮説は棄却される。Wald test も Score (logrank) test も同じ帰無仮説を別の統計量を使って検定したものである。
以下は、愚直にすべての切片の係数を通常の最尤法で推定した場合と、条件付き最尤法で推定した場合の x の係数を比較した結果である。
# 比較のために通常の最尤法でも推定してみる
glm2 <- glm(y ~ -1 + factor(id) + x, data=d1, family="binomial")
library(texreg)
screenreg(list(glm2, cl1), omit.coef = "factor", custom.note = "%stars. \n 切片の係数は省略")##
## ========================================
## Model 1 Model 2
## ----------------------------------------
## x 0.78 *** 0.62 ***
## (0.11) (0.10)
## ----------------------------------------
## AIC 1175.81 511.22
## BIC 2162.26
## Log Likelihood -386.90
## Deviance 773.81
## Num. obs. 1000 1000
## R^2 0.04
## Max. R^2 0.42
## Num. events 742
## Missings 0
## ========================================
## *** p < 0.001, ** p < 0.01, * p < 0.05.
## 切片の係数は省略やはり条件付き最尤法のほうが、真の傾きである 0.5 に近いことがわかるだろう（この d1 というデータは x の真の傾きが 0.5 になるように作ってある）。
2.2 Hybrid: ハイブリッド・モデル
ハイブリッド・モデルはすでに線形の固定効果モデルで紹介したものと同じように、時変変数の固体内偏差をとり、それをランダム効果モデルに投入すればよい。それゆえモデルを式で表すと、 \[ \log \frac{\hat p_{it}}{1- \hat p_{it}} = \beta_0 + \beta_1 (X_{it} - \bar X_{i \cdot}) + \beta_2 Z_i + \alpha_i \qquad (3) \] となる。ただし、\(\bar X_{i \cdot}\) は時変変数の固体内平均、\(Z_i\)は時定変数である。上の式は時変変数と時定変数が1つずつ説明変数として投入されているが、複数あっても考え方は同じである。\(\beta_1\) が固定効果モデルの係数と一致する。時定変数もモデルに投入してよいが、投入してもしなくても、時変変数の固体内偏差の係数に変化はない。
R で計算する場合、lme4 パッケージの glmer() 関数で、
glmer(モデル式, data=データフレーム名, family=binomial) # モデル式の最後に “+ (1 | 個体ID)” をつける
とすればよい。以下は計算例。
d1$x.m <- ave(d1$x, d1$id) # 個体内平均の計算 pdata.frame で Between()を使った場合と同じ
d1$x.d <- d1$x - d1$x.m # 固体内偏差の計算
library(lme4)
hybrd1 <- glmer(y ~ x.d + (1 | id), data=d1, family=binomial) # 固定効果モデルと同じ
hybrd2 <- update(hybrd1, ~. + x.m) # 時定変数もおまけで追加
screenreg(list(cl1, hybrd1, hybrd2), single.row = TRUE)##
## ====================================================================================
## Model 1 Model 2 Model 3
## ------------------------------------------------------------------------------------
## x 0.62 (0.10) ***
## (Intercept) 1.23 (0.10) *** 1.24 (0.09) ***
## x.d 0.58 (0.09) *** 0.57 (0.09) ***
## x.m 0.58 (0.08) ***
## ------------------------------------------------------------------------------------
## AIC 511.22 1097.86 1048.47
## R^2 0.04
## Max. R^2 0.42
## Num. events 742
## Num. obs. 1000 1000 1000
## Missings 0
## BIC 1112.58 1068.10
## Log Likelihood -545.93 -520.23
## Num. groups: id 200 200
## Variance: id.(Intercept) 0.53 0.16
## Variance: Residual 1.00 1.00
## ====================================================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05上の結果を見ると、条件付き最尤法でもハイブリッド・モデルでもほぼ同じ結果が得られるのがわかる。Model 3 の固体内偏差と固体内平均の係数を比較すると、若干、固体内平均の係数が大きくなる。TCUH と固体内平均に相関があるほど、この乖離は大きくなる。
3 Ordered: 順序ロジットと多項ロジット
従属変数が複数のカテゴリからなる場合、順序ロジットか多項ロジットというモデルでの推定が考えられる。順序ロジットは条件付き最尤法ができず、多項ロジットはできるが、データの準備が面倒なので、ここではハイブリッド・モデルを使った方法だけを紹介する。
順序ロジット・モデルにも幾つか集類があるが、もっともよく使われているのは、 \[ \log \frac{\mathrm{P}(Y_i \leq j)}{\mathrm{P}(Y_i > j)} =\beta_{0j} - \beta_1 X_i \] というものである。ただし、\(\mathrm{P}(A)\) で A が真である確率であり、被説明変数が \(J\) 個のカテゴリからなるとすると、\(j=1, \; \ldots, \; J-1\) で、\(j\) は被説明変数のカテゴリを指している。上のモデルでは説明変数は一つだが、複数あってもよい。
例えば、被説明変数の値が、「下」、「中の下」、「中の上」、「上」の4カテゴリからなり、一つ の説明変数で予測するとすると、 \[ \log \frac{\mathrm{P}(Y_i \leq \mathrm{下})}{\mathrm{P}(Y_i > \mathrm{下})} \, =\beta_{0\mathrm{下}} \,- \beta_1 X_i \] \[ \log \frac{\mathrm{P}(Y_i \leq \mathrm{中の下})}{\mathrm{P}(Y_i > \mathrm{中の下})} \, =\beta_{0\mathrm{中の下}} \, - \beta_1 X_i \] \[ \log \frac{\mathrm{P}(Y_i \leq \mathrm{中の上})}{\mathrm{P}(Y_i > \mathrm{中の上})} \, =\beta_{0\mathrm{中の上}} \, - \beta_1 X_i \] という3つの式をたてて、係数を同時に推定する。注意すべき点は、\(Y_i\)が小さな値をとる確率を予測しているが、説明変数の係数の符号がマイナスなので、けっきょく\(\beta_1\)がプラスの値をとるならば、説明変数の値が大きいほど被説明変数の値も大きくなりやすい、と解釈できる。また、切片は一つ一つの式で異なるが、説明変数の傾きは共通である（これを平行仮定ということがある）。これらのパラメータの制約をゆるめたり、強めたりすることも、ソフトウェアによっては可能である。
3.1 clm(), polr()
まず、横断データ用の通常の順序ロジット用の関数はいくつかあるが、ordinal パッケージの clm() と MASS パッケージの polr() を紹介しておく。
clm(モデル式, data=データフレーム名) # 被説明変数は因子でなければならない
polr(モデル式, data=データフレーム名) # 被説明変数は因子でなければならない
以下は計算例。
# 架空のデータの作成
y01 <- rbinom(N*T, 4, p)
d2 <- data.frame(id, t, x, y=y01)
d2$x.m <- ave(d2$x, d2$id) # 個体内平均の計算 pdata.frame で Between()を使った場合と同じ
d2$x.d <- d2$x - d2$x.m # 固体内偏差の計算. pdata.frame で Within() を使った場合と同じ
library(ordinal)
ord1 <- clm(factor(y) ~ x, data=d2) # プールデータで分析
summary(ord1)## formula: factor(y) ~ x
## data: d2
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 1000 -1101.20 2212.41 5(0) 5.79e-07 3.6e+01
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## x 0.89371 0.05288 16.9 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 0|1 -4.78037 0.25336 -18.868
## 1|2 -3.12613 0.13827 -22.609
## 2|3 -1.40936 0.08541 -16.502
## 3|4 0.52756 0.07278 7.249link が logit になっているが、通常の順序ロジット・モデルの推定値であることを示している。これは link=“probit” のような引数を指定すればプロビットでも推定できる。threshold とはふつう切片とか定数とか呼んでいるもののことであるが、順序ロジットの場合、被説明変数の数を \(k\) とすれば \(k - 1\) 個の切片を推定することになる。この切片の値に対して制約をかける（例えば、等間隔にするなど）ことができるが、デフォルトでは特に制約をかけない。この切片に制約をかけないモデルを “flexible” とこのプログラムでは呼んでいる。“nobs” は number of observations, logLik は対数尤度 (log likelihood)、niter 以下は計算のプロセスのログのようである。Threshold の “0 | 1” は被説明変数が 0 以下か、1以上かを予測する式の切片であり、以下も同様。
以下は polr() での計算例。
library(MASS)
ord2 <- polr(factor(y) ~ x, data=d2) # プールデータで分析
summary(ord2)## Call:
## polr(formula = factor(y) ~ x, data = d2)
##
## Coefficients:
## Value Std. Error t value
## x 0.8937 0.05288 16.9
##
## Intercepts:
## Value Std. Error t value
## 0|1 -4.7804 0.2534 -18.8680
## 1|2 -3.1261 0.1383 -22.6087
## 2|3 -1.4094 0.0854 -16.5021
## 3|4 0.5276 0.0728 7.2488
##
## Residual Deviance: 2202.406
## AIC: 2212.406library(texreg)
screenreg(ord2)##
## ============================
## Model 1
## ----------------------------
## x 0.89 ***
## (0.05)
## ----------------------------
## AIC 2212.41
## BIC 2236.94
## Log Likelihood -1101.20
## Deviance 2202.41
## Num. obs. 1000
## ============================
## *** p < 0.001, ** p < 0.01, * p < 0.05summary() では p 値が出力されないが、screenreg() では * で有意かどうかは示される。
3.2 clmm()
ランダム効果付きの順序ロジットの関数としては、clmm を紹介しておく。
clmm(モデル式, data=データフレーム名) # 被説明変数は因子、モデル式の最後に “+ (1 | 個体ID)”
以下は計算例
cl1 <- clmm(factor(y) ~ x.d + (1 | id), data = d2)
summary(cl1)## Cumulative Link Mixed Model fitted with the Laplace approximation
##
## formula: factor(y) ~ x.d + (1 | id)
## data: d2
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 1000 -1193.73 2399.46 378(1137) 2.32e-03 2.6e+01
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.089 1.044
## Number of groups: id 200
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## x.d 0.74600 0.07545 9.887 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 0|1 -4.7335 0.2699 -17.536
## 1|2 -3.0891 0.1621 -19.057
## 2|3 -1.3990 0.1152 -12.143
## 3|4 0.5356 0.1045 5.125“Random Effects:” の欄の数値がランダム効果の分散と標準偏差、Number of Groups は、個体数。パネルデータの分析では、ランダム効果で個人の時間的に不変の観察されない異質性を表すが、一般的にはサンプルをいくつかのグループに分けて、そのグループに共通の残差をランダム効果で表すので、Group という言い方をするほうが多いだろう。
3.3 Multinominal Logit
被説明変数のカテゴリに順序がない場合は、以下のような多項ロジットと呼ばれるモデルが、よく使われている。 \[ \log \frac{\mathrm{P}(Y_i = j)}{\mathrm{P}(Y_i =1)} =\beta_{0j} + \beta_{1j} X_i \] このモデルでは1番目のカテゴリを基準カテゴリとして、その他の値をとる確率を\(X_i\)の値で予測している。例えば、被説明変数がリンゴ、ミカン、バナナのうちどれが一番好きか選んでもらった結果であるとする。多項ロジットモデルは、最初のリンゴを基準カテゴリとすると \[ \log \frac{\mathrm{P}(Y_i = \mathrm{ミカン})}{\mathrm{P}(Y_i =\mathrm{リンゴ})} \, =\beta_{0\mathrm{ミカン}}\, + \beta_{1\mathrm{ミカン}} X_i \] \[ \log \frac{\mathrm{P}(Y_i = \mathrm{バナナ})}{\mathrm{P}(Y_i =\mathrm{リンゴ})} =\beta_{0\mathrm{バナナ}} + \beta_{1\mathrm{バナナ}} X_i \] という式をたて、リンゴと比べたミカンとバナナの選ばれやすさを予測するモデルである。多項ロジットモデルでは、切片も傾きも特に制約をかけずに推定するのが一般的である。
3.3.1 multinom()
横断データ用の多項ロジットの場合は、nnet パッケージの mulitinom() 関数
multinom(モデル式, data=データフレーム名) # 被説明変数は因子でなければならない
library(nnet)
ml1 <- multinom(factor(y) ~ x, data=d2)## # weights: 15 (8 variable)
## initial value 1609.437912
## iter 10 value 1105.767119
## final value 1090.806960
## convergedsummary(ml1)## Call:
## multinom(formula = factor(y) ~ x, data = d2)
##
## Coefficients:
## (Intercept) x
## 1 2.014693 0.4492537
## 2 3.893295 1.0448151
## 3 4.920016 1.7029887
## 4 4.885367 2.2534252
##
## Std. Errors:
## (Intercept) x
## 1 0.5908531 0.2545404
## 2 0.5619355 0.2469941
## 3 0.5589975 0.2515546
## 4 0.5592132 0.2559868
##
## Residual Deviance: 2181.614
## AIC: 2197.614screenreg(ml1)##
## ======================================================================
## Model 1 Model 2 Model 3 Model 4
## ----------------------------------------------------------------------
## (Intercept) 2.01 *** 3.89 *** 4.92 *** 4.89 ***
## (0.59) (0.56) (0.56) (0.56)
## x 0.45 1.04 *** 1.70 *** 2.25 ***
## (0.25) (0.25) (0.25) (0.26)
## ----------------------------------------------------------------------
## AIC 2197.61 2197.61 2197.61 2197.61
## BIC 2236.88 2236.88 2236.88 2236.88
## Log Likelihood -1090.81 -1090.81 -1090.81 -1090.81
## Deviance 2181.61 2181.61 2181.61 2181.61
## Num. obs. 1000 1000 1000 1000
## ======================================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05summary() では Z 得点も p 値も計算されないが、screenreg() を使えば検定結果を * であらわしてくれる。
3.3.2 ランダム効果付きの多項ロジット
R の場合、以下の関数がランダム効果付きの多項ロジット・モデルに近いことをやってくれるが、どれも微妙に異なり、パネルデータのハイブリッド・モデルの推定には使いにくい。
- mlogit パッケージの mlogit()、
- mclogit パッケージの mclogit()
- mixcat パッケージの npmlt()
- ordinal パッケージの clmm2()
ランダム効果付きの多項ロジット・モデルは、SAS, HLM, Stata のような商用ソフトを使えば比較的簡単にできる。
4 Count Data: 度数が被説明変数の場合
被説明変数が、回数や度数、人数を表す場合（例えば転職の回数や友人数）、正規分布を仮定した線形回帰ではなく、被説明変数がポアソン分布や負二項分布をすると仮定した回帰モデルが用いられることがある。よく用いられるのは、 \[ \log\mathrm{E}(Y_i) =\beta_0 + \beta_1 X_i \] というモデルである。ただし、\(\mathrm{E}(Y_i)\) は \(Y_i\) の期待値である。ポアソン分布も負二項分布も 0 以上の整数の値をとる確率変数の分布なので、正規分布を仮定するよりは、度数データの分析に適していると考えられる。
4.1 Poisson: ポアソン分布
\(x\) という確率変数が \(j\)（\(j\) は 0 以上の整数）という値をとる確率を \(\mathrm{P}(x = j)\) とすると、ポアソン分布は以下の式で定義される。 \[ \mathrm{P}(x = j) = \frac{\lambda^j \exp(-\lambda)}{j!} \] ポアソン分布する確率変数 (\(x\)) の期待値も分散も \(\lambda\) になり、生起確率は非常に低いが試行回数は非常に多い二項分布に近似することが知られている。それゆえ低確率だが試行回数の多い出来事の回数（例えば死亡事故）はポアソン分布に従うと言われている。
負二項分布の式は省略するが、分布形はポアソン分布を横に押し広げたようなかたちである（図4 を参照）。
図4: ポアソン分布と負二項分布の確率（密度？）関数
ポアソン分布は、期待値と分散の大きさが必ず一致する分布であるが、実際のデータでは、残差（ロジスティック回帰分析やポアソン回帰分析のような一般化線形モデルでは、あまり「残差」という言い方はしないが、モデルから予測される \(Y_i\) の期待値と実際の \(Y_i\) の差を残差とここでは呼んでいる）の分散は \(Y_i\) の期待値よりも大きくなることもある。そのような状況は overdispersion と呼ばれ、標準誤差が過小に推定されてしまうことが知られている。そこで overdispersion が生じている場合は、\(Y_i\) が負二項分布をしていると仮定すると問題が解消すると言われている。負二項分布は正規分布のように、期待値と分散は別々のパラメータで指定される（つまり両者は一致しない）ので、データから残差のバラつきを推定することが可能であり、overdispersion の問題は回避できる。
どちらの分布を仮定してもパラメータの推定値は同じで標準誤差が違うだけであり、負二項回帰分析には overdispersion の心配がないのでポアソン回帰は最初から使わず、負二項分布を仮定して回帰分析すればいいという考え方もある。ただし、負二項回帰分析のほうが、残差のバラつき (dispersion parameter) を推定するために一つ余分にパラメータを推定するし、やや計算に時間がかかる気がする。またデータがポアソン分布していると、エラーが出てうまく推定できないこともあった。両方やってみて、ほぼ同じ結果ならばポアソン回帰のほうが簡潔でよいかもしれない。判断に迷う場合は、尤度比検定をすればよい。dispersion parameter をゼロに固定すると、負二項分布とポアソン分布は一致するので、ポアソン回帰と負二項回帰は入れ子の関係にあり、尤度比検定ができる。
4.2 glm(), family=poisson
R でポアソン回帰分析をする場合は、以下のように glm() で “family=poisson” を引数に指定する。
glm(モデル式, data = データフレーム名, family = poisson)
以下は計算例。
head(d4) # データの中身をのぞいてみる## x y
## 1 0.01 0
## 2 0.02 0
## 3 0.03 0
## 4 0.04 0
## 5 0.05 0
## 6 0.06 0plot(y ~ x, data = d4) # 分布を確認
pois1 <- glm(y ~ x, data = d4, family=poisson)
summary(pois1)##
## Call:
## glm(formula = y ~ x, family = poisson, data = d4)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.7334 -1.0313 -0.5561 0.6151 4.7713
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.91077 0.10002 -19.10 <2e-16 ***
## x 0.97461 0.02426 40.17 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 3452.18 on 499 degrees of freedom
## Residual deviance: 962.75 on 498 degrees of freedom
## AIC: 2034.1
##
## Number of Fisher Scoring iterations: 5library(car)
Boot(pois1) # ブートストラップで標準誤差の歪みを確認##
## CASE RESAMPLING BOOTSTRAP FOR CENSORED DATA
##
##
## Call:
## boot::boot(data = data.frame(update(object, model = TRUE)$model),
## statistic = boot.f, R = R, .fn = f)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* -1.910773 -4.293382e-04 0.13207082
## t2* 0.974608 1.987638e-05 0.03540865glm() の summary() の結果はロジスティック回帰分析の場合と同じである。上の架空のデータは切片が -2 、傾きが 1 で overdispersion が生じるように人工的に作ってあるのだが、係数は正しく推定されているのがわかる。overdispersion がある場合、標準誤差が過小に推定されてしまうが、そのことを確認するためにブートストラップ法で標準誤差を計算しなおしている（ブートストラップ法は残差の分布に関する仮定に依存しない頑健な推定値を返してくれるが、計算に時間が掛かる）。“t1*" が1番目の係数（この場合は切片）、“t2*" が2番目の係数（この場合は x の傾き）の推定結果である。これを見ると、やはりブートストラップで推定した標準誤差のほうがやや大きいのがわかる。
4.3 glm.nb()
負二項分布回帰をする場合は、glm.nb() 関数を使う。引数は以下のように lm() と同じである。
glm.nb(モデル式, data = データフレーム名)
以下は計算例。
nb1 <- glm.nb(y ~ x, data=d4)
summary(nb1)##
## Call:
## glm.nb(formula = y ~ x, data = d4, init.theta = 3.923611939,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.7679 -0.8759 -0.4212 0.4432 2.5935
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.96722 0.12610 -15.60 <2e-16 ***
## x 0.99021 0.03398 29.14 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(3.9236) family taken to be 1)
##
## Null deviance: 1751.51 on 499 degrees of freedom
## Residual deviance: 483.31 on 498 degrees of freedom
## AIC: 1818.2
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 3.924
## Std. Err.: 0.558
##
## 2 x log-likelihood: -1812.220summary() の出力は glm() の summary() の出力とほとんど同じであるが、最後に “Theta” とその標準誤差（Std. Err）が出力される。Theta は残差のバラつきの指標で dispersion parameter と呼ばれる。x の係数は約 1 で、おおむね正解が推定されており、その点はポアソン回帰分析の場合と同じである。ポイントは x の標準誤差で、ポアソン回帰のそれよりも大きく、ポアソン回帰のブートストラップ推定の結果に近いことがわかるだろう。
ポアソン回帰と負二項回帰の適合度の尤度比検定をするならば、両者の逸脱度の差が帰無仮説（dispersion parameter = 0、つまりポアソン回帰モデルでよい) のもとで自由度 1 のカイ二乗分布に従うことを利用する。逸脱度は、
deviance(lm(), glm(), glm.nb() などの結果)
で取り出せる。以下は計算例。
dev1 <- deviance(pois1) - deviance(nb1)
1 - pchisq(dev1, df=1) # pchisq() で dev1 の値と自由度 1 に対応する左側の確率を返す## [1] 0この例の場合、p 値は 0.001 未満なので、負二項回帰モデルのほうが採択される。
5 Zero Inflated Model
度数データでは、ゼロが特別な意味を持ち、ポアソン分布や負二項分布の予測よりも高い頻度で観察されることがある。例えば、年賀状をまったく書かない人と、基本的には書く人がいるとする。まったく書かない人は常にゼロであるが、、書く人は状況によって書く枚数が変動するとする。そうすると、一枚も年賀状を書かなかった人の中には、まったく書かない人と、基本的には書くが何らかの事情でたまたまゼロ枚だった人が混じっていることになる。このような場合、Zero Inflated Model が用いられる。モデルは全く書かない人と基本的には書く人を識別する二項ロジット・モデルと、基本的には書く人に関して枚数を予測するモデルを同時推定するものである。以下の pscl パッケージの zeroinfl() 関数で、計算できる。
zeroinfl(モデル式, data=データフレーム名) # 一段階目とに段階目の説明変数が違う場合は、| のあとに一段回目の説明変数を指定する
library(pscl)
example(zeroinfl)##
## zernfl> ## data
## zernfl> data("bioChemists", package = "pscl")
##
## zernfl> ## without inflation
## zernfl> ## ("art ~ ." is "art ~ fem + mar + kid5 + phd + ment")
## zernfl> fm_pois <- glm(art ~ ., data = bioChemists, family = poisson)
##
## zernfl> fm_qpois <- glm(art ~ ., data = bioChemists, family = quasipoisson)
##
## zernfl> fm_nb <- glm.nb(art ~ ., data = bioChemists)
##
## zernfl> ## with simple inflation (no regressors for zero component)
## zernfl> fm_zip <- zeroinfl(art ~ . | 1, data = bioChemists)
##
## zernfl> fm_zinb <- zeroinfl(art ~ . | 1, data = bioChemists, dist = "negbin")
##
## zernfl> ## inflation with regressors
## zernfl> ## ("art ~ . | ." is "art ~ fem + mar + kid5 + phd + ment | fem + mar + kid5 + phd + ment")
## zernfl> fm_zip2 <- zeroinfl(art ~ . | ., data = bioChemists)
##
## zernfl> fm_zinb2 <- zeroinfl(art ~ . | ., data = bioChemists, dist = "negbin")5.1 Fixed Effects Poisson Regression
ポアソン回帰や負二項回帰で固定効果モデルを推定するには、R の場合2つの方法がある。
- 一つは glm() や glm.nb() を使い通常の最尤法で愚直にすべての個体の切片を推定する。
- もう一つは、ハイブリッド・モデルで推定する。
ポアソン回帰や負二項回帰の場合、通常の最尤法で愚直にすべての個体の切片を推定しても一致推定量が得られることが知られている。
以下は通常の最尤法で愚直にすべての個体の切片を推定した計算例。
head(d3)## id t x x.m x.d y1
## 1 1 1 -2.108865 -2.144809 0.03594322 0
## 2 1 2 -2.106576 -2.144809 0.03823257 0
## 3 1 3 -1.511525 -2.144809 0.63328377 0
## 4 1 4 -2.229997 -2.144809 -0.08518874 0
## 5 1 5 -2.767079 -2.144809 -0.62227082 0
## 6 2 1 -2.452040 -2.176638 -0.27540146 0plot(y1 ~ x, data=d3)
pois2 <- glm(y1 ~ -1 + factor(id) + x, data=d3, family=poisson)
nb2 <- glm.nb(y1 ~ -1 + factor(id) + x, data=d3)
screenreg(list(pois2, nb2), omit.coef="id")##
## ========================================
## Model 1 Model 2
## ----------------------------------------
## x 0.97 *** 1.00 ***
## (0.05) (0.06)
## ----------------------------------------
## AIC 1779.27 1717.71
## BIC 2765.73 2709.08
## Log Likelihood -688.64 -656.86
## Deviance 674.22 462.66
## Num. obs. 1000 1000
## ========================================
## *** p < 0.001, ** p < 0.01, * p < 0.055.2 Hybrid Poisson
ハイブリッド・モデルも、ロジスティック回帰分析の場合と同じように推定できる。
glmer(モデル式, data = データフレーム名, family = poisson) # 説明変数は固体内偏差、モデル式の最後に +(1 | 個体ID)
glmer.nb(モデル式, data = データフレーム名) # 説明変数は固体内偏差、モデル式の最後に +(1 | 個体ID)
以下はハイブリッド・モデルでの計算例。
library(lme4)
poisr1 <- glmer(y1 ~ x.d +(1 | id), data=d3, family=poisson)
summary(poisr1)## Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
## Family: poisson ( log )
## Formula: y1 ~ x.d + (1 | id)
## Data: d3
##
## AIC BIC logLik deviance df.resid
## 1943.0 1957.7 -968.5 1937.0 997
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9557 -0.4176 -0.2322 -0.1163 10.0865
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 3.582 1.893
## Number of obs: 1000, groups: id, 200
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.97806 0.18699 -10.58 <2e-16 ***
## x.d 0.97036 0.04551 21.32 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## x.d -0.158nbr1 <- glmer.nb(y1 ~ I(x.d/100) +(1 | id), data=d3) # 警告メッセージが出るので説明変数を100で割った
summary(nbr1)## Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
## Family: Negative Binomial(1.665) ( log )
## Formula: y1 ~ I(x.d/100) + (1 | id)
## Data: ..2
##
## AIC BIC logLik deviance df.resid
## 1822.9 1842.5 -907.4 1814.9 996
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5731 -0.5004 -0.2979 -0.1382 6.3009
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.9574 1.3991
## Residual 0.5965 0.7723
## Number of obs: 1000, groups: id, 200
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) -1.9377 0.1717 -11.28 <2e-16 ***
## I(x.d/100) 98.7057 6.8928 14.32 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## I(x.d/100) -0.230以上、駆け足で離散変数の固定効果モデルを紹介した。ランダム効果モデルについては、特に触れていないが、ハイブリッド・モデルで、説明変数を固体内偏差ではなく、そのまま使うだけで推定できる。通常の線形モデルの場合と同様、固定効果モデルとほとんど同じ推定値が得られるならば、ランダム効果モデルのほうが誤差の少ない推定が可能である。
ここでは触れなかったが、pglm というパッケージも公開されており、これを使えばもっと簡単に離散変数の固定効果モデルが推定できるかもしれない。ただ現時点ではどうやって計算しているのかまったく書かれていないので、今回は触れなかった。今後の発展を期待したい。