In egeminiani/GJRM: Generalised Joint Regression Modelling

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  message = FALSE
)
#suppressPackageStartupMessages(library(ggplot2))
#theme_set(theme_light())

GJRM: Generalised Joint Regression Modelling

Authors: Giampiero Marra, Rosalba Radice
License: GPL (>= 2)

The GJRM package offers numerous routines for fitting various joint (and univariate) regression models, with several types of covariate effects, in the presence of equations' errors association, endogeneity, non-random sample selection or partial observability.

Installation

You can install the GJRM package from CRAN:

install.packages("GJRM")

Joint models with binary margins: the `gjrm` function

Let us consider a simple example that shows how to use the gjrm function. First, we 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)

Classic bivariate probit model

We can use the gjrm function to estimate a classic bivariate probit model 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 calculated using the AIC and BIC functions.

AIC(out)
BIC(out)

Bivariate probit model with splines

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")

Bivariate probit model with splines and varying dependence parameter

A bivariate probit model 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)