II Econometric Analysis Using R

Also available in Stata and Python versions

### Chapter 15. Binary Response Model

#### Example 15.1

library(wooldridge)
library(AER)
library(stargazer)
library(pglm)
library(sandwich)
library(mfx)
library(haven)
library(dplyr)
#library(lme4)
library(GJRM)

Married Women´s Labor Force Participation

lpm <- lm(inlf ~ nwifeinc + exper + expersq + educ + age + kidslt6 + kidsge6, data=mroz)
lpmr <- coeftest(lpm, vcovHC(lpm,  type="HC1"))
stargazer(lpm, lpmr, column.labels=c("LPM", "LPM_robust"), no.space=TRUE, type="text")
##
## =======================================================
##                             Dependent variable:
##                     -----------------------------------
##                              inlf
##                               OLS           coefficient
##                                                test
##                               LPM               LPM
##                               (1)               (2)
## -------------------------------------------------------
## nwifeinc                   -0.003**          -0.003**
##                             (0.001)           (0.002)
## exper                      0.039***          0.039***
##                             (0.006)           (0.006)
## expersq                    -0.001***         -0.001***
##                            (0.0002)          (0.0002)
## educ                       0.038***          0.038***
##                             (0.007)           (0.007)
## age                        -0.016***         -0.016***
##                             (0.002)           (0.002)
## kidslt6                    -0.262***         -0.262***
##                             (0.034)           (0.032)
## kidsge6                      0.013             0.013
##                             (0.013)           (0.014)
## Constant                   0.586***          0.586***
##                             (0.154)           (0.152)
## -------------------------------------------------------
## Observations                  753
## R2                           0.264
## Residual Std. Error    0.427 (df = 745)
## F Statistic         38.218*** (df = 7; 745)
## =======================================================
## Note:                       *p<0.1; **p<0.05; ***p<0.01

#### Example 15.2

Married Women´s Labor Force Participation

LPM <- lm(inlf ~ nwifeinc + educ + exper + expersq + age + kidslt6 + kidsge6, data = mroz)
Logit <- glm(inlf ~ nwifeinc + educ + exper + expersq + age + kidslt6 + kidsge6, family=binomial(link="logit"), data = mroz)
Probit <- glm(inlf ~ nwifeinc + educ + exper + expersq + age + kidslt6 + kidsge6, family=binomial(link="probit"), data = mroz)

stargazer(LPM, Logit, Probit,  no.space=TRUE, type="text", title = "Table 15.1 LPM, Logit, and Probit Estimates of Labor Force Participation: (inlf)")
##
## Table 15.1 LPM, Logit, and Probit Estimates of Labor Force Participation: (inlf)
## ===============================================================
##                                 Dependent variable:
##                     -------------------------------------------
##                                        inlf
##                               OLS           logistic   probit
##                               (1)              (2)       (3)
## ---------------------------------------------------------------
## nwifeinc                   -0.003**         -0.021**  -0.012**
##                             (0.001)          (0.008)   (0.005)
## educ                       0.038***         0.221***  0.131***
##                             (0.007)          (0.043)   (0.025)
## exper                      0.039***         0.206***  0.123***
##                             (0.006)          (0.032)   (0.019)
## expersq                    -0.001***        -0.003*** -0.002***
##                            (0.0002)          (0.001)   (0.001)
## age                        -0.016***        -0.088*** -0.053***
##                             (0.002)          (0.015)   (0.008)
## kidslt6                    -0.262***        -1.443*** -0.868***
##                             (0.034)          (0.204)   (0.118)
## kidsge6                      0.013            0.060     0.036
##                             (0.013)          (0.075)   (0.044)
## Constant                   0.586***           0.425     0.270
##                             (0.154)          (0.860)   (0.508)
## ---------------------------------------------------------------
## Observations                  753              753       753
## R2                           0.264
## Log Likelihood                              -401.765  -401.302
## Akaike Inf. Crit.                            819.530   818.604
## Residual Std. Error    0.427 (df = 745)
## F Statistic         38.218*** (df = 7; 745)
## ===============================================================
## Note:                               *p<0.1; **p<0.05; ***p<0.01

