In egeminiani/GJRM: Generalised Joint Regression Modelling

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

Generate data with one endogenous binary variable and continuous outcome.

library(GJRM)
set.seed(0)
n <- 1000
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 <- -0.25 - 1.25*y1 + f2(x2) + u[,2]
dataSim <- data.frame(y1, y2, x1, x2)

Recursive Model:

rc <- resp.check(y2, margin = "N", print.par = TRUE, loglik = TRUE)
AIC(rc)
BIC(rc)

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

Semiparametric Recursive Model:

eq.mu.1 <- y1 ~ x1 + s(x2)
eq.mu.2 <- y2 ~ y1 + s(x2)
eq.sigma <- ~ 1
eq.theta <- ~ 1
fl <- list(eq.mu.1, eq.mu.2, eq.sigma, eq.theta)

out <- gjrm(fl, data = dataSim,
            margins = c("probit","N"), 
            gamlssfit = TRUE,
            Model = "B")
conv.check(out)
summary(out)
post.check(out)
jc.probs(out, 1, 1.5, intervals = TRUE)[1:4,]
AT(out, nm.end = "y1")
AT(out, nm.end = "y1", type = "univariate")