Carry out Bivariate Probit Model with unbalanced dependent variables in R

Question

I am trying to understand to carryout a Bivariate Probit Model to measure differential levels of of individual victimisation. As such, I am trying to model the outcomes of whether an individual is exposed to crime or not and if this has an effect on an individual being a one off or repeat victim. Following on the approach from "Osborn, D.R., Ellingworth, D., Hope, T. and Trickett, A., 1996. Are repeatedly victimized households different?. Journal of Quantitative Criminology, 12, pp.223-245". See page 233. Here is the data

pdu <- structure(list(polyvictimization9 = c(0L, 1L, 1L, 0L, 0L, 0L, 
0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 
0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 
0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 
0L, 1L, 1L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 
0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), 
    yrsaddr_1 = c(1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 
    0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 
    0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 
    1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 
    0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L
    ), yrsaddr_2 = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 
    1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 
    0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 
    1L), cartot_1 = c(1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 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, 1L, 0L, 0L, 0L, 1L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 
    1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 
    0L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    1L), polyvictimization10 = c(NA, 0L, 1L, NA, NA, NA, NA, 
    0L, 0L, NA, NA, NA, NA, NA, 0L, NA, NA, NA, NA, NA, NA, NA, 
    NA, NA, NA, NA, 0L, NA, NA, NA, 0L, NA, NA, NA, NA, 0L, 0L, 
    NA, NA, NA, NA, NA, 0L, 0L, 0L, NA, NA, NA, NA, 0L, NA, 0L, 
    NA, NA, NA, NA, NA, 0L, 0L, 0L, NA, NA, NA, NA, NA, NA, NA, 
    NA, 0L, 1L, NA, 0L, 0L, NA, NA, 1L, 0L, NA, 0L, NA, NA, 0L, 
    NA, 1L, NA, 0L, NA, NA, 0L, NA, 1L, 0L, 1L, NA, NA, NA, NA, 
    NA, NA, NA, NA), cartot_2 = c(0L, 1L, 0L, 1L, 0L, 0L, 1L, 
    0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 
    0L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 
    1L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 
    1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 
    0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 
    1L, 1L, 0L, 0L)), class = "data.frame", row.names = c(NA, 
-101L))

From the data, the 2 binary dependent variables are polyvictimization9 and polyvictimization10. polyvictimization9 (n = 101) is categorised as 0 = Non victim and 1 = victims (single victim and repeat victims). polyvictimization10 (n = 31) is obtain by recoding polyvictimization9 into 0 = single victim 1 = repeat victim and NA represents non victims).

To execute the bivariate probit model, I used the library(GJRM). The first model is

• polyvictimization9 ~ yrsaddr_1 + yrsaddr_2 + cartot_1

The second model is

• polyvictimization10 ~ cartot_1 + cartot_2

The two models have some similar IVs and distinct IVs as well.

Here is the code

    bl_pdu <- gjrm(list(polyvictimization9 ~ yrsaddr_1 + yrsaddr_2 + cartot_1, 
                    polyvictimization10 ~ cartot_1 + cartot_2),
               data = pdu, 
               margins = c("probit", "probit"),
               model = "B")

This code throws up the error below

    Error in if (check.taus && any(abs(tau) > 0.99999)) stop("some tau is too close to -1 or 1") : 
     missing value where TRUE/FALSE needed
    In addition: Warning message:
    In cor(res1, res2, method = "kendall") : the standard deviation is zero

Please, I would appreciate any help.