#### Example 15.3 pp.587

Testing Exogeniety of Education in the Women´s LFP Model

v2 <- resid(lm(educ ~ nwifeinc + exper + expersq + age + kidslt6 + kidsge6 + motheduc + fatheduc + huseduc, data=mroz))
summary(Probit <- glm(inlf ~ nwifeinc + educ + exper + expersq + age + kidslt6 + kidsge6 + v2, family=binomial(link="probit"), data = mroz))
##
## Call:
## glm(formula = inlf ~ nwifeinc + educ + exper + expersq + age +
##     kidslt6 + kidsge6 + v2, family = binomial(link = "probit"),
##     data = mroz)
##
## Deviance Residuals:
##     Min       1Q   Median       3Q      Max
## -2.2341  -0.9134   0.4252   0.8611   2.3711
##
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.6209080  0.6507213   0.954  0.33999
## nwifeinc    -0.0102850  0.0053524  -1.922  0.05466 .
## educ         0.1035744  0.0406773   2.546  0.01089 *
## exper        0.1262473  0.0190483   6.628 3.41e-11 ***
## expersq     -0.0019431  0.0006023  -3.226  0.00126 **
## age         -0.0543806  0.0086517  -6.286 3.27e-10 ***
## kidslt6     -0.8630815  0.1185043  -7.283 3.26e-13 ***
## kidsge6      0.0313809  0.0443001   0.708  0.47872
## v2           0.0433656  0.0502575   0.863  0.38821
## ---
## 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:  801.85  on 744  degrees of freedom
## AIC: 819.85
##
## Number of Fisher Scoring iterations: 4

#### Example 15.3 pp.589

Endogeniety of Nonwife income in the Women´s LFP Model

v2 <- resid(nwife <- lm(nwifeinc ~  huseduc + exper + expersq + educ + age + kidslt6 + kidsge6, data = mroz))
coeftest(nwife, vcovHC(nwife,  type="HC1"))
##
## t test of coefficients:
##
##                Estimate  Std. Error t value  Pr(>|t|)
## (Intercept) -1.4720e+01  3.6161e+00 -4.0708 5.186e-05 ***
## huseduc      1.1782e+00  1.7019e-01  6.9224 9.588e-12 ***
## exper       -3.1299e-01  1.3616e-01 -2.2987  0.021798 *
## expersq     -4.7756e-04  3.9865e-03 -0.1198  0.904678
## educ         6.7470e-01  2.5156e-01  2.6821  0.007479 **
## age          3.4015e-01  5.4086e-02  6.2891 5.440e-10 ***
## kidslt6      8.2627e-01  1.0441e+00  0.7913  0.428994
## kidsge6      4.3553e-01  2.8859e-01  1.5091  0.131686
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(Probit <- glm(inlf ~ nwifeinc+ exper + expersq + educ  + age + kidslt6 + kidsge6 + v2, family=binomial(link="probit"), data = mroz))
##
## Call:
## glm(formula = inlf ~ nwifeinc + exper + expersq + educ + age +
##     kidslt6 + kidsge6 + v2, family = binomial(link = "probit"),
##     data = mroz)
##
## Deviance Residuals:
##     Min       1Q   Median       3Q      Max
## -2.2523  -0.9078   0.4204   0.8566   2.2803
##
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.0171183  0.5380339   0.032  0.97462
## nwifeinc    -0.0368639  0.0183848  -2.005  0.04495 *
## exper        0.1163118  0.0193869   6.000 1.98e-09 ***
## expersq     -0.0019458  0.0005999  -3.244  0.00118 **
## educ         0.1702142  0.0377615   4.508 6.56e-06 ***
## age         -0.0449529  0.0101351  -4.435 9.19e-06 ***
## kidslt6     -0.8444319  0.1197268  -7.053 1.75e-12 ***
## kidsge6      0.0477912  0.0449431   1.063  0.28761
## v2           0.0267092  0.0191539   1.394  0.16318
## ---
## 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:  800.61  on 744  degrees of freedom
## AIC: 818.61
##
## Number of Fisher Scoring iterations: 4
probitmfx(Probit, data = mroz, atmean = F)
## Call:
## probitmfx(formula = Probit, data = mroz, atmean = F)
##
## Marginal Effects:
##                dF/dx   Std. Err.       z     P>|z|
## nwifeinc -0.01105760  0.00547445 -2.0199 0.0433983 *
## exper     0.03488859  0.00541166  6.4469 1.141e-10 ***
## expersq  -0.00058367  0.00017639 -3.3089 0.0009365 ***
## educ      0.05105699  0.01090377  4.6825 2.834e-06 ***
## age      -0.01348394  0.00293039 -4.6014 4.196e-06 ***
## kidslt6  -0.25329353  0.03248353 -7.7976 6.310e-15 ***
## kidsge6   0.01433531  0.01345082  1.0658 0.2865334
## v2        0.00801162  0.00572509  1.3994 0.1616970
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#### Example 15.4

