3
$\begingroup$

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?

$\endgroup$

closed as off-topic by gung Apr 12 at 14:23

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question appears to be off-topic because EITHER it is not about statistics, machine learning, data analysis, data mining, or data visualization, OR it focuses on programming, debugging, or performing routine operations within a statistical computing platform. If the latter, you could try the support links we maintain." – gung
If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\begingroup$ 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. $\endgroup$ – gung Apr 12 at 14:24
2
$\begingroup$

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
Link function for mu.1: probit 
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
Link function for mu.2: probit 
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
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Jens Friedrich Apr 15 at 9:45

Not the answer you're looking for? Browse other questions tagged or ask your own question.