離散変数の固定効果モデルの例題と練習問題
太郎丸 博
2015年6月22日
1 Example 例題
plm パッケージの Labor Supply データを使って、賃金が子供の数に及ぼす影響を推測するために、固定効果モデルを推定しなさい。
2 Example Answer 例解
2.1 Question 問題
子供数が賃金を高める可能性がある一方、逆に賃金が高いと子供を生み育てやすい、という可能性も考えられる。以下では、LaborSupply データを使って賃金が子供数に及ぼす影響を検討していく。
2.2 Data データと分析方法
LaborSupply データは 1979～1988年の間に行われたパネル調査の結果（有効対象者数が 532人のバランスのとれたパネルデータ）であるが、対象者は男性で 22～60 歳である。
同居している子供数がデータに含まれているが、子供数は図1の右側のグラフのように、30歳代後半までは増加するがその後減少する。これは大学に入学したり、仕事についたり、結婚した子供が親と別居するようになるからであると考えられる。これは親の賃金が減少したからだとは考えにくい。賃金は生むか生まないかを決定するときに影響があるかもしれないが、その後賃金が下がったからといって、自由に子供の数を減らすことはできない。それゆえ、賃金が影響するのは、子供が増えるときであって、減るときではないだろう。このように考えるならば、固定効果モデルで推定して意味があるのは、子供の平均的な数が増加する 22～34歳程度のあいだであろう。
以下ではポアソン回帰と負二項回帰の固定効果モデルを使って、賃金が子供数に及ぼす影響を検討する。主なコントロール変数は年齢で、労働時間や障害の有無も最初のモデルには投入したが、まったく光華がないので、その後のモデルには投入していない。分析はすべてのサンプルと35歳未満のサンプルで行った。
library(plm)
data("LaborSupply")
summary(LaborSupply[, 1:5])## lnhr lnwg kids age disab
## Min. :2.770 Min. :-0.260 Min. :0.000 Min. :22.00 Min. :0.0000
## 1st Qu.:7.580 1st Qu.: 2.370 1st Qu.:1.000 1st Qu.:32.00 1st Qu.:0.0000
## Median :7.650 Median : 2.640 Median :2.000 Median :38.00 Median :0.0000
## Mean :7.657 Mean : 2.609 Mean :1.556 Mean :38.92 Mean :0.0609
## 3rd Qu.:7.780 3rd Qu.: 2.860 3rd Qu.:2.000 3rd Qu.:45.00 3rd Qu.:0.0000
## Max. :8.560 Max. : 4.690 Max. :6.000 Max. :60.00 Max. :1.0000table(LaborSupply$year)##
## 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988
## 532 532 532 532 532 532 532 532 532 532par(mfrow=c(1,2))
plot(jitter(kids) ~ jitter(year), data = LaborSupply)
lm1 <- lm(kids ~ year + I(year^2), LaborSupply)
pred1 <- predict(lm1, data.frame(year=1979:1988))
lines(1979:1988, pred1, col="red") # 二次曲線は予測値を計算してから lines() で描画
plot(jitter(kids) ~ jitter(age), data=LaborSupply)
lm2 <- lm(kids ~ age + I(age^2), LaborSupply)
pred2 <- predict(lm2, data.frame(age=22:60))
lines(22:60, pred2, col="red")
図1: 調査年と子供数、年齢と子供数の散布図（OLS で 予測した二次曲線つき）
2.3 Results 分析結果
下の表の最初のモデルは、すべてのサンプルでポアソン回帰の固定効果モデルを推定した結果である。対数時給の係数は予測通り正の値を示しているが、有意ではない。サンプルを35歳未満に限定し、まったく効果のなかった障害の有無と労働時間はモデルから取り除いたのが、2番目のモデルである。このモデルでは賃金の係数は高まり、有意になっている。同じモデルを負二項回帰モデルでも推定したが、dispersion parameter の推定が収束しなかった。
ポアソン回帰の結果だけでは、標準誤差の正しさに不安が残るので、ブートストラップで賃金の係数の 95% 信頼区間を計算した。confint() 関数は method=“boot” を引数に指定すればブートストラップ法で区間推定してくれるが、パラメトリック・ブートストラップなので、本当にポアソン分布しているかどうか不安な状態ではあまり意味がない。固定効果モデルの結果を Boot() で推定してみると、賃金の係数の標準誤差は 0.056 で有意になるが、これも抽出単位が個人ではない（パーソン・イヤーが抽出単位になる）ので、あまり意味が無いと思う。そこで、自分で個人をリサンプリングするプログラムを書いてブートストラップ推定してみたところ、賃金の係数の95% 信頼区間は、0.101～0.348 で有意である。標準誤差は 0.075 でむしろポアソン回帰の結果よりも小さくなった。ポアソン分布よりもむしろ残差のバラつきは小さいのかもしれない。固定効果モデルで、年齢もモデルに投入しているので、かなり残差は小さいはずである。
以上を総合すると、少なくとも35歳未満の時期には、賃金が子供の数を増やす効果は、統計的に有意であるといえる。
pois1 <- glm(kids ~ -1 + factor(id) + age + I(age^2) + lnwg + disab + lnhr,
data=LaborSupply, family = poisson)
d0 <- subset(LaborSupply, age < 35) # 35歳未満にデータを限定
valid.id <- unique(d0$id)[table(d0$id) > 1] # 2時点以上データがある個人IDを取り出す
d1 <- d0[is.element(d0$id, valid.id), ] # 2時点以上データがある個人だけのデータを作る
# table(d1$id)
table(d1$year)##
## 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988
## 285 285 267 237 204 166 142 119 93 63pois2 <- update(pois1, ~-1 + factor(id) + age + I(age^2) + lnwg, data = d1)
library(MASS)
nb1 <- glm.nb(kids ~ -1 + factor(id) + age + I(age^2) + lnwg, data=d1)## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit
## reached## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit
## reachednb1$theta## [1] 75019.58nb1$SE.theta## [1] 208050.9# ハイブリッド・モデルのためのスクリプト
d1$age.sq <- (d1$age/10) ^ 2
d1$m.age <- ave(d1$age/10, d1$id)
d1$m.age.sq <- ave(d1$age.sq, d1$id)
d1$m.lnwg <- ave(d1$lnwg, d1$id)
d1$d.age <- d1$age/10 - d1$m.age
d1$d.age.sq <- d1$age.sq - d1$m.age.sq
d1$d.lnwg <- d1$lnwg - d1$m.lnwg
# head(d1, 20)
library(lme4)
poisr1 <- glmer(kids ~ d.age + d.age.sq + d.lnwg + (1 | id),
data=d1, family=poisson)
library(texreg)
screenreg(list(pois1, pois2, nb1), omit.coef = "id", digits=3,
custom.model.names = c("ポアソン全サンプル", "ポアソン35歳未満","負二項35歳未満"))##
## ===========================================================
## ポアソン全サンプル ポアソン35歳未満 負二項35歳未満
## -----------------------------------------------------------
## age 0.536 *** 0.647 *** 0.647 ***
## (0.024) (0.131) (0.131)
## I(age^2) -0.007 *** -0.009 *** -0.009 ***
## (0.000) (0.002) (0.002)
## lnwg 0.091 0.236 * 0.236 *
## (0.065) (0.115) (0.115)
## disab 0.034
## (0.067)
## lnhr 0.002
## (0.049)
## -----------------------------------------------------------
## AIC 12757.061 5047.942 5049.980
## BIC 16290.107 6645.785 6653.352
## Log Likelihood -5841.531 -2234.971 -2234.990
## Deviance 1321.767 372.144 372.139
## Num. obs. 5320 1861 1861
## ===========================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05# 負二項回帰の推定が収束しなかったので、ポアソン回帰の係数をブートストラップで再推定
# (cf1 <- confint(poisr1, method="boot")) # 推定はできたがあまり意味なし
# confint(poisr1, method="boot", type="semiparametric", use.u=TRUE) # 意味はあるがエラーが出る
# (cf2 <- Boot(pois2)) # あまり意味が無いがやってみた
##### ブートストラップのためのプログラム ######
# b <- numeric(1000) # 計算結果の係数を格納するベクトル
# for(i in 1:1000){ # 1000回繰り返しのブートストラップ
# id1 <- sample(unique(d1$id), 310, replace=TRUE) # 個人を抽出
# d2 <- NULL # 空のオブジェクト、リサンプリングしたデータフレーム
# for(j in 1:310){
# d2 <- rbind(d2, d1[d1$id == id1[j], ]) # 抽出した個人をデータに加える
# }
# poisr2 <- update(poisr1, data=d2) # リサンプリングしたデータ d2 で poisr1 と同じ分析
# b[i] <- fixef(poisr2)[4] # 賃金の係数だけ取り出す
# }
# sd(b)
# quantile(b, c(0.05, 0.95))思い起してみるに、問題やデータの性質を考えると、イベント・ヒストリー分析のほうが適切なので、ここまで頑張ってブートストラップのプログラムを書くぐらいならば、イベント・ヒストリー分析したほうが生産的だった。まあ、ブートストラップの勉強にはなったが。
3 Exercise 練習問題
plm パッケージの Males データを使って、賃金が結婚に及ぼす影響を固定効果モデルで検討しなさい。