Women’s Labor Force Participation and Having More than Two Children

df <- read.csv("labsup.csv")
df$agesq <- df$age*df\$age
summary(LPM_OLS <- lm(worked ~ morekids + nonmomi + educ + age + agesq + black + hispan, data=df))
##
## Call:
## lm(formula = worked ~ morekids + nonmomi + educ + age + agesq +
##     black + hispan, data = df)
##
## Residuals:
##     Min      1Q  Median      3Q     Max
## -0.9455 -0.5002  0.2280  0.3924  0.9112
##
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.4486450  0.1649612  -2.720 0.006538 **
## morekids    -0.1091198  0.0054742 -19.933  < 2e-16 ***
## nonmomi     -0.0011675  0.0001354  -8.625  < 2e-16 ***
## educ         0.0206475  0.0009053  22.808  < 2e-16 ***
## age          0.0562704  0.0112458   5.004 5.65e-07 ***
## agesq       -0.0007829  0.0001932  -4.052 5.10e-05 ***
## black        0.0176482  0.0338829   0.521 0.602469
## hispan      -0.1285859  0.0339362  -3.789 0.000152 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4714 on 31849 degrees of freedom
## Multiple R-squared:  0.08175,    Adjusted R-squared:  0.08154
## F-statistic:   405 on 7 and 31849 DF,  p-value: < 2.2e-16
summary(Probit <- glm(worked ~ morekids + nonmomi + educ + age + agesq + black + hispan, family=binomial(link="probit"), data=df))
##
## Call:
## glm(formula = worked ~ morekids + nonmomi + educ + age + agesq +
##     black + hispan, family = binomial(link = "probit"), data = df)
##
## Deviance Residuals:
##     Min       1Q   Median       3Q      Max
## -2.0869  -1.1769   0.7224   0.9868   2.0069
##
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.4964755  0.4520026  -5.523 3.33e-08 ***
## morekids    -0.2986554  0.0150534 -19.840  < 2e-16 ***
## nonmomi     -0.0031412  0.0003719  -8.446  < 2e-16 ***
## educ         0.0554282  0.0025044  22.132  < 2e-16 ***
## age          0.1479387  0.0308659   4.793 1.64e-06 ***
## agesq       -0.0020364  0.0005309  -3.836 0.000125 ***
## black        0.0412891  0.0915437   0.451 0.651966
## hispan      -0.3586388  0.0917648  -3.908 9.30e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
##     Null deviance: 43130  on 31856  degrees of freedom
## Residual deviance: 40436  on 31849  degrees of freedom
## AIC: 40452
##
## Number of Fisher Scoring iterations: 4
probitmfx(Probit, data=df, atmean = F)
## Call:
## probitmfx(formula = Probit, data = df, atmean = F)
##
## Marginal Effects:
##                dF/dx   Std. Err.        z     P>|z|
## morekids -0.10946602  0.00550420 -19.8877 < 2.2e-16 ***
## nonmomi  -0.00113875  0.00013437  -8.4745 < 2.2e-16 ***
## educ      0.02009384  0.00088671  22.6611 < 2.2e-16 ***
## age       0.05363076  0.01117731   4.7982 1.601e-06 ***
## agesq    -0.00073825  0.00019231  -3.8389 0.0001236 ***
## black     0.01502517  0.03343145   0.4494 0.6531200
## hispan   -0.13215805  0.03374329  -3.9166 8.982e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## dF/dx is for discrete change for the following variables:
##
## [1] "morekids" "black"    "hispan"
summary(LPM_IV <- ivreg(worked ~ morekids + nonmomi + educ + age + agesq + black + hispan | samesex + nonmomi + educ + age + agesq + black + hispan, data=df))
##
## Call:
## ivreg(formula = worked ~ morekids + nonmomi + educ + age + agesq +
##     black + hispan | samesex + nonmomi + educ + age + agesq +
##     black + hispan, data = df)
##
## Residuals:
##     Min      1Q  Median      3Q     Max
## -0.9713 -0.4784  0.1867  0.3907  0.9460
##
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.4549690  0.1658195  -2.744 0.006077 **
## morekids    -0.2008320  0.0964721  -2.082 0.037372 *
## nonmomi     -0.0012600  0.0001671  -7.541 4.78e-14 ***
## educ         0.0175522  0.0033755   5.200 2.01e-07 ***
## age          0.0603517  0.0120811   4.996 5.90e-07 ***
## agesq       -0.0008178  0.0001975  -4.140 3.48e-05 ***
## black        0.0168118  0.0340432   0.494 0.621424
## hispan      -0.1308112  0.0341655  -3.829 0.000129 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4735 on 31849 degrees of freedom
## Multiple R-Squared: 0.07365, Adjusted R-squared: 0.07345
## Wald test: 345.9 on 7 and 31849 DF,  p-value: < 2.2e-16
# Bivariate probit
eq1 <- worked ~ morekids + nonmomi + educ + age + agesq + black + hispan
eq2 <- morekids ~ samesex + nonmomi + educ + age + agesq + black + hispan
f.list <- list(eq1, eq2)
mr <- c("probit", "probit")
# library(GJRM)
summary(biprobit <- gjrm(f.list, data=df, Model="B", margins=mr))
##
## COPULA:   Gaussian
## MARGIN 1: Bernoulli
## MARGIN 2: Bernoulli
##
## EQUATION 1
## Link function for mu.1: probit
## Formula: worked ~ morekids + nonmomi + educ + age + agesq + black + hispan
##
## Parametric coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.4753259  0.4496307  -5.505 3.69e-08 ***
## morekids    -0.7024746  0.2040666  -3.442 0.000577 ***
## nonmomi     -0.0034902  0.0003950  -8.836  < 2e-16 ***
## educ         0.0405659  0.0085401   4.750 2.03e-06 ***
## age          0.1632221  0.0312418   5.224 1.75e-07 ***
## agesq       -0.0021523  0.0005277  -4.078 4.53e-05 ***
## black        0.0367334  0.0910000   0.404 0.686460
## hispan      -0.3614829  0.0912099  -3.963 7.40e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## EQUATION 2
## Link function for mu.2: probit
## Formula: morekids ~ samesex + nonmomi + educ + age + agesq + black + hispan
##
## Parametric coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.5725583  0.4514336  -3.483 0.000495 ***
## samesex      0.1446575  0.0144320  10.023  < 2e-16 ***
## nonmomi     -0.0027063  0.0003685  -7.345 2.06e-13 ***
## educ        -0.0907147  0.0024968 -36.333  < 2e-16 ***
## age          0.1190243  0.0307613   3.869 0.000109 ***
## agesq       -0.0010280  0.0005284  -1.946 0.051693 .
## black       -0.0277802  0.0921479  -0.301 0.763053
## hispan      -0.0690522  0.0922843  -0.748 0.454306
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## n = 31857  theta = 0.254(0.0345,0.432)  tau = 0.164(0.022,0.284)
## total edf = 17

