# How to estimate a bivariate probit (biprobit) model in R with a different set of explanatory variables? [closed]

I'm trying to estimate a bivariate probit model (also called biprobit model) in R where the set of explanatory variables is different for both binary outcomes. Thus, my setting is:

\begin{align} Y_1^* &= X_1 + Y_2 + u_1 \\ Y_1\ &= Y_1^* > 0 \\[8pt] Y_2^* &= X_1 + Z_1 + u_2 \\ Y_2\ &= Y_2^* > 0 \end{align}

where the $$u$$s display error terms (which are multivariate normal distributed), the $$Y$$s are the dependent binary variables, $$X_1$$ and $$Z_1$$ are different explanatory variables. (The $$Y^*$$s are the latent continuous variables.)

In Stata this setting is called "Seemingly unrelated bivariate probit regression" and the command would be biprobit (Y1= X1 Y2) (Y2= X1 Z1)

However, I did not find a solution in R by using Google and https://rseek.org/. I thought the zeligverse package would be a solution --> http://docs.zeligproject.org/articles/zeligchoice_bprobit.html However, it does just work if the set of explanatory variables is the same for both equations, as in the following MWE:

library(zeligverse)
data(sanction)
summary(zelig(cbind(import, export) ~ coop + cost + target,
model="bprobit", data=sanction))


From the documentation I would specify two separate equations like the following

fml <- list(mu1=import ~ coop + cost, mu2=export ~ cost + target)
summary(zelig(formula=fml, model="bprobit", data=sanction))


Then the error Error in formula.default(object, env = baseenv()) : invalid formula appears which has also been discussed in https://github.com/IQSS/Zelig/issues/240 and in particular in https://groups.google.com/forum/#!topic/zelig-statistical-software/n5CQnXeQvAM

Thus, I looked for another package and thought the biprobit function of the mets package would be an alternative --> https://www.rdocumentation.org/packages/mets/versions/1.2.5/topics/biprobit However, I do not understand the syntax (and do also not find adeqaute examples) such that I can specify two equations with a different set of explanatory variables.

Any suggestions?

• Questions that are only about software (e.g. error messages, code or packages, etc.) are generally off topic here. If you have a substantive machine learning or statistical question, please edit to clarify. – gung - Reinstate Monica Apr 12 '19 at 14:24

We can use the GJRM package to estimate bivariate probit models in R.

In Stata we would do

. use sanction, clear
(Written by R)

. biprobit(import = coop cost) (export = cost target), nolog


yielding

Seemingly unrelated bivariate probit            Number of obs     =         78
Wald chi2(4)      =      43.13
Log likelihood = -76.136346                     Prob > chi2       =     0.0000

------------------------------------------------------------------------------
|      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
import       |
coop |   .2549989   .1849965     1.38   0.168    -.1075876    .6175853
cost |   .9459759   .2687549     3.52   0.000     .4192259    1.472726
_cons |  -2.564098   .5364901    -4.78   0.000    -3.615599   -1.512597
-------------+----------------------------------------------------------------
export       |
cost |   1.531134    .344884     4.44   0.000     .8551743    2.207095
target |  -.3305822   .2657153    -1.24   0.213    -.8513747    .1902103
_cons |  -1.816743   .5750259    -3.16   0.002    -2.943773   -.6897134
-------------+----------------------------------------------------------------
/athrho |  -.1438342   .2451419    -0.59   0.557    -.6243035    .3366351
-------------+----------------------------------------------------------------
rho |  -.1428504   .2401395                     -.5541173      .32447
------------------------------------------------------------------------------
LR test of rho=0: chi2(1) = .349229                       Prob > chi2 = 0.5545

.


and in R we would do accordingly

treat.eq <- import ~ coop + cost
out.eq <- export ~ cost + target
f.list <- list(treat.eq, out.eq)
mr <- c("probit", "probit")

library(GJRM)
bvp <- gjrm(f.list, data=sanction, Model="B", margins=mr)
summary(bvp)


yielding

COPULA:   Gaussian
MARGIN 1: Bernoulli
MARGIN 2: Bernoulli

EQUATION 1
Formula: import ~ coop + cost

Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  -2.5641     0.5365  -4.779 1.76e-06 ***
coop          0.2550     0.1850   1.378 0.168080
cost          0.9460     0.2688   3.520 0.000432 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

EQUATION 2
Formula: export ~ cost + target

Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  -1.8167     0.5750  -3.159  0.00158 **
cost          1.5311     0.3449   4.440 9.01e-06 ***
target       -0.3306     0.2657  -1.244  0.21345
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

n = 78  theta = -0.143(-0.543,0.29)  tau = -0.0913(-0.366,0.187)
total edf = 7


The values are pretty much identical.

