In egeminiani/GJRM: Generalised Joint Regression Modelling

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

Let us generate the data from a joint model with binary margins:

library(GJRM)

set.seed(0)
n     <- 400
Sigma <- matrix(0.5, 2, 2); diag(Sigma) <- 1
u     <- rMVN(n, rep(0,2), 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)
dataSim <- data.frame(y1, y2, x1, x2, x3)

A classic bivariate probit model can be easily estimated through the gjrm function as follows.

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

The convergence status can be examined through the conv.check function.

conv.check(out)

An overview of the parameter estimates is available through the summary method.

summary(out)

Information criteria are readily found using the functions AIC and BIC.

AIC(out)
BIC(out)

A bivariate probit model with Splines can be estimated in the following way.

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

We can have a look at the estimated smooth function plots. The red lines indicate the true curves.

x2    <- sort(x2)
f1.x2 <- f1(x2)[order(x2)] - mean(f1(x2))
f2.x2 <- f2(x2)[order(x2)] - mean(f2(x2))
f3.x3 <- rep(0, length(x3))

par(mfrow=c(2,2),mar=c(4.5,4.5,2,2))

plot(out, eq = 1, select = 1, seWithMean = TRUE, scale = 0)
lines(x2, f1.x2, col = "red")

plot(out, eq = 1, select = 2, seWithMean = TRUE, scale = 0)
lines(x3, f3.x3, col = "red")

plot(out, eq = 2, select = 1, seWithMean = TRUE, scale = 0)
lines(x2, f2.x2, col = "red")

plot(out, eq = 2, select = 2, seWithMean = TRUE, scale = 0)
lines(x3, f3.x3, col = "red")

A bivariate probit with Splines and a varying dependence parameter can be fitted as follows.

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

outD     <- gjrm(fl, data = dataSim,
                 margins = c("probit", "probit"),
                 Model = "B")
conv.check(outD)
summary(outD)
# outD$theta

And similarly, the estimated smooth functions are plotted as follows.

par(mfrow=c(2,2),mar=c(4.5,4.5,2,2))

plot(outD, eq = 1, seWithMean = TRUE)
plot(outD, eq = 2, seWithMean = TRUE)
plot(outD, eq = 3, seWithMean = TRUE)