In egeminiani/GJRM: Generalised Joint Regression Modelling

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

Example using Gumbel copula and normal-gamma margins

library(GJRM)

set.seed(1)
n <- 2000

rh <- 0.5      

sigmau <- matrix(c(1, rh, rh, 1), 2, 2)
u      <- rMVN(n, rep(0,2), sigmau)

sigmac <- matrix(rh, 3, 3); diag(sigmac) <- 1
cov    <- rMVN(n, rep(0,3), sigmac)
cov    <- pnorm(cov)

bi <- round(cov[,1]); x1 <- cov[,2]; x2 <- cov[,3]

f11 <- function(x) -0.7*(4*x + 2.5*x^2 + 0.7*sin(5*x) + cos(7.5*x))
f12 <- function(x) -0.4*( -0.3 - 1.6*x + sin(5*x))  
f21 <- function(x) 0.6*(exp(x) + sin(2.9*x)) 

ys <-  0.58 + 2.5*bi + f11(x1) + f12(x2) + u[, 1] > 0
y  <- exp(-0.68 - 1.5*bi + f21(x1) + u[, 2])
yo <- y*(ys > 0)

dataSim <- data.frame(ys, yo, bi, x1, x2)


out <- gjrm(list(ys ~ bi + s(x1) + s(x2), 
                 yo ~ bi + s(x1)), 
                 data = dataSim, BivD = "G0", 
                 margins = c("probit", "GA"),
                 Model = "BSS")
conv.check(out)
post.check(out)
summary(out)

ATE <- NA
n.m <- 10
res <- imputeSS(out, n.m)

for(i in 1:n.m){
  dataSim$yo[dataSim$ys == 0] <- res[[i]]
  outg <- gamlss(list(yo ~ bi + s(x1)), margin = "GA", data = dataSim)
  out$gamlss <- outg
  ATE[i] <- AT(out, nm.end = "bi", type = "univariate")$res[2]
  print(i)
}

AT(out, nm.end = "bi")
mean(ATE)