## Multiple equation models: Estimation and marginal effects using mlexp

We continue with the series of posts where we illustrate how to obtain correct standard errors and marginal effects for models with multiple steps. In this post, we estimate the marginal effects and standard errors for a hurdle model with two hurdles and a lognormal outcome using **mlexp**. **mlexp** allows us to estimate parameters for multiequation models using maximum likelihood. In the last post (Multiple equation models: Estimation and marginal effects using gsem), we used **gsem** to estimate marginal effects and standard errors for a hurdle model with two hurdles and an exponential mean outcome.

We exploit the fact that the hurdle-model likelihood is separable and the joint log likelihood is the sum of the individual hurdle and outcome log likelihoods. We estimate the parameters of each hurdle and the outcome separately to get initial values. Then, we use **mlexp** to estimate the parameters of the model and **margins** to obtain marginal effects.

**Starting points: model and initial values**

We model the amount spent on dental care. The first hurdle is whether an individual spends or does not spend on dental care. The second hurdle is the individual level of insurance coverage. Different levels of coverage lead to different spending on dental care.

We assume probit and ordered probit models for the two hurdles. In contrast to the previous post, we use a lognormal distribution to model the amount spent. With these distributional assumptions, we use maximum likelihood rather than quasi-likelihood, as in the previous post, to estimate the parameters of the model.

We obtain initial values by **mlexp** for each hurdle and the lognormal outcome and store the log-likelihood expression for each step in a local macro. These local macros are summed together in the final use of **mlexp** to get parameter estimates and standard errors.

**spend** is a binary outcome for whether an individual spends money on dental care. This is the first hurdle variable. We store the log-likelihood expression in the local macro **probit**. See Appendix for more information. Then, we use **mlexp** to estimate the parameters of the hurdle. The point estimates are stored in the matrix **binit**.

