In egeminiani/GJRM: Generalised Joint Regression Modelling

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

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

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     <- -1.55 + 2*x1 + f1(x2) + u[,1]
y2     <- -0.25 - 1.25*x1 + f2(x2) + u[,2]
dataSim <- data.frame(y1, y2, x1, x2, x3)

resp.check(y1, "N")
resp.check(y2, "N")

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

The order above is the one to follow when using more than two equations.

out <- gjrm(fl, data = dataSim, margins = c("N", "N"),
            Model = "B")
conv.check(out)
post.check(out)
summary(out)
AIC(out)
BIC(out)
jc.probs(out, 1.4, 2.3, intervals = TRUE)[1:4,]