Average Marginal Effects from Chamberlain-Mundlak Device CRE

[Dimitriy V. Masterov]

I am trying to calculate the average marginal effects for the Chamberlain-Mundlak Correlated Random Effects probit model.

I can do this in several ways (see below). The one that is giving me trouble is the panel random effects probit where the RHS variables are augmented with $\bar x_i$, the average of $x_{it}$ for each panel. I thought the AMEs from this model, appropriately rescaled, would agree with all the other equivalent ways of doing it, but they do not. If anything, the unscaled version seems best.
I would like to get some intuition on why this is not the case (or what I am doing wrong).

Moreover, my scaled index function coefficients in (a) (-.08865639 and -.02240692) do not agree with the pooled probit ones in (3) (-.1173749 and -.0288098).

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)

The relevant part of the results is:

. 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
------------------------------------------------------------------------------

. 
. 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

I would have expected (1) to be off since it sets the FE to zero. I would like to figure out how to match (6) with (2)-(5).

[Dimitriy V. Masterov]

I contacted Stata Technical Support, who informed me that this lies outside the scope. I would really like to know the answer to this.

[Dimitriy V. Masterov]

Cross-posted here.

[Sourabh Paul]

Hi Dimitriy,
I was also replicating JMW's slides on Mundlak's CRE models and found your post here. You need to scale the xb() by multiplying 1/sqrt(1 + e(sigma_u)^2) as well. You only scaled the coefficient and not the density function.
* Try this:
margins, expression(normalden(xb()*(1/sqrt(1 + e(sigma_u)^2)))*(_b[kids])*(1/sqrt(1 + e(sigma_u)^2))) force

*Also, as described in JMW, you can find it manually like this:
predict xbhat, xb
gen xbhata = xbhat/sqrt(1 + e(sigma_u)^2)
gen scale = normalden(xbhata)
qui su scale
di r(mean)*_b[kids]/sqrt(1 + e(sigma_u)^2)
di r(mean)*_b[lhinc]/sqrt(1 + e(sigma_u)^2)

[Dimitriy V. Masterov]

I think you are right on the derivation, but I am still not quite matching the other specifications. If anything, it's pretty close to the AMEs without any rescaling.

I am getting:

