In egeminiani/GJRM: Generalised Joint Regression Modelling

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

Generate data with one endogenous continuous exposure and binary 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    <- -0.25 - 2*x1 + f2(x2) + u[,2]
y2    <- ifelse(-0.25 - 0.25*y1 + f1(x2) + u[,1] > 0, 1, 0)
dataSim <- data.frame(y1, y2, x1, x2)
eq.mu.1 <- y2 ~ y1 + s(x2)
eq.mu.2 <- y1 ~ x1 + 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"),
            Model = "B")
conv.check(out)
summary(out)
post.check(out)
AT(out, nm.end = "y1")
AT(out, nm.end = "y1", type = "univariate")
RR(out, nm.end = "y1", rr.plot = TRUE)
RR(out, nm.end = "y1", type = "univariate")
OR(out, nm.end = "y1", or.plot = TRUE)
OR(out, nm.end = "y1", type = "univariate")