#### Example 15.5

Panel Data Models for Women’s Labor Force Participation

#library(haven)
lfpP <- pdata.frame(lfp, index = c("id", "period"))

summary(FE_Linear <- plm(lfp ~ kids + lhinc + educ + black + age + agesq + factor(period), data=lfpP, model = "within"))
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = lfp ~ kids + lhinc + educ + black + age + agesq +
##     factor(period), data = lfpP, model = "within")
##
## Balanced Panel: n = 5663, T = 5, N = 28315
##
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max.
## -0.8410962 -0.0079438  0.0015233  0.0080867  0.8438588
##
## Coefficients:
##                   Estimate Std. Error t-value  Pr(>|t|)
## kids            -0.0388976  0.0060283 -6.4525 1.123e-10 ***
## lhinc           -0.0089439  0.0039088 -2.2882  0.022138 *
## factor(period)2 -0.0042799  0.0040162 -1.0657  0.286590
## factor(period)3 -0.0108953  0.0040151 -2.7136  0.006662 **
## factor(period)4 -0.0123002  0.0040174 -3.0618  0.002203 **
## factor(period)5 -0.0176797  0.0040152 -4.4032 1.072e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares:    1036.8
## Residual Sum of Squares: 1033.6
## R-Squared:      0.0031198
## F-statistic: 11.8119 on 6 and 22646 DF, p-value: 2.8458e-13
summary(Probit_Pooled <- glm(lfp ~ kids + lhinc + educ + black + age + agesq + factor(period), data=lfpP, family=binomial(link="probit")))
##
## Call:
## glm(formula = lfp ~ kids + lhinc + educ + black + age + agesq +
##     factor(period), family = binomial(link = "probit"), data = lfpP)
##
## Deviance Residuals:
##     Min       1Q   Median       3Q      Max
## -3.0040  -1.2284   0.6943   0.8688   2.7039
##
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)
## (Intercept)     -1.064e+00  1.373e-01  -7.751 9.13e-15 ***
## kids            -1.989e-01  7.549e-03 -26.348  < 2e-16 ***
## lhinc           -2.111e-01  1.330e-02 -15.870  < 2e-16 ***
## educ             7.969e-02  3.218e-03  24.763  < 2e-16 ***
## black            2.209e-01  3.348e-02   6.600 4.11e-11 ***
## age              1.449e-01  6.158e-03  23.533  < 2e-16 ***
## agesq           -1.991e-03  7.576e-05 -26.284  < 2e-16 ***
## factor(period)2 -1.242e-02  2.535e-02  -0.490   0.6241
## factor(period)3 -3.252e-02  2.530e-02  -1.285   0.1987
## factor(period)4 -4.610e-02  2.530e-02  -1.822   0.0684 .
## factor(period)5 -5.778e-02  2.525e-02  -2.288   0.0221 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
##     Null deviance: 35418  on 28314  degrees of freedom
## Residual deviance: 33113  on 28304  degrees of freedom
## AIC: 33135
##
## Number of Fisher Scoring iterations: 4
probitmfx(Probit_Pooled, data=lfp, atmean = F)
## Call:
## probitmfx(formula = Probit_Pooled, data = lfp, atmean = F)
##
## Marginal Effects:
##                       dF/dx   Std. Err.        z     P>|z|
## kids            -6.6018e-02  2.4187e-03 -27.2948 < 2.2e-16 ***
## lhinc           -7.0054e-02  4.3617e-03 -16.0611 < 2.2e-16 ***
## educ             2.6447e-02  1.0359e-03  25.5304 < 2.2e-16 ***
## black            6.9883e-02  1.0013e-02   6.9792 2.969e-12 ***
## age              4.8097e-02  1.9857e-03  24.2220 < 2.2e-16 ***
## agesq           -6.6085e-04  2.4247e-05 -27.2544 < 2.2e-16 ***
## factor(period)2 -4.1304e-03  8.4420e-03  -0.4893   0.62466
## factor(period)3 -1.0839e-02  8.4694e-03  -1.2798   0.20062
## factor(period)4 -1.5392e-02  8.4963e-03  -1.8116   0.07004 .
## factor(period)5 -1.9322e-02  8.5061e-03  -2.2716   0.02311 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## dF/dx is for discrete change for the following variables:
##
## [1] "black"           "factor(period)2" "factor(period)3" "factor(period)4"
## [5] "factor(period)5"
#library(dplyr)
lfp <- lfp %>% group_by(id) %>%
mutate(lfp1 = lfp[1], lfp_1 = lag(lfp, n = 1),
kidsbar=mean(kids), lhincbar=mean(lhinc))
#lfp <- ddply(lfp, ~id, lfp1=lfp[1], ... ) #if using plyr, instead of dplyr
lfpP <- pdata.frame(lfp, index = c("id", "period"))
summary(RE_Probit <- glm(lfp ~ kids + lhinc + kidsbar + lhincbar + educ + black + age + agesq + factor(period), data=lfpP, family=binomial("probit")))
##
## Call:
## glm(formula = lfp ~ kids + lhinc + kidsbar + lhincbar + educ +
##     black + age + agesq + factor(period), family = binomial("probit"),
##     data = lfpP)
##
## Deviance Residuals:
##     Min       1Q   Median       3Q      Max
## -2.3160  -1.2207   0.6885   0.8683   2.6359
##
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)
## (Intercept)     -7.261e-01  1.428e-01  -5.084 3.70e-07 ***
## kids            -1.174e-01  3.773e-02  -3.111  0.00187 **
## lhinc           -2.881e-02  2.520e-02  -1.143  0.25293
## kidsbar         -8.569e-02  3.849e-02  -2.227  0.02597 *
## lhincbar        -2.502e-01  2.948e-02  -8.486  < 2e-16 ***
## educ             8.413e-02  3.268e-03  25.742  < 2e-16 ***
## black            2.031e-01  3.357e-02   6.049 1.45e-09 ***
## age              1.516e-01  6.216e-03  24.397  < 2e-16 ***
## agesq           -2.067e-03  7.641e-05 -27.053  < 2e-16 ***
## factor(period)2 -1.357e-02  2.538e-02  -0.535  0.59291
## factor(period)3 -3.320e-02  2.533e-02  -1.311  0.18997
## factor(period)4 -3.903e-02  2.534e-02  -1.540  0.12344
## factor(period)5 -5.524e-02  2.529e-02  -2.185  0.02890 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
##     Null deviance: 35418  on 28314  degrees of freedom
## Residual deviance: 33033  on 28302  degrees of freedom
## AIC: 33059
##
## Number of Fisher Scoring iterations: 4
RE_APE<- probitmfx(RE_Probit, data=lfpP, atmean = F)
RE_APE
## Call:
## probitmfx(formula = RE_Probit, data = lfpP, atmean = F)
##
## Marginal Effects:
##                       dF/dx   Std. Err.        z     P>|z|
## kids            -3.8852e-02  1.2484e-02  -3.1120  0.001858 **
## lhinc           -9.5363e-03  8.3409e-03  -1.1433  0.252902
## kidsbar         -2.8364e-02  1.2736e-02  -2.2272  0.025937 *
## lhincbar        -8.2811e-02  9.7246e-03  -8.5156 < 2.2e-16 ***
## educ             2.7849e-02  1.0467e-03  26.6061 < 2.2e-16 ***
## black            6.4344e-02  1.0116e-02   6.3607 2.009e-10 ***
## age              5.0195e-02  1.9947e-03  25.1642 < 2.2e-16 ***
## agesq           -6.8428e-04  2.4342e-05 -28.1113 < 2.2e-16 ***
## factor(period)2 -4.4999e-03  8.4318e-03  -0.5337  0.593564
## factor(period)3 -1.1038e-02  8.4571e-03  -1.3051  0.191850
## factor(period)4 -1.2987e-02  8.4719e-03  -1.5329  0.125302
## factor(period)5 -1.8420e-02  8.4891e-03  -2.1698  0.030022 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## dF/dx is for discrete change for the following variables:
##
## [1] "black"           "factor(period)2" "factor(period)3" "factor(period)4"
## [5] "factor(period)5"
#library(pglm)
summary(CRE_Probit <- pglm(lfp ~ kids + lhinc + kidsbar + lhincbar + educ + black + age + agesq + factor(period), data=lfpP, model="random", effect = "individual", family=binomial("probit")))  
## --------------------------------------------
## Maximum Likelihood estimation
## Newton-Raphson maximisation, 7 iterations
## Return code 2: successive function values within tolerance limit
## Log-Likelihood: -8717.271
## 14  free parameters
## Estimates:
##                   Estimate Std. error t value  Pr(> t)
## (Intercept)     -3.2619293  0.7839086  -4.161 3.17e-05 ***
## kids            -0.3705209  0.0670346  -5.527 3.25e-08 ***
## lhinc           -0.0969143  0.0452282  -2.143 0.032130 *
## kidsbar         -0.3280459  0.0863456  -3.799 0.000145 ***
## lhincbar        -0.8722061  0.0817363 -10.671  < 2e-16 ***
## educ             0.2931048  0.0212183  13.814  < 2e-16 ***
## black            1.0279369  0.1612094   6.376 1.81e-10 ***
## age              0.5564136  0.0381874  14.571  < 2e-16 ***
## agesq           -0.0074894  0.0004929 -15.193  < 2e-16 ***
## factor(period)2 -0.0362766  0.0477579  -0.760 0.447498
## factor(period)3 -0.1119849  0.0478704  -2.339 0.019318 *
## factor(period)4 -0.1211022  0.0478184  -2.533 0.011324 *
## factor(period)5 -0.1815886  0.0473095  -3.838 0.000124 ***
## sigma            4.4862132  0.0783991  57.223  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------
probitmfx(CRE_Probit, data=lfpP, atmean = F)
## Call:
## probitmfx(formula = CRE_Probit, data = lfpP, atmean = F)
##
## Marginal Effects:
##                       dF/dx   Std. Err.        z     P>|z|
## kids            -3.8852e-02  1.2484e-02  -3.1120  0.001858 **
## lhinc           -9.5363e-03  8.3409e-03  -1.1433  0.252902
## kidsbar         -2.8364e-02  1.2736e-02  -2.2272  0.025937 *
## lhincbar        -8.2811e-02  9.7246e-03  -8.5156 < 2.2e-16 ***
## educ             2.7849e-02  1.0467e-03  26.6061 < 2.2e-16 ***
## black            6.4344e-02  1.0116e-02   6.3607 2.009e-10 ***
## age              5.0195e-02  1.9947e-03  25.1642 < 2.2e-16 ***
## agesq           -6.8428e-04  2.4342e-05 -28.1113 < 2.2e-16 ***
## factor(period)2 -4.4999e-03  8.4318e-03  -0.5337  0.593564
## factor(period)3 -1.1038e-02  8.4571e-03  -1.3051  0.191850
## factor(period)4 -1.2987e-02  8.4719e-03  -1.5329  0.125302
## factor(period)5 -1.8420e-02  8.4891e-03  -2.1698  0.030022 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## dF/dx is for discrete change for the following variables:
##
## [1] "black"           "factor(period)2" "factor(period)3" "factor(period)4"
## [5] "factor(period)5"
summary(FE_Logit <- pglm(lfp ~ kids + lhinc, data=lfpP, effect = "time", family=binomial("logit")))
## --------------------------------------------
## Maximum Likelihood estimation
## Newton-Raphson maximisation, 7 iterations
## Return code 2: successive function values within tolerance limit
## Log-Likelihood: -9025.528
## 4  free parameters
## Estimates:
##             Estimate Std. error t value  Pr(> t)
## (Intercept)  5.63516    0.56080  10.048  < 2e-16 ***
## kids        -0.66352    0.04743 -13.989  < 2e-16 ***
## lhinc       -0.24869    0.07131  -3.488 0.000487 ***
## sigma        8.66079    0.19558  44.283  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------

