Linear Probability , Logit and Probit Model


女性が働いている確率を説明するに当たって、outcome は働いているか、働いていないかというbinary responseであるという特徴から、タイトルの３つのモデルを考えます。
使うデータセットの説明は以下の通りです。

Wooldridge の Introductory Econometrics の公開されているデータセットの中のmrozというデータの内容で、

```{r}
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になる確率に一致する。

```{r}
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)

```{r}

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
```{r}
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の範囲でしか動かないという特徴があります。

```{r}
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 の算出

```{r}
logLik(logitres) 
```

```
'log Lik.' -401.7652 (df=8)
```

### McFadden's pseudo R-squared
```{r}
1 - logitres$deviance/logitres$null.deviance
```

```
[1] 0.2196814
```

### prediction(just same extreme women as LPM)

```{r}
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は正規分布の累積確率分布で予測します。

```{r}
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

```{r}
logLik(probitres) 
```

```
'log Lik.' -401.3022 (df=8)
```

### McFadden's pseudo R-squared

```{r}
1 - probitres$deviance/probitres$null.deviance
```

```
[1] 0.2205805
```

同じく働く確率が0~1の間に収まっているかの検証。

```{r}
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$$


```{r}
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でありそう。