女性が働いている確率を説明するに当たって、outcome は働いているか、働いていないかというbinary responseであるという特徴から、タイトルの３つのモデルを考えます。
使うデータセットの説明は以下の通りです。
Wooldridge の Introductory Econometrics の公開されているデータセットの中のmrozというデータの内容で、
library(foreign)
mroz<-read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/mroz.dta")
というようにしてRで読み込めます。
使う変数は
inlf: in the labor force (働いているかどうか)
nwifeinc: non wife income (自分で稼ぐ分以外の収入。夫の収入など)
educ: education (教育年数)
exper: experience(就業経験)
age: age (年齢)
kidslt6: kids less than 6 (6歳未満の子供の数)
kidsge6: kids greater than or equal to 6 (6歳より上の歳の子供の数)
LPM , Logit , Probit regression (後者２つはMaximum Liklihood estimation)
LPM (Linear Probability Model)
$$inlf_i=\beta_0+\beta_1nwifeinc_i+\beta_2educ_i+\beta_3exper_i+\beta_4exper^2_i+\beta_5age_i+\beta_6kidslt6_i+\beta_7kidsge6_i+u_i$$
inlfがbinary(0,1変数)の時、OLSの回帰直線に当たるE(inlf)はinlfが1になる確率に一致する。
linprob <- lm(inlf~nwifeinc+educ+exper+I(exper^2)+age+kidslt6+kidsge6,data=mroz)
t-test using heteroscedasticity-robust SE(homoskedasticには構造上なり得ないので、今回はheteroskedsasticity-robust se を使ってt-test)
library(lmtest);library(car)
coeftest(linprob,vcov=hccm)
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.58551922 0.15358032 3.8125 0.000149 ***
nwifeinc -0.00340517 0.00155826 -2.1852 0.029182 *
educ 0.03799530 0.00733982 5.1766 2.909e-07 ***
exper 0.03949239 0.00598359 6.6001 7.800e-11 ***
I(exper^2) -0.00059631 0.00019895 -2.9973 0.002814 **
age -0.01609081 0.00241459 -6.6640 5.183e-11 ***
kidslt6 -0.26181047 0.03215160 -8.1430 1.621e-15 ***
kidsge6 0.01301223 0.01366031 0.9526 0.341123
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1`
どの変数のcoefficientも5%では有意。
例として２つの極端な値で働く確率を予測してみます。
prediction for 2 extreme women
xpred <- list(nwifeinc=c(100,0),educ=c(5,17),exper=c(0,30),
age=c(20,52),kidslt6=c(2,0),kidsge6=c(0,0))
predict(linprob,xpred,type = "response")
1 2
-0.4104582 1.0428084
LPM の問題としては、回帰線は確率を表すのに1を超えたり0を下回ったりするところ。observaitionが端の方では合わない。この欠点を次の２つで克服しにいきます。
Logit model
Logistic distributionの cumurative distribution function で確率を予測します。
$$\frac{1}{1+\exp^-(\beta_0+\beta_1nwifeinc_i+\beta_2educ_i+\beta_3exper_i+\beta_4exper^2_i+\beta_5age_i+\beta_6kidslt6_i+\beta_7kidsge6_i)}$$
この際のパラメータの推定方法はMaximum Likelihood Estimation です。
logistic 分布の累積確率分布や次のprobit model で使う 正規分布の累積確率分布は0~1の範囲でしか動かないという特徴があります。
summary(logitres<-glm(inlf~nwifeinc+educ+exper+I(exper^2)+age+kidslt6+kidsge6,
family=binomial(link=logit),data=mroz))
Call:
glm(formula = inlf ~ nwifeinc + educ + exper + I(exper^2) + age +
kidslt6 + kidsge6, family = binomial(link = logit), data = mroz)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.1770 -0.9063 0.4473 0.8561 2.4032
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.425452 0.860365 0.495 0.62095
nwifeinc -0.021345 0.008421 -2.535 0.01126 *
educ 0.221170 0.043439 5.091 3.55e-07 ***
exper 0.205870 0.032057 6.422 1.34e-10 ***
I(exper^2) -0.003154 0.001016 -3.104 0.00191 **
age -0.088024 0.014573 -6.040 1.54e-09 ***
kidslt6 -1.443354 0.203583 -7.090 1.34e-12 ***
kidsge6 0.060112 0.074789 0.804 0.42154
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1029.75 on 752 degrees of freedom
Residual deviance: 803.53 on 745 degrees of freedom
AIC: 819.53
Number of Fisher Scoring iterations: 4
Log likelihood value の算出
logLik(logitres)
'log Lik.' -401.7652 (df=8)
McFadden's pseudo R-squared
1 - logitres$deviance/logitres$null.deviance
[1] 0.2196814
prediction(just same extreme women as LPM)
predict(logitres, xpred,type = "response")
1 2
0.005218002 0.950049117
ちゃんと極端な説明変数の値でも、LPMと違って0~1の間に収まっていることが確認できる。
Probit model
$$\Phi(\beta_0+\beta_1nwifeinc_i+\beta_2educ_i+\beta_3exper_i+\beta_4exper^2_i+\beta_5age_i+\beta_6kidslt6_i+\beta_7kidsge6_i)$$
logit model は働く確率をlogistic 分布の累積確率関数で予測しましたが、probitは正規分布の累積確率分布で予測します。
summary(probitres<-glm(inlf~nwifeinc+educ+exper+I(exper^2)+age+kidslt6+kidsge6,
family=binomial(link=probit),data=mroz))
Call:
glm(formula = inlf ~ nwifeinc + educ + exper + I(exper^2) + age +
kidslt6 + kidsge6, family = binomial(link = probit), data = mroz)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.2156 -0.9151 0.4315 0.8653 2.4553
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.2700736 0.5080782 0.532 0.59503
nwifeinc -0.0120236 0.0049392 -2.434 0.01492 *
educ 0.1309040 0.0253987 5.154 2.55e-07 ***
exper 0.1233472 0.0187587 6.575 4.85e-11 ***
I(exper^2) -0.0018871 0.0005999 -3.145 0.00166 **
age -0.0528524 0.0084624 -6.246 4.22e-10 ***
kidslt6 -0.8683247 0.1183773 -7.335 2.21e-13 ***
kidsge6 0.0360056 0.0440303 0.818 0.41350
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1029.7 on 752 degrees of freedom
Residual deviance: 802.6 on 745 degrees of freedom
AIC: 818.6
Number of Fisher Scoring iterations: 4
Log likelihood value
logLik(probitres)
'log Lik.' -401.3022 (df=8)
McFadden's pseudo R-squared
1 - probitres$deviance/probitres$null.deviance
[1] 0.2205805
同じく働く確率が0~1の間に収まっているかの検証。
predict(probitres,xpred,type = "response")
1 2
0.001065043 0.959869044
大丈夫でした。LPMの弱点は克服できています。
Likelihood Ratio Test for probit model
exper and age are irrelevant 説の検証
$$ H_0:\beta_3=\beta_4=\beta_5=0\ \ \ \ \ H_1:at\ least\ one\ of \ them \ \neq0$$
restr <- glm(inlf~nwifeinc+educ+ kidslt6+kidsge6,
family=binomial(link=logit),data=mroz)
lrtest(restr,probitres)
Likelihood ratio test
Model 1: inlf ~ nwifeinc + educ + kidslt6 + kidsge6
Model 2: inlf ~ nwifeinc + educ + exper + I(exper^2) + age + kidslt6 +
kidsge6
#Df LogLik Df Chisq Pr(>Chisq)
1 5 -464.92
2 8 -401.30 3 127.25 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
$H_0$は有意水準0.1%でも棄却できるので、experもageもrelevantでありそう。