I am trying to calculate the average marginal effects for the Chamberlain-Mundlak Correlated Random Effects probit model. The ultimate goal is to get something equivalent to the AME from the fixed effects panel logit. The problem with the latter is that it eliminates the FE before estimation, making it impossible to include them in the AME.
As I understand it, one possible way to do this is to fit a panel random effects probit where the RHS variables are augmented with $\bar x_i$, the average of $x_{it}$ for each panel. This is the Chamberlain-Mundlak CRE. The inclusion of the mean terms should capture the correlation between the unobserved heterogeneity and the covariates that renders the random effect model inconsistent.
If I understood Jeff Wooldridge's notes (p. 40), the AME would be the derivative of RE probit model's
$$\frac{1}{N} \sum_{i=1}^N \Phi(x_{it}\hat \beta_a + \hat \psi_a + \bar x_i\hat\xi_a)$$
with respect to $x$, and $\hat \beta_a=\frac{\hat \beta}{\sqrt(1+\hat \sigma^2)}$.
That would mean that
$$AME =\frac{1}{N} \sum_{i=1}^N \phi(x_{it}\hat \beta_a + \hat \psi_a + \bar x_i\hat\xi_a) \cdot \frac{\beta_a}{\sqrt(1+\hat \sigma^2)}$$
This means that I just need to scale the ordinary average marginal effects from the RE panel probit like JW does on page 52 with the index function coefficients.
However, the estimates I get are pretty different than the other equivalent methods he suggests that do not require rescaling, like the pooled probit, xtgee, and GLM, or even FE OLS (all of which are very similar). I am struggling to figure out why that is the case.
Here's my Stata code:
set more off
estimates drop _all
copy https://mitpress.mit.edu/sites/default/files/titles/content/wooldridge/statafiles.zip statafiles.zip, replace
unzipfile statafiles.zip, replace
use "lfp.dta", clear
xtset id period // balanced panel
/* Calculate means of covariates for each id */
foreach var of varlist kids lhinc {
egen `var'bar = mean(`var'), by(id)
}
/* (1) FE Panel Logit */
qui xtlogit lfp kids lhinc i.period, fe nolog
margins, dydx(kids lhinc) post
estimates store m1, title("Panel Logit FE") // these set the FE to its average of zero
/* (2) FE Linear Panel Model */
qui xtreg lfp kids lhinc i.period, fe cluster(id)
margins, dydx(kids lhinc) post
estimates store m2, title("Panel FE")
/* (3) Pooled Probit (Inefficient Relative to CRE) */
probit lfp kids lhinc kidsbar lhincbar educ i.black c.age##c.age i.period, cluster(id) nolog
margins, dydx(kids lhinc) post
estimates store m3, title("Pooled Probit")
/* (4) GEE may be more efficient */
qui xtgee lfp kids lhinc kidsbar lhincbar educ i.black c.age##c.age i.period, ///
fam(bin) link(probit) corr(exch) robust nolog
margins, dydx(kids lhinc) post
estimates store m4, title("Panel GEE")
/* (5) GLM */
qui glm lfp kids lhinc kidsbar lhincbar educ i.black c.age##c.age i.period, ///
fam(bin) link(probit) cluster(id) nolog
margins, dydx(kids lhinc) post
estimates store m5, title("GLM")
/* (6) Chamberlain-Mundlak Device Style Correlated Random Effects Probit */
/* ("Style" means also adding other time-invariant variables) */
qui xtprobit lfp kids lhinc kidsbar lhincbar educ i.black c.age##c.age i.period, re nolog
/* (a) Scaled coefficients to compare with pooled index function probit ones in (3) */
di (1/sqrt(1 + e(sigma_u)^2))*_b[kids]
di (1/sqrt(1 + e(sigma_u)^2))*_b[lhinc]
local factor = 1/sqrt(1 + e(sigma_u)^2) // scale factor
/* (b) CRE MEs Attempt #1 */
margins, expression(normalden(xb())*(_b[kids])*(1/sqrt(1 + e(sigma_u)^2))) force
margins, expression(normalden(xb())*(_b[lhinc])*(1/sqrt(1 + e(sigma_u)^2))) force
/* (c) CRE MEs Attempt #2 */
margins, dydx(kids lhinc) post
estimates store m6, title("CRE")
/* (d) CRE MEs Attempt #3 */
nlcom (kids:_b[kids]*`factor') (lhinc:_b[lhinc]*`factor')
esttab *, se mtitle modelwidth(15)
This yields:
esttab *, se mtitle modelwidth(15)
------------------------------------------------------------------------------------------------------------------------------
(1) (2) (3) (4) (5) (6)
Panel Logit FE Panel FE Pooled Probit Panel GEE GLM CRE
------------------------------------------------------------------------------------------------------------------------------
kids -0.0507* -0.0389*** -0.0389*** -0.0373*** -0.0389*** -0.0302***
(0.0254) (0.00917) (0.00892) (0.00932) (0.00892) (0.00533)
lhinc -0.0145*** -0.00894 -0.00954* -0.00917 -0.00954* -0.00762*
(0.00150) (0.00459) (0.00475) (0.00491) (0.00475) (0.00357)
------------------------------------------------------------------------------------------------------------------------------
N 5275 28315 28315 28315 28315 28315
------------------------------------------------------------------------------------------------------------------------------
Standard errors in parentheses
* p<0.05, ** p<0.01, *** p<0.001
(1) is very different since the FE logit ignores the FE. (2)-(5) are all very similar. The non-scaled CRE AMEs are shown in (6).
I tried 3 ways of scaling these (all which seem more off than the unscaled version):
. qui xtprobit lfp kids lhinc kidsbar lhincbar educ i.black c.age##c.age i.period, re nolog
.
. /* (a) Scaled coefficients to compare with pooled index function probit ones in (3) */
. di (1/sqrt(1 + e(sigma_u)^2))*_b[kids]
-.08865639
. di (1/sqrt(1 + e(sigma_u)^2))*_b[lhinc]
-.02240692
. local factor = 1/sqrt(1 + e(sigma_u)^2) // scale factor
.
. /* (b) CRE MEs Attempt #1 */
. margins, expression(normalden(xb())*(_b[kids])*(1/sqrt(1 + e(sigma_u)^2))) force
(note: expression is a function of possibly stochastic quantities other than e(b))
Predictive margins Number of obs = 28,315
Model VCE : OIM
Expression : normalden(xb())*(_b[kids])*(1/sqrt(1 + e(sigma_u)^2))
------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_cons | -.0061709 .0011396 -5.41 0.000 -.0084045 -.0039372
------------------------------------------------------------------------------
. margins, expression(normalden(xb())*(_b[lhinc])*(1/sqrt(1 + e(sigma_u)^2))) force
(note: expression is a function of possibly stochastic quantities other than e(b))
Predictive margins Number of obs = 28,315
Model VCE : OIM
Expression : normalden(xb())*(_b[lhinc])*(1/sqrt(1 + e(sigma_u)^2))
------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_cons | -.0015596 .0007353 -2.12 0.034 -.0030008 -.0001185
------------------------------------------------------------------------------
.
. /* (c) CRE MEs Attempt #2 */
. margins, dydx(kids lhinc) post
Average marginal effects Number of obs = 28,315
Model VCE : OIM
Expression : Pr(lfp=1), predict(pr)
dy/dx w.r.t. : kids lhinc
------------------------------------------------------------------------------
| Delta-method
| dy/dx Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
kids | -.0301539 .0053261 -5.66 0.000 -.0405928 -.0197149
lhinc | -.0076211 .0035696 -2.13 0.033 -.0146174 -.0006247
------------------------------------------------------------------------------
. estimates store m6, title("CRE")
.
. /* (d) CRE MEs Attempt #3 */
. nlcom (kids:_b[kids]*`factor') (lhinc:_b[lhinc]*`factor')
kids: _b[kids]*.2233101065698412
lhinc: _b[lhinc]*.2233101065698412
------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
kids | -.0067337 .0011894 -5.66 0.000 -.0090648 -.0044025
lhinc | -.0017019 .0007971 -2.13 0.033 -.0032642 -.0001395
------------------------------------------------------------------------------