. local probit ln(cond(1.spend,normal({spend: x1 x2 x4 _cons}), > 1-normal({spend:}))) . mlexp (`probit') initial: log likelihood = -6931.4718 alternative: log likelihood = -5926.3721 rescale: log likelihood = -5926.3721 Iteration 0: log likelihood = -5926.3721 Iteration 1: log likelihood = -4789.4607 Iteration 2: log likelihood = -4776.9361 Iteration 3: log likelihood = -4776.9332 Iteration 4: log likelihood = -4776.9332 Maximum likelihood estimation Log likelihood = -4776.9332 Number of obs = 10,000 ---------------------------------------------------------------------------- | Coef. Std. Err. z P>|z| [95% Conf. Interval] -----------+---------------------------------------------------------------- spend | x1 | .5138169 .0160374 32.04 0.000 .4823841 .5452498 x2 | -.5276993 .0224236 -23.53 0.000 -.5716488 -.4837498 x4 | .4822064 .0185374 26.01 0.000 .4458736 .5185391 _cons | .5568866 .02899 19.21 0.000 .5000672 .613706 ---------------------------------------------------------------------------- . matrix binit = e(b)

**insurance** is a three-level ordered outcome indicating insurance level. This is the variable for the second hurdle. We store the log-likelihood expression in the local macro **oprobit** and use **mlexp** as before.

. local oprobit ln(cond(insurance==1,normal(-{insurance: x3 x4}+{cut1}), > cond(insurance==2,normal({cut2}-{insurance:})- > normal({cut1}-{insurance:}), > 1-normal({cut2}-{insurance:})))) . mlexp (`oprobit') initial: log likelihood = -(could not be evaluated) feasible: log likelihood = -23924.936 rescale: log likelihood = -19788.939 rescale eq: log likelihood = -11884.962 Iteration 0: log likelihood = -11884.962 Iteration 1: log likelihood = -10261.611 Iteration 2: log likelihood = -10227.115 Iteration 3: log likelihood = -10226.993 Iteration 4: log likelihood = -10226.993 Maximum likelihood estimation Log likelihood = -10226.993 Number of obs = 10,000 ---------------------------------------------------------------------------- | Coef. Std. Err. z P>|z| [95% Conf. Interval] -----------+---------------------------------------------------------------- insurance | x3 | .2926129 .0100288 29.18 0.000 .2729569 .312269 x4 | -.2754986 .0144941 -19.01 0.000 -.3039066 -.2470906 -----------+---------------------------------------------------------------- /cut1 | -1.270488 .0255784 -49.67 0.000 -1.320621 -1.220355 /cut2 | -.2612825 .0235687 -11.09 0.000 -.3074763 -.2150887 ---------------------------------------------------------------------------- . matrix binit = binit, e(b)

Now, we use **mlexp** to estimate the parameters of the lognormal outcome. **spent** corresponds to the amount spent on dental care. We use factor-variable notation to use a different intercept for each level of insurance, \(\beta_{ins,1}\ldots,\beta_{ins,3}\). The covariates are specified in equation **spent**, and the constant intercepts are specified in **spent_int**. The log-likelihood expression is saved in the local macro **lognormal**. We restrict estimation to the positive sample and use **mlexp** to estimate the outcome parameters.

. local lognormal -.5*((ln(spent)-{spent: x4 x5 x6} - > {spent_int: ibn.insurance})/ > (exp({lnsigma})))^2 > -ln((spent*exp({lnsigma})*sqrt(2*_pi))) . mlexp (`lognormal') if spend initial: log likelihood = -16596.787 alternative: log likelihood = -16544.473 rescale: log likelihood = -15515.652 rescale eq: log likelihood = -14206.308 Iteration 0: log likelihood = -14206.308 Iteration 1: log likelihood = -13818.45 Iteration 2: log likelihood = -13520.664 Iteration 3: log likelihood = -13519.085 Iteration 4: log likelihood = -13519.084 Maximum likelihood estimation Log likelihood = -13519.084 Number of obs = 7,228 ---------------------------------------------------------------------------- | Coef. Std. Err. z P>|z| [95% Conf. Interval] -----------+---------------------------------------------------------------- spent | x4 | .2953067 .0121253 24.35 0.000 .2715415 .3190718 x5 | -.2917276 .0088451 -32.98 0.000 -.3090637 -.2743915 x6 | .2891539 .0227041 12.74 0.000 .2446548 .3336531 -----------+---------------------------------------------------------------- spent_int | insurance | 1 | .1924971 .0239883 8.02 0.000 .1454809 .2395134 2 | .2950143 .0214191 13.77 0.000 .2530336 .3369951 3 | .454055 .0202381 22.44 0.000 .414389 .4937209 -----------+---------------------------------------------------------------- /lnsigma | -.2953913 .0083172 -35.52 0.000 -.3116926 -.2790899 ---------------------------------------------------------------------------- . matrix binit = binit, e(b)

**Joint Estimation and marginal effects**

Now, we use **mlexp** to estimate the parameters of the joint model. The joint log likelihood is specified as the sum of the individual log likelihoods. We merely add up the local macros that we created in the last section. The matrix **binit** contains the point estimates from the individual steps. We specify this matrix in **from()** to give **mlexp** good starting values. We specify **vce(robust)** so that we can use **margins** to estimate the marginal effects over the population of covariates using **margins** and **vce(unconditional)**.

mlexp (`probit' + `oprobit' + cond(spend,(`lognormal'),0)), > vce(robust) from(binit) Iteration 0: log pseudolikelihood = -28523.01 Iteration 1: log pseudolikelihood = -28523.01 Maximum likelihood estimation Log pseudolikelihood = -28523.01 Number of obs = 10,000 ---------------------------------------------------------------------------- | Robust | Coef. Std. Err. z P>|z| [95% Conf. Interval] -----------+---------------------------------------------------------------- spend | x1 | .513817 .0159257 32.26 0.000 .4826032 .5450307 x2 | -.5276994 .0224874 -23.47 0.000 -.5717739 -.4836249 x4 | .4822064 .0185804 25.95 0.000 .4457894 .5186234 _cons | .5568866 .0289611 19.23 0.000 .5001239 .6136494 -----------+---------------------------------------------------------------- insurance | x3 | .2926129 .0100343 29.16 0.000 .2729461 .3122798 x4 | -.2754986 .0144074 -19.12 0.000 -.3037366 -.2472605 -----------+---------------------------------------------------------------- cut1 | _cons | -1.270488 .0254071 -50.01 0.000 -1.320285 -1.220691 -----------+---------------------------------------------------------------- cut2 | _cons | -.2612825 .0236404 -11.05 0.000 -.3076167 -.2149482 -----------+---------------------------------------------------------------- spent | x4 | .2953067 .012 24.61 0.000 .271787 .3188263 x5 | -.2917276 .0088412 -33.00 0.000 -.309056 -.2743991 x6 | .2891539 .0224396 12.89 0.000 .2451731 .3331347 -----------+---------------------------------------------------------------- spent_int | insurance | 1 | .1924972 .0242748 7.93 0.000 .1449194 .240075 2 | .2950143 .021299 13.85 0.000 .253269 .3367597 3 | .4540549 .0202593 22.41 0.000 .4143473 .4937625 -----------+---------------------------------------------------------------- /lnsigma | -.2953912 .0084227 -35.07 0.000 -.3118993 -.278883 ----------------------------------------------------------------------------

The point estimates from the individual steps match the joint estimates. Now, we obtain marginal effects of the outcome covariates on the conditional mean of the outcome.

The conditional mean of expenditure on the independent variables is given by

\begin{eqnarray*}

E\left(\text{spent}|X\right) = \Phi(X_p{\boldsymbol \beta}_p)

\left[\begin{matrix}\Phi(\kappa_1-X_o{\boldsymbol \beta}_o)\exp{\beta}_{{ins},1} + \cr

\left\{\Phi(\kappa_2-X_o{\boldsymbol \beta}_o) – \Phi(\kappa_1-X_o{\boldsymbol \beta}_o)\right\} \exp{ \beta}_{{ins},2} + \cr \left\{1-\Phi(\kappa_2-X_o{\boldsymbol \beta}_o)\right\} \exp{ \beta}_{{ins},3}\end{matrix}\right] \exp\left(X_e{\boldsymbol \beta}_e + .5\sigma^2\right)\cr

\end{eqnarray*}

where \(\kappa_{1}\) and \(\kappa_{2}\) are the cutoff points from the insurance model, which represents insurance level. We use the subscripts p, o, and e to emphasize that the covariates and coefficients related to the probit, ordered probit, and exponential mean are different.

Now, we use **margins** to estimate the marginal effects. We use the **expression()** option to write an expression for the expected value of amount spent.

. margins, expression(normal(xb(spend))*( > exp(_b[spent_int:1.insurance])* > normal(_b[/cut1]-xb(insurance))+ > exp(_b[spent_int:2.insurance])* > (normal(_b[/cut2]-xb(insurance)) - > normal(_b[/cut1]-xb(insurance)))+ > exp(_b[spent_int:3.insurance])*( > 1-normal(_b[/cut2]-xb(insurance))))* > exp(xb(spent)+.5*exp(_b[/lnsigma])^2)) > dydx(x4 x5 x6) vce(unconditional) Average marginal effects Number of obs = 10,000 Expression : normal(xb(spend))*( exp(_b[spent_int:1.insurance])* normal(_b[/cut1]-xb(insurance))+ exp(_b[spent_int:2.insurance])* (normal(_b[/cut2]-xb(insurance)) - normal(_b[/cut1]-xb(insurance)))+ exp(_b[spent_int:3.insurance])*( 1-normal(_b[/cut2]-xb(insurance))))* exp(xb(spent)+.5*exp(_b[/lnsigma])^2) dy/dx w.r.t. : x4 x5 x6 ---------------------------------------------------------------------------- | Unconditional | dy/dx Std. Err. z P>|z| [95% Conf. Interval] -----------+---------------------------------------------------------------- x4 | .9525996 .0311351 30.60 0.000 .8915759 1.013623 x5 | -.6329119 .0224351 -28.21 0.000 -.6768838 -.58894 x6 | .6273282 .0499249 12.57 0.000 .5294772 .7251792 ----------------------------------------------------------------------------

**Final considerations**

We illustrated how to use **mlexp** to obtain the estimates and standard errors for a multiple hurdle model and its marginal effects. In subsequent posts, we obtain these results for other multistep models using other Stata tools.

**Appendix 1**

Below is the code used to produce the data.

clear set seed 11134 set obs 10000 // Generating exogenous variables generate x1 = rnormal() generate x2 = int(3*rbeta(2,3)) generate x3 = rchi2(1)-2 generate x4 = ln(rchi2(4)) generate x5 = rnormal() generate x6 = rbeta(2,3)>.6 // Generating unobservables generate ep = rnormal() // for probit generate eo = rnormal() // for ordered probit generate e = rnormal()*.75 // for lognormal equation // Generating linear predictions generate xbp = .5*(1 + x1 - x2 + x4) generate xbo = .3*(1 + x3 - x4) generate yotemp = xbo + eo generate insurance = yotemp replace insurance = 1 if yotemp < -1 replace insurance = 2 if yotemp> -1 & yotemp<0 replace insurance = 3 if yotemp> 0 generate xbe = .3*(- x5 + x6 + x4 + .5*insurance) // Generating outcomes generate spend = xbp + ep > 0 generate yexp = exp(xbe + e) generate spent = spend*yexp

**Appendix 2**

Macros are a programming feature provided by Stata. They can be used to make complicated expressions easy to write. Everywhere a punctuated macro name appears in a command, the macro contents are substituted for the macro name. So a complicated expression can be stored in a macro with a short name, and the expression can be used repeatedly by only typing the punctuated short name. See Programming an estimation command in Stata: Global macros versus local macros for a discussion of macros.

In this blog, we use local macros to store expressions used in calculating the log likelihood of our model. For the probit model, we stored the log-likelihood expression in the local macro **probit**.

. local probit > ln(cond(1.spend,normal({spend: x1 x2 x4 _cons}),1-normal({spend:})))

When we type the punctuated **`probit’** later, Stata uses the expression stored inside **probit**. We can see the expression by displaying it.

. display "`probit'" ln(cond(1.spend,normal({spend: x1 x2 x4 _cons}),1-normal({spend:})))