In egeminiani/GJRM: Generalised Joint Regression Modelling

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

Joint models with three binary margins

Generate data. Correlation between the two equations 0.5 - Sample size 400

library(GJRM)

set.seed(0)
n <- 400
Sigma <- matrix(0.5, 3, 3); diag(Sigma) <- 1
u     <- rMVN(n, rep(0,3), Sigma)
x1    <- round(runif(n)); x2 <- runif(n); x3 <- runif(n)
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*x1 + f2(x2) + u[,2] > 0, 1, 0)
y3    <- ifelse(-0.75 + 0.25*x1 + u[,3] > 0, 1, 0)
dataSim <- data.frame(y1, y2, y3, x1, x2)

f.l <- list(y1 ~ x1 + s(x2),
            y2 ~ x1 + s(x2),
            y3 ~ x1)

out <- gjrm(f.l, data = dataSim, Model = "T",
            margins = c("probit", "probit", "probit"))

out1 <- gjrm(f.l, data = dataSim, Chol = TRUE, Model = "T",
             margins = c("probit", "probit", "probit"))

conv.check(out)
summary(out)
plot(out, eq = 1)
plot(out, eq = 2)
AIC(out)
BIC(out)

out <- gjrm(f.l, data = dataSim, Model = "T",
            margins = c("probit","logit","cloglog"))

out1 <- gjrm(f.l, data = dataSim, Chol = TRUE, Model = "T",
             margins = c("probit","logit","cloglog"))

conv.check(out)
summary(out)

plot(out, eq = 1)
plot(out, eq = 2)

AIC(out)
BIC(out)

f.l <- list(y1 ~ x1 + s(x2),
            y2 ~ x1 + s(x2),
            y3 ~ x1, 
            ~ 1, ~ 1, ~ 1)

out1 <- gjrm(f.l, data = dataSim, Chol = TRUE, Model = "T",
             margins = c("probit", "probit", "probit"))

f.l <- list(y1 ~ x1 + s(x2),
            y2 ~ x1 + s(x2),
            y3 ~ x1,
            ~ 1, ~ s(x2), ~ 1)

out2 <- gjrm(f.l, data = dataSim, Chol = TRUE, Model = "T",
             margins = c("probit", "probit", "probit"))

f.l <- list(y1 ~ x1 + s(x2),
            y2 ~ x1 + s(x2),
            y3 ~ x1,
            ~ x1, ~ s(x2), ~ x1 + s(x2))

out2 <- gjrm(f.l, data = dataSim, Chol = TRUE, Model = "T",
             margins = c("probit", "probit", "probit"))

f.l <- list(y1 ~ x1 + s(x2),
            y2 ~ x1 + s(x2),
            y3 ~ x1,
            ~ x1, ~ x1, ~ s(x2))

out2 <- gjrm(f.l, data = dataSim, Chol = TRUE, Model = "T",
             margins = c("probit", "probit", "probit"))

f.l <- list(y1 ~ x1 + s(x2),
            y2 ~ x1 + s(x2),
            y3 ~ x1,
            ~ x1, ~ x1 + x2, ~ s(x2))

out2 <- gjrm(f.l, data = dataSim, Chol = TRUE, Model = "T",
             margins = c("probit", "probit", "probit"))

f.l <- list(y1 ~ x1 + s(x2),
            y2 ~ x1 + s(x2),
            y3 ~ x1,
            ~ x1 + x2, ~ x1 + x2, ~ x1 + x2)

out2 <- gjrm(f.l, data = dataSim, Chol = TRUE, Model = "T",
             margins = c("probit", "probit", "probit"))

jcres1 <- jc.probs(out2, 1, 1, 1, type = "joint", cond = 0,
                   intervals = TRUE, n.sim = 100)

nw <- data.frame( x1 = 0, x2 = seq(0, 1, length.out = 100) )
jcres2 <- jc.probs(out2, 1, 1, 1, newdata = nw, type = "joint",
                   cond = 0, intervals = TRUE, n.sim = 100)