plm パッケージの Males データを使って結婚が賃金に及ぼす影響を固定効果モデルを使って分析しなさい。
既婚男性と未婚男性を比較すると既婚男性のほうが賃金が高いことは横断データで繰り返し確認されてきた事実であり、結婚が男性の賃金を高める効果は関連する諸変数を統制しても有意である。これは扶養家族手当の類を制度化している会社もあるし、結婚することでモチベーションが高まって、労働生産性が高まり、賃金を上昇させるということかもしれない。しかしこれは、時間によっては変化しない観察されない異質性によって引き起こされた擬似相関かもしれない。例えば、人柄の良さや容姿の美しさが賃金と結婚しやすさに影響を及ぼしているのかもしれない。人柄も容姿も長期的には変化するだろうが、短期的には変化しないと仮定して差し支えないだろうから、固定効果モデルでこれらの効果を統制して、結婚の賃金に対する効果を推定してみよう。
データは plm パッケージの Males というサンプル・データで、ヘルプ・ファイルによれば Vella and Verveek (1998) が用いていたデータらしい。彼らのデータは National Longitudinal Survey (Youth Sample) で、米国で 1980年までに学校を卒業して働き始めた若い男性に関して 1980~1987年のあいだ継続して調査されたものである。545人\(\times 8\)時点のバランスのとれたパネルデータである。
被説明変数は時給の対数 (wage) で、説明変数は労働経験年数 (exper)、団体交渉で賃金が決まっているか (union)、結婚しているかどうか (married)、健康に問題があるか (health)、職業9分類 (occupation) である。プール・データの記述統計は以下の通りである。
library(plm)
data(Males)
names(Males) # 変数名を確認
## [1] "nr" "year" "school" "exper" "union" "ethn" "married"
## [8] "health" "wage" "industry" "occupation" "residence"
d0 <- pdata.frame(Males[, c(-3, -6, -10, -12)]) # 使わない変数を削除
summary(d0[,-1:-2]) # 個人IDの記述統計は要らないし、年次は後で詳しく見るので除外
## exper union married health wage
## Min. : 0.000 no :3296 no :2446 no :4286 Min. :-3.579
## 1st Qu.: 4.000 yes:1064 yes:1914 yes: 74 1st Qu.: 1.351
## Median : 6.000 Median : 1.671
## Mean : 6.515 Mean : 1.649
## 3rd Qu.: 9.000 3rd Qu.: 1.991
## Max. :18.000 Max. : 4.052
##
## occupation
## Craftsmen, Foremen_and_kindred :934
## Operatives_and_kindred :881
## Service_Workers :509
## Clerical_and_kindred :486
## Professional, Technical_and_kindred:453
## Laborers_and_farmers :401
## (Other) :696
xtabs(~year, Males)
## year
## 1980 1981 1982 1983 1984 1985 1986 1987
## 545 545 545 545 545 545 545 545
賃金の個人別の変化をプロットしたのが図1 である。
par(mfrow=c(1,2), mar=c(3,3,2, 2), cex.main=0.7, mgp=c(2, 0.7, 0)) # 図の余白や文字の大きさの調整
plot(wage ~ year, data=d0, main="対数時給")
abline(h=mean(d0$wage), col="red") # abline(h=x) で高さxの水平線を引く
plot(d0$year, diff(d0$wage), main="対数時給の変化") # diff() で一時点前との差分をとってくれる(最初の時点はNA)
abline(h=0, col="red")
text(2, 4.5, paste("1年間の賃金変化の平均 =",
round(mean(diff(d0$wage), na.rm=T),2)
)
)
ほんの少しずつだが、平均的な対数時給が上昇しているのがわかる。図2は婚姻状態の移動をモザイクプロットで示したものである。
par(mfrow=c(3,3), mar=c(3, 3, 1.5, 2), mgp=c(1.5, 0.5, 0))
for(i in 1:7){
plot( married[year==1980+i] ~ married[year==1979+i], data=d0,
xlab="前年に結婚していたか",
ylab="当年に結婚していたか", main=paste(1980+i, "年"), cex.lab=0.7)
}
1980年のときには結婚していたのは2割程度であるが、次第に結婚している者が増え、1987年には6割程度が結婚しているのがわかる。若干の離婚があるものの、結婚がその数を上回っているので。図3は婚姻状態と賃金の関係を箱ひげ図で示したものと、両方の変数の固体内偏差をとって散布図にしたものである。
par(mfrow=c(1,2))
plot(wage ~ married, data=d0, main="婚姻状態と対数時給")
t.test(wage ~ married, data=d0) # 平均値の差の検定
##
## Welch Two Sample t-test
##
## data: wage by married
## t = -14.126, df = 4329.9, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.2508787 -0.1897290
## sample estimates:
## mean in group no mean in group yes
## 1.552436 1.772740
d0$marriage <- as.numeric(d0$married == "yes") # Within は数値ベクトルしか受け付けないので数値のダミー変数に変換
plot(Within(d0$marriage), Within(d0$wage), main="固体内偏差の散布図")
固体内偏差をとっても取らなくても、結婚しているほうが賃金が高い傾向が見られる。
固定効果モデルで結婚の効果を分析した結果が下の表1である。
levels(d0$occupation) <- c("Professionals", "Managers", "Sales", "Clerks", "Craftsmen", "Operatives", "Laborers", "Farm", "Service") # 因子のカテゴリ名が長すぎるので短い名前に変更
pl1 <- plm(wage ~ married, data=d0, model="within")
pl2 <- update(pl1, ~exper + I(exper^2) + union + health + occupation)
pl3 <- update(pl2, ~.+ married)
# 個人内平均を使ったクロスセクションの回帰分析もおまけで
d0$school <- Males$school
d0$ethn <- Males$ethn
pl4 <- update(pl3, ~.+school+ethn, model="between")
library(stargazer)
stargazer(pl1, pl2, pl3, pl4, type="text", title="表1", keep.stat=c("n", "adj.rsq")) # stargazer は2バイト文字をほとんど受け付けないらしい
##
## 表1
## ===========================================================
## Dependent variable:
## --------------------------------------
## wage
## (1) (2) (3) (4)
## -----------------------------------------------------------
## marriedyes 0.243*** 0.045** 0.109***
## (0.018) (0.018) (0.040)
##
## school 0.078***
## (0.012)
##
## ethnblack -0.097**
## (0.047)
##
## ethnhisp 0.032
## (0.041)
##
## exper 0.120*** 0.115*** -0.059
## (0.008) (0.009) (0.050)
##
## I(exper2) -0.004*** -0.004*** 0.006*
## (0.001) (0.001) (0.003)
##
## unionyes 0.084*** 0.083*** 0.341***
## (0.019) (0.019) (0.048)
##
## healthyes -0.014 -0.011 -0.098
## (0.047) (0.047) (0.229)
##
## occupationManagers -0.014 -0.012 0.0001
## (0.032) (0.032) (0.103)
##
## occupationSales -0.064* -0.061 -0.024
## (0.038) (0.038) (0.122)
##
## occupationClerks -0.081*** -0.080*** -0.181*
## (0.031) (0.031) (0.102)
##
## occupationCraftsmen -0.031 -0.029 -0.032
## (0.030) (0.030) (0.084)
##
## occupationOperatives -0.032 -0.028 -0.183**
## (0.031) (0.031) (0.087)
##
## occupationLaborers -0.042 -0.039 -0.404***
## (0.034) (0.034) (0.109)
##
## occupationFarm -0.061 -0.058 -0.698***
## (0.066) (0.066) (0.187)
##
## occupationService -0.045 -0.045 -0.461***
## (0.034) (0.034) (0.089)
##
## Constant 0.884***
## (0.225)
##
## -----------------------------------------------------------
## Observations 4,360 4,360 4,360 545
## Adjusted R2 0.041 0.156 0.157 0.295
## ===========================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
モデル1を見ると、他の変数をコントロールしなければ、\(\exp(0.243)=1.28\)倍だけ結婚している者のほうが賃金が高いことがわかる。モデル3をみると、他の変数をコントロールしても\(\exp(0.045)=1.05\)倍だけ結婚している者のほうが賃金が高い。それゆえ、時間によっては変化しない観察されない異質性とモデルに投入した変数をコントロールしても、結婚は賃金を高めることがあるのがわかった。
おまけで個人内平均をデータとして使って横断的な回帰分析した結果も表示したが、この場合、結婚の係数は 0.109 であるから、この 0.109 のうち、0.064 は時間によって変化しない観察されない異質性によって生じていることがわかる。
plm パッケージの Wages データを使って結婚が女性の賃金に及ぼす影響を固定効果モデルを使って分析しなさい。 ただし、個体IDと時点の変数がなぜかないので、自分で作りなさい。