In egeminiani/GJRM: Generalised Joint Regression Modelling

knitr::opts_chunk$set(
  warning = FALSE,
  message = FALSE,
  collapse = TRUE,
  comment = "#>"
)

Let us generate data with one endogenous variable and exclusion restriction.

library(GJRM)
set.seed(0)
n     <- 400
Sigma <- matrix(0.5, 2, 2); diag(Sigma) <- 1
u     <- rMVN(n, rep(0,2), Sigma)
cov   <- rMVN(n, rep(0,2), Sigma)
cov   <- pnorm(cov)
x1    <- round(cov[,1]); x2 <- cov[,2]
f1    <- function(x) cos(pi*2*x) + sin(pi*x)
f2    <- function(x) x+exp(-30*(x-0.5)^2)
y1    <- ifelse(-1.55 + 2*x1 + f1(x2) + u[,1] > 0, 1, 0)
y2    <- ifelse(-0.25 - 1.25*y1 + f2(x2) + u[,2] > 0, 1, 0)
dataSim <- data.frame(y1, y2, x1, x2)

The hypothesis of absence of endogeneity can be tested through the LM.bpm function.

LM.bpm(list(y1 ~ x1 + s(x2), y2 ~ y1 + s(x2)), dataSim, Model = "B")

A classic recursive bivariate probit is fitted as follows.

out <- gjrm(list(y1 ~ x1 + x2,
                 y2 ~ y1 + x2),
            data = dataSim,
            margins = c("probit", "probit"),
            Model = "B")
conv.check(out)
summary(out)
AIC(out)
BIC(out)

A flexible recursive bivariate probit model can be fitted as follows.

out <- gjrm(list(y1 ~ x1 + s(x2),
                 y2 ~ y1 + s(x2)),
            data = dataSim,
            margins = c("probit", "probit"),
            Model = "B")

conv.check(out)
summary(out)
AIC(out)
BIC(out)

Testing the hypothesis of absence of endogeneity post estimation:

gt.bpm(out)

Treatment effect, risk ratio and odds ratio with CIs are obtained as:

mb(y1, y2, Model = "B")
AT(out, nm.end = "y1", hd.plot = TRUE)
RR(out, nm.end = "y1")
OR(out, nm.end = "y1")
AT(out, nm.end = "y1", type = "univariate")
re.imp <- imputeCounter(out, m = 10, "y1")
re.imp$AT

Try a Clayton copula model...

outC <- gjrm(list(y1 ~ x1 + s(x2),
                  y2 ~ y1 + s(x2)),
             data = dataSim, BivD = "C0",
             margins = c("probit", "probit"),
             Model = "B")
conv.check(outC)
summary(outC)
AT(outC, nm.end = "y1")
re.imp <- imputeCounter(outC, m = 10, "y1")
re.imp$AT

Try a Joe copula model...

outJ <- gjrm(list(y1 ~ x1 + s(x2),
                  y2 ~ y1 + s(x2)),
             data = dataSim, BivD = "J0",
             margins = c("probit", "probit"),
             Model = "B")
conv.check(outJ)
summary(outJ)
AT(outJ, "y1")
re.imp <- imputeCounter(outJ, m = 10, "y1")
re.imp$AT
VuongClarke(out, outJ)

Recursive bivariate probit modelling with unpenalized splines can be achieved as follows.

outFP <- gjrm(list(y1 ~ x1 + s(x2, bs = "cr", k = 5),
                   y2 ~ y1 + s(x2, bs = "cr", k = 6)),
              fp = TRUE, data = dataSim,
              margins = c("probit", "probit"),
              Model = "B")
conv.check(outFP)
summary(outFP)