Estimation of marginal effects for a multivariate probit

[Malcolm Wardlaw]

Can anyone help me on calculating the marginal effects conditional marginal effects for a multivariate probit model? As an illustration, let's take the following example.

clear    
insheet using "http://people.stern.nyu.edu/wgreene/Econometrics/healthcare.csv", comma 
biprobit (doctor = age female educ) (hospital =  age female educ) 
margins, dydx(*) predict(pmarg1) atmeans
margins, dydx(*) predict(pcond1) atmeans

This reports the marginal effect of the covariates on doctor and the marginal effect of the covariates on doctor conditional on hospital = 1. I've also manually calculated the conditional effects as the sum of the direct and indirect effect according to Green (1996) using the normal() and binormal() functions at the mean values.

Now I want to extend the model to a multivariate setting. I can recover the model parameters using something like mvprobit or cmp or using John Mullahy's bvpmp() mata function:

mvprobit (doctor = age female educ) (hospital =  age female educ) (public = age female educ)
***
cmp (doctor = age female educ) (hospital =  age female educ) (public = age female educ), ind($cmp_probit $cmp_probit $cmp_probit)
***
mata
yn="doctor hospital public"
xn="age female educ"
bv1=bvpmvp(yn,xn,"",1,.001,1)
end

However, as it well known, neither command supports the recovery of the marginal or conditional effects in any case. What I'd like to do is a) recover the marginal effect of the covariates on the probability of any given outcome doctor=1 or hospital=1 or public=1, and b) recover the marginal effects of x conditional on a vector of other y outcomes, similar to pcond1 but allowing for say, hospital=1 and public = 0.

John Mullahy (2017) "Marginal effects in multivariate probit models" Empirical Economics, extends Green's results to the multivariate case, but I'm having extreme difficulty figuring out how to implement it. Has anyone successfully implemented the APEs or the conditional marginal effects from this paper? I'd be happy just getting manually calculating point estimates given the model parameters without trying to recover the standard errors.

[John Mullahy]

Malcolm: I would recommend using the mvnormalderiv(.) functions in Mata. The analytic code in my EE paper is admittedly awkward to implement, and still needs to call a multivariate normal joint distribution in any event. I've recently been using mvnormalderiv(.)with success, and it sure cleans up the programming algebra. See

help mf_mvnormal

mvnormalderiv(.) will return an M-vector of derivatives w.r.t. the points of evaluation. These can be converted into the PEs using the chain rule and then summing. PEs conditional on other y-outcomes follows essentially the same logic.

[Malcolm Wardlaw]

Thanks John. I'm still on Stata14, so I'll need to upgrade.

I don't want to ask too much, but could you maybe give me a little bit of pseudo-code just to put me on the right track. I have a hard time bridging the analytic representations with the code implementation whenever I'm dealing with matrix heavy operations.

[John Mullahy]

Malcolm: Here is some quickly written code (thus perhaps bug-laden) that gives some idea of how I'd approach this

/*   Suppose there are M=4 outcomes and N=1000 obs.           */
/*   Suppose interest is in PEs of Pr(y1=1,y2=1,y3=1,y4=1|x)  */
 
x=J(1000,1,1),uniform(1000,2)    /*  x matrix=[1,x1,x2]       */
b=uniform(3,4)                   /*  ests. of betas           */
r=.75*I(4)+J(4,4,.25)            /*  est. of correl matrix    */
 
/*  initialize dxb and dr     */
dxb=.
dr=.
 