Data

sanction <- structure(list(mil = c(1L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 1L,
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L,
0L, 0L, 0L, 1L), coop = c(4L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 2L,
2L, 1L, 1L, 3L, 3L, 3L, 1L, 4L, 3L, 1L, 3L, 4L, 1L, 1L, 1L, 1L,
3L, 1L, 1L, 1L, 4L, 1L, 1L, 1L, 1L, 2L, 3L, 1L, 1L, 1L, 2L, 2L,
2L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,
2L, 2L, 3L, 3L, 1L, 3L, 2L, 3L, 2L, 3L, 1L, 1L, 3L, 2L, 2L, 3L,
2L, 1L, 4L, 1L, 3L), target = c(3L, 3L, 3L, 3L, 3L, 3L, 2L, 3L,
1L, 3L, 2L, 2L, 1L, 3L, 2L, 2L, 2L, 3L, 1L, 3L, 3L, 2L, 2L, 2L,
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 3L, 2L, 1L, 2L, 3L,
2L, 3L, 3L, 3L, 1L, 2L, 3L, 1L, 3L, 2L, 2L, 3L, 1L, 1L, 2L, 2L,
2L, 3L, 3L, 2L, 2L, 3L, 1L, 2L, 3L, 1L, 3L, 2L, 2L, 1L, 3L, 2L,
1L, 1L, 3L, 3L, 2L, 1L), import = c(1L, 0L, 1L, 1L, 1L, 0L, 0L,
0L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 1L,
0L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 0L,
0L, 1L, 0L, 0L, 1L, 0L, 1L), export = c(1L, 1L, 0L, 1L, 1L, 1L,
1L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L,
0L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 0L,
0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 1L,
1L, 0L, 1L, 0L, 1L, 0L, 0L, 1L), cost = c(4L, 3L, 2L, 2L, 2L,
2L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 4L, 3L, 1L, 2L, 2L, 1L, 3L, 2L,
1L, 4L, 1L, 1L, 2L, 2L, 2L, 1L, 3L, 1L, 1L, 1L, 1L, 2L, 1L, 1L,
1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 1L,
1L, 1L, 2L, 2L, 2L, 2L, 3L, 1L, 3L, 3L, 1L, 2L, 1L, 3L, 2L, 2L,
2L, 3L, 1L, 2L, 1L, 2L, 2L, 1L, 2L), num = c(15L, 4L, 1L, 1L,
1L, 1L, 3L, 3L, 2L, 1L, 1L, 1L, 8L, 7L, 21L, 1L, 7L, 4L, 1L,
120L, 7L, 1L, 1L, 1L, 1L, 32L, 1L, 1L, 1L, 150L, 1L, 1L, 1L,
5L, 2L, 10L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 3L, 2L, 2L, 9L, 7L, 1L, 10L,
2L, 8L, 2L, 13L, 1L, 1L, 4L, 1L, 8L, 14L, 2L, 1L, 13L, 1L, 10L
), ncost = structure(c(2L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 4L, 1L,
4L, 1L, 1L, 2L, 3L, 4L, 1L, 1L, 4L, 3L, 1L, 4L, 2L, 4L, 4L, 1L,
1L, 1L, 4L, 3L, 4L, 4L, 4L, 4L, 1L, 4L, 4L, 4L, 4L, 1L, 1L, 1L,
1L, 1L, 4L, 4L, 1L, 4L, 4L, 4L, 1L, 1L, 4L, 4L, 4L, 1L, 1L, 1L,
1L, 3L, 4L, 3L, 3L, 4L, 1L, 4L, 3L, 1L, 1L, 1L, 3L, 4L, 1L, 4L,
1L, 1L, 4L, 1L), .Label = c("little effect", "major loss", "modest loss",
"net gain"), class = "factor")), class = "data.frame", row.names = c("1",
"2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13",
"14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24",
"25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35",
"36", "37", "38", "39", "40", "41", "42", "43", "44", "45", "46",
"47", "48", "49", "50", "51", "52", "53", "54", "55", "56", "57",
"58", "59", "60", "61", "62", "63", "64", "65", "66", "67", "68",
"69", "70", "71", "72", "73", "74", "75", "76", "77", "78"))

readstata13::save.dta13(sanction, "sanction.dta")  # save as Stata *.dta

• Wonderful ! That looks really great. Thank you a lot. May I ask if you can also tell me how to get the marginal effects from the bivariate probit regression using this GJRM package. In Stata it would be just margins, dydx(*) pred(pmarg1) for the Average Partial Effects (APE) of the outcome equation. Is there any build-in version for the GJRM package or do I have to calculate them manually? The PE function of the GJRM package seems not appropriate. – Jens Friedrich Apr 15 '19 at 9:45