Note: Something is wrong here. CRE_Probit and FE_Logit are not exactly the same as in the textbook or the Stata version here..

#### Example 15.6

Dynamic Women’s LFP Equation

lfp <- lfp %>% group_by(id) %>%
mutate(kids2 = kids[2], kids3=kids[3], kids4 = kids[4], kids5 = kids[5],
lhinc2 = lhinc[2], lhinc3=lhinc[3], lhinc4 = lhinc[4], lhinc5 = lhinc[5])
lfpP <- pdata.frame(lfp, index = c("id", "period"))

# summary(RE_Probit <- pglm(lfp ~ lfp1 + lfp_1 + kids + lhinc + kidsbar + lhincbar + kids2 + kids3 + kids4 + kids5 + lhinc2+ lhinc3 + lhinc4 + lhinc5 + educ + black + age + agesq + per3 + per4 + per5, data=lfpP, model="random", family=binomial(link="probit")))
# probitmfx(RE_Probit, data=lfpP, atmean = F)

summary(Probit_Pooled <- glm(lfp ~ lfp_1 + kids + lhinc + kidsbar + lhincbar + educ + black + age + agesq, data=lfp, family=binomial(link="probit")))
##
## Call:
## glm(formula = lfp ~ lfp_1 + kids + lhinc + kidsbar + lhincbar +
##     educ + black + age + agesq, family = binomial(link = "probit"),
##     data = lfp)
##
## Deviance Residuals:
##     Min       1Q   Median       3Q      Max
## -2.8325  -0.4037   0.2895   0.3316   2.6913
##
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.0808656  0.2326162  -8.945  < 2e-16 ***
## lfp_1        2.8732966  0.0269968 106.431  < 2e-16 ***
## kids        -0.0408845  0.0629982  -0.649    0.516
## lhinc       -0.0596502  0.0404712  -1.474    0.141
## kidsbar     -0.0209056  0.0641473  -0.326    0.745
## lhincbar    -0.0763569  0.0473610  -1.612    0.107
## educ         0.0306517  0.0053277   5.753 8.75e-09 ***
## black        0.0741814  0.0545141   1.361    0.174
## age          0.0866667  0.0101908   8.504  < 2e-16 ***
## agesq       -0.0011241  0.0001253  -8.969  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
##     Null deviance: 28409  on 22651  degrees of freedom
## Residual deviance: 10664  on 22642  degrees of freedom
##   (5663 observations deleted due to missingness)
## AIC: 10684
##
## Number of Fisher Scoring iterations: 5
probitmfx(Probit_Pooled, data=lfp, atmean = F)
## Call:
## probitmfx(formula = Probit_Pooled, data = lfp, atmean = F)
##
## Marginal Effects:
##                dF/dx   Std. Err.        z     P>|z|
## lfp_1     8.3684e-01  4.4786e-03 186.8545 < 2.2e-16 ***
## kids     -5.0871e-03  7.8389e-03  -0.6490    0.5164
## lhinc    -7.4221e-03  5.0366e-03  -1.4736    0.1406
## kidsbar  -2.6012e-03  7.9817e-03  -0.3259    0.7445
## lhincbar -9.5008e-03  5.8944e-03  -1.6118    0.1070
## educ      3.8139e-03  6.6489e-04   5.7361 9.689e-09 ***
## black     9.1341e-03  6.6477e-03   1.3740    0.1694
## age       1.0784e-02  1.2761e-03   8.4504 < 2.2e-16 ***
## agesq    -1.3987e-04  1.5706e-05  -8.9055 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## dF/dx is for discrete change for the following variables:
##
## [1] "lfp_1" "black"