/*   compute derivatives w.r.t. xb vector                */
mvnormalderiv(x*b,vech(r)',dxb,dr)
 
/*   compute N-vector of PEs w.r.t. x1 using chain rule  */
pe1=dxb*b[2,.]'
 
/*   compute N-vector of PEs w.r.t. x2 using chain rule  */
pe2=dxb*b[3,.]'

Once you have the N-vectors of PEs you can get the APEs by simply computing the means of the column vectors.

Note that this is the case of Pr(1,1,1,1). Other probabilities (e.g. Pr(1,0,0,1), etc.) as well as the conditional-on-some-y's probabilities you mention in #1 require a little more programming. The mvnormalderiv() functions that allow lower and upper limits facilitate the computation of the former, whereas for the conditional-on-some-y's probabilities using mvnormal() alongside mvnormalderiv() should work.

Hope this gives some useful ideas.

[Malcolm Wardlaw]

Thanks a ton John. This is super helpful.

[Jason ZY Liu]

@John Mullahy Hi Dr. Mullahy, thank you for pointing out a direction! I am working your code and I guess your method works fine for marginal effect for joint probability. But if I am only interested in the marginal effect of marginal probability in mvprobit (M=3), how can I get that? Thanks!

[John Mullahy]

By "marginal effect of marginal probability" do you mean dPr(y(m)=1|x)/dx? If so that's just a standard univariate probit.

[Jason ZY Liu]

Yes, I mean the dPr(y(m)=1|x)/dx. I thought the correlations of error terms will affect the derivation.

[Jason ZY Liu]

@John Mullahy, I wonder is there any way to calculate the conditional probability after mvprobit? mvprobit only provides joint probability and marginal probability and it seems to be impossible to calculate the conditional probability if M=3? Am I right? Thanks

[John Mullahy]

There are many possible conditional probabilities to consider with M=3. E.g. Pr(y1|y2,y3), Pr(y1,y2|y3), Pr(y1|y3), etc. etc. Which is/are of particular interest?

[Franci Taale]

I am conducting a research which involves estimating the marginal effects/elasticities of independent parameters in a churdle model. I would like to do this for: (1) participation, (2) consumption, (3) conditional level of consumption, and (4) unconditional level of consumption. I thought this task was something common and already solved elsewhere, but I cannot lay my hands on any suitable commands for execution. Can someone come to my aid?

[Jess Rudnick]

I am following up on this thread to ask for further help of implementing marginal effects, as well as any on insight on the possibility of generating predicted probabilities, after estimating a a multivariate probit (mvprobit) model.

Using the user-written mvprobit and mvppred commands (Cappellari and Jenkins 2003), I have successfully estimated an mvprobit model with 8 dependent binary variables, and 17 covariates (10 are binary dummy variables, 7 continuous). Here is an example of the model structure, though some covariates have been left out for simplicity's sake:

mvprobit (Practice1=Crop Water Irrigation FarmSize Educ Income Years Info) ///
(Practice2=Crop Water Irrigation FarmSize Educ Income Years Info) ///
(Practice3=Crop Water Irrigation FarmSize Educ Income Years Info) ///
(Practice4=Crop Water Irrigation FarmSize Educ Income Years Info) ///
(Practice5=Crop Water Irrigation FarmSize Educ Income Years Info) ///
(Practice6=Crop Water Irrigation FarmSize Educ Income Years Info) ///
(Practice7=Crop Water Irrigation FarmSize Educ Income Years Info) ///
(Practice8=Crop Water Irrigation FarmSize Educ Income Years Info) /// , dr(35)

(1) Per the discussion between John Mullahy and Malcolm Wardlaw above, I am interested in calculating the marginal effect of a 1 unit change in a covariate (e.g. FarmSize) on the probability of success of each Y. I understand the necessity of determining the state of the conditional probabilities (i.e. setting constant the outcome of the other dependent variables), and thanks to John Mullahy's paper (2017) "Marginal effects in multivariate probit models" know that there are 2^8 possible combinations of outcome states, thus I would choose a few scenarios which suit my empirical case best. For example, I would like to calculate dPr( y₁=1 |x)/dx, given Pr(y₁=1| y₂=1, y₃=1, y₄=1, y₅=1, y₆=1, y₇=1, y₈=1), or any other set combination of outcomes. But I do not understand how to implement this in practice.

(2) Additionally, I am interested in calculating the change in predicted probability, for a single dependent variable, given two "scenarios" in covariate values and can implement this manually for a single univariate probit. However, I am not sure how the correlated errors between the equations in the multivariate probit influence those calculations. For example, I would like to adapt the following univariate concept to the multivariate space, and calculate for each dependent variable:

 margins, at(x1 = 0 x2=3 x3=5) atmeans

Thank you for the help!!

[John Mullahy]

Jess:

1. You should be able to do this manually by a combination of calls to Mata's mvnormal(.) function using the parameter estimates retrieved from mvprobit and particular values of the covariates to formulate the arguments to mvnormal(.). E.g. one would compute the conditional probability as a ratio of the joint probability to the joint-marginal probability, Pr(y1|y2,y3,x)=Pr(y1,y2,y3|x)/Pr(y2,y3|x), etc. For discrete changes in covariates this would just need to be evaluated at two different levels of the arguments; for derivatives you'll need chain-rule results that depend in part on a call to Mata's mvnormalderiv(.) function. There may be other approaches but this is the one that occurs to me.

2. The margins results for each dependent variable considered separately, i.e. dPr(y1|x)/dx, don't depend on the correlation structure. (Moreover since the mvprobit results use a normalized correlation matrix (diagonal elements=1) for modeling the covariance structure they are using the same normalization as does standard univariate probit.) You can just use standard univariate probit followed by margins.

[Tuan Nguyen Anh]

Dear @John Mullahy, I would like to attempt the code provided by@Malcolm Wardlaw at the beginning of the thread that was used one of your command bvpmvp() simultaneously with mvprobit. However, I could not apply the command into my code because I could not install the package. Can you please guide me through the installation path? I just want to compare it with mvprobit as I found in this material: https://journals.sagepub.com/doi/pdf...867X1601600107
Thank you in advance.

[John Mullahy]

Hello Tuan. There's nothing to install. bvpmvp() is just a Mata program that calls on some other user-written Mata functions. If you don't have access to the code then send me a private Statalist message with your email address and I'll email it to you.