In egeminiani/GJRM: Generalised Joint Regression Modelling

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

Generate data with partial observability.

library(GJRM)

set.seed(0)
n     <- 10000
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)
y1    <- ifelse(-1.55 + 2*x1 + x2 + u[,1] > 0, 1, 0)
y2    <- ifelse( 0.45 - x3 + u[,2] > 0, 1, 0)
y     <- y1*y2
dataSim <- data.frame(y, x1, x2, x3)

Bivariate probit with Partial Observability.

out <- gjrm(list(y ~ x1 + x2,
                 y ~ x3),
            data = dataSim, 
            Model = "BPO",
            margins = c("probit", "probit"))
conv.check(out)
summary(out)

Let us have a look at the first ten estimated probabilities for the four events from object out.

cbind(out$p11, out$p10, out$p00, out$p01)[1:10,]

The case with smooth function (more computationally intensive)

f1 <- function(x) cos(pi*2*x) + sin(pi*x)
y1 <- ifelse(-1.55 + 2*x1 + f1(x2) + u[,1] > 0, 1, 0)
y2 <- ifelse( 0.45 - x3 + u[,2] > 0, 1, 0)
y  <- y1*y2
dataSim <- data.frame(y, x1, x2, x3)

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

Let us plot the estimated and true functions.

x2 <- sort(x2); f1.x2 <- f1(x2)[order(x2)] - mean(f1(x2))
plot(out, eq = 1, scale = 0); lines(x2, f1.x2, col = "red")

For additional information, see also ?war.

egeminiani/GJRM documentation built on Sept. 1, 2020, 6:41 p.m.