. /* (b) CRE MEs Attempt #1 */
. margins, expression(normalden(xb()*(1/sqrt(1 + e(sigma_u)^2)))*(_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()*(1/sqrt(1 + e(sigma_u)^2)))*(_b[kids])*(1/sqrt(1 + e(sigma_u)^2))

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _cons |  -.0288166   .0050865    -5.67   0.000    -.0387861   -.0188472
------------------------------------------------------------------------------

. margins, expression(normalden(xb()*(1/sqrt(1 + e(sigma_u)^2)))*(_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()*(1/sqrt(1 + e(sigma_u)^2)))*(_b[lhinc])*(1/sqrt(1 + e(sigma_u)^2))

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _cons |  -.0072831   .0034109    -2.14   0.033    -.0139683   -.0005979
------------------------------------------------------------------------------

[Sourabh Paul]

Yes, you are right. The xtprobit (mundlak) result itself does not match with JMW (pp 51 of the slides). For example, the coefficient of kids is -0.3174 whereas my estimate is -0.397. Also, my sigma_u is 4.365 and that of JMW is 2.383. I thought it could be due to the version (13) I am using (the default routine for estimation may be different). I tried to use the following option with xtprobit and it is much closer to Wooldidge:

xtprobit lfp kids lhinc kidsbar lhincbar educ i.black c.age##c.age i.period, re intm(agh) intp(12)

Please post here if can successfully replicate JMW. Thanks.

[Jeff Wooldridge]

Hi Dimitriy, Sourabh: I suspect I didn't have something quite right in those lecture notes. I've learned some things since then, and when I have a little time I'll look for some more current code. One issue is that it used to be easiest to insert the mean value of the heterogeneity but then average across the covariates. That's not the same as the average partial effects computed in other ways. I need to do a little reviewing.

It is true that the estimation of the random effects probit changed across Stata versions. In fact, the reported estimates in Chapter 15 of my book are no longer what one gets in Stata 14.

Jeff

[Sourabh Paul]

Hi Jeff: Thank you so much for your comments. Looking forward to your revised notes.

May I also request you to consider incorporating a section on survey weights in Nonlinear Panel models? Whether this could be implemented using generalized linear latent and mixed models (gllamm) and whether Average Partial Effects can be estimated using the approach mentioned in your text. One way out is to use Pooled Probit model with weights and directly derive the marginal effects. Will it be robust? How efficient will it be?

Thanks a lot!

Sourabh

[Dimitriy V. Masterov]

Thanks, Prof. Wooldridge! I tried running the code under version control, and this did not give me your estimates.

[Sasha Gulabivala]

Hi,

Thank you for the insights above.

I am also working through Prof Wooldridge's LFP example before I estimate the APEs for my own research project, which will use the CRE Probit, MLE estimation.

In the textbook Econometric Analysis of Cross Section and Panel Data (Chapter 15, Page 485) , Wooldridge writes that:

We only need to estimate betac = beta / ( 1 + sigmac^2)^1/2 to estimate the APEs, for either continuous or discrete explanatory variables

(1) I wondered why then, to obtain the APEs, it is suggested by Sourabh Paul above that the density function should be scaled (I understand that the coefficient should be scaled):

You need to scale the xb() by multiplying 1/sqrt(1 + e(sigma_u)^2) as well. You only scaled the coefficient and not the density function.
* Try this:
margins, expression(normalden(xb()*(1/sqrt(1 + e(sigma_u)^2)))*(_b[kids])*(1/sqrt(1 + e(sigma_u)^2))) force

(2) Also, why has the square root of (1/sqrt(1 + e(sigma_u)^2)) not been taken, in accordance with the textbook quote above?


Thanks in advance

[Sasha Gulabivala]

Regarding (2), sorry, please ignore the second question - I did not see the "sqrt" initially. I now understand that it has been included.
However, regarding (1), I am still puzzled as to why the density function needs to be scaled? Thanks.

[Sasha Gulabivala]

Jeff Wooldridge please would you be able to suggest the Stata code for computing AMEs for a CRE Probit model (xtprobit, re)?

Do you think

margins, expression(normalden(xb()*(1/sqrt(1 + e(sigma_u)^2)))*(_b[kids])*(1/sqrt(1 + e(sigma_u)^2))) force

should work?

Many thanks

[Rose Simmons]

Hello Dimitriy V. Masterov , Sourabh Paul , Jeff Wooldridge , Saeesha Gulabivala

Having read the Wooldridge textbook (2002), it looks like the correlated random effects model would be interesting and good for my panel data (where the y variable is discrete, either 1 or 0). I did not want to use FE as it eliminates my time-invariant variables. I also did not want to use Logit RE as I wanted to calculate meaningful marginal effects. So I am considering using Probit CRE.

I wondered whether there are any Stata tests I would need to conduct before running a CRE probit model, so that I know whether I can use it for my data?

Any advice would be much appreciated

[Sasha Gulabivala]

Jeff Wooldridge , the examples in your textbook use kids and lhinc, which are both continuous, and I think the code suggested by Sourabh Paul for calculating AMEs would work for continuous variables:

  
 margins, expression(normalden(xb()*(1/sqrt(1 + e(sigma_u)^2)))*(_b[kids])*(1/sqrt(1 + e(sigma_u)^2))) force

However, would you be able to suggest the Stata command for discrete variables please?

[Yulia Shenderovich]

Hello Dimitriy V. Masterov, Jeff Wooldridge , Sasha Gulabivala , Sourabh Paul , Rose Simmons

Thank you for this very useful discussion! I have been reading this thread, Prof. Wooldridge's textbook, and also this recent paper 10.1080/13504851.2018.1441498 (James R. Bland & Amanda C. Cook (2019) Random effects probit and logit: understanding predictions and marginal effects, Applied Economics Letters, 26:2, 116-123) to understand how to compute marginal effects in random effects probit models. Based on this thread, I had three key follow-up questions.

Question 1. Calculating marginal effects for a binary predictor (as asked above). Bland & Cook provide the following code: margins, dydx(x1) expression(normal(xb()/sqrt(1+ exp(_b[lnsig2u:_cons])))). When I try this with the lfp.dta dataset cited here, this gives very similar results as does a simple dydx(x1), as you can see in the output below.

I run the regression similar as the posts above but use “de-meaned” Xit observations to generate the within- and between-person effects (as per Bell, A., Fairbrother, M., & Jones, K. (2019). Fixed and random effects models: making an informed choice. Quality & Quantity, 53(2), 1051-1074. https://doi.org/10.1007/s11135-018-0802-x)

Using Stata/SE 14.2 and the LFP.dta data example available here on the MIT Press website (http://mitp-content-server.mit.edu:1...tata_files.zip), I produce the following:

 
. /*1. Create a binary predictor that changes within person over time, for illustration - husband's income */
. /*The cut-off is at the mean - below mean is "low income" and above "high", just for illustration */
. gen lhinc_bin = lhinc
. recode lhinc_bin (min/7.80503 = 0) (7.80504/max = 1)
(lhinc_bin: 28315 changes made)
 
. tab lhinc_bin
 
  lhinc_bin |      Freq.     Percent        Cum.
------------+-----------------------------------
          0 |     13,344       47.13       47.13
          1 |     14,971       52.87      100.00
------------+-----------------------------------
      Total |     28,315      100.00
 
. /* 2. Generate person-specific means and deviations for time-varying predictors*/
. foreach var of varlist kids lhinc_bin {
  2. egen `var'bar = mean(`var'), by(id) //generate mean
  3. gen d_`var'=`var'-`var'bar //gen deviations from means
  4. }
 
. /* 3. Using the regression model from above to explain labor force participation, a binary outcome*/
. xtset id period
       panel variable:  id (strongly balanced)
        time variable:  period, 1 to 5
                delta:  1 unit
 
. qui xtprobit lfp d_kids d_lhinc_bin kidsbar lhinc_binbar educ i.black c.age##c.age i.period, re nolog
.
. /* 4a. Approach 1 to compute average marginal effects for a binary predictor */
. margins, dydx (d_lhinc_bin)
 
Average marginal effects                        Number of obs     =     28,315
Model VCE    : OIM
 
Expression   : Pr(lfp=1), predict(pr)
dy/dx w.r.t. : d_lhinc_bin
 
------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
 d_lhinc_bin |   -.006587   .0044022    -1.50   0.1346    -.0152151    .0020411
------------------------------------------------------------------------------
 
.
. /* 4b. Approach 2 to compute average marginal effects for a binary predictor, based on Bland and Cook */
. margins, dydx (d_lhinc_bin) expression(normal(xb()/sqrt(1+ exp(_b[lnsig2u:_cons]))))
 
Average marginal effects                        Number of obs     =     28,315
Model VCE    : OIM
 
Expression   : normal(xb()/sqrt(1+ exp(_b[lnsig2u:_cons])))
dy/dx w.r.t. : d_lhinc_bin
 
------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
 d_lhinc_bin |  -.0063035   .0042127    -1.50   0.1346    -.0145603    .0019532
------------------------------------------------------------------------------

Which of these two (or another) way of calculating marginal effects for a binary predictor in a probit model with random effects would you recommend?

Question 2. I was also curious about interpreting these effects in prose. My understanding is that we would interpret the partial effect of d_lhnic_bin (the within-person effect) to indicate that participants had 0.006 lower probability of being in the workforce when their husband’s income became “high”, at the sample grand mean of all the other variables (acknowledging that this result is non-significant). What interpretation would you recommend?

Question 3. Furthermore, do you have any advice on the interpretation for marginal effects of the person-average predictor, e.g. lhinc_binbar? Since this represents between-person difference, I was thinking this would be something like “0.090 average lower probability of being in the workforce for a person whose husband’s income was high”

 

. margins, dydx (lhinc_binbar)
 
Average marginal effects                        Number of obs     =     28,315
Model VCE    : OIM
 
Expression   : Pr(lfp=1), predict(pr)
dy/dx w.r.t. : lhinc_binbar
 
------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
lhinc_binbar |  -.0899362   .0105402    -8.53   0.0000    -.1105945   -.0692779
------------------------------------------------------------------------------

Thank you very much!
Yulia