In egeminiani/GJRM: Generalised Joint Regression Modelling

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

Survival Models

library(GJRM)

set.seed(0)
n  <- 2000
c  <- runif(n, 3, 8)
u  <- runif(n, 0, 1)
z1 <- rbinom(n, 1, 0.5)
z2 <- runif(n, 0, 1)
t  <- rep(NA, n)
beta_0 <- -0.2357
beta_1 <- 1

f <- function(t, beta_0, beta_1, u, z1, z2){
  S_0 <- 0.7 * exp(-0.03*t^1.9) + 0.3*exp(-0.3*t^2.5)
  exp(-exp(log(-log(S_0))+beta_0*z1 + beta_1*z2))-u
}

for (i in 1:n){
  t[i] <- uniroot(f, c(0, 8), tol = .Machine$double.eps^0.5,
                  beta_0 = beta_0, beta_1 = beta_1, u = u[i],
                  z1 = z1[i], z2 = z2[i], extendInt = "yes" )$root
}

delta1 <- ifelse(t < c, 1, 0)
u1     <- apply(cbind(t, c), 1, min)
dataSim <- data.frame(u1, delta1, z1, z2)

c <- runif(n, 4, 8)
u <- runif(n, 0, 1)
z <- rbinom(n, 1, 0.5)
beta_0 <- -1.05
t <- rep(NA, n)
f <- function(t, beta_0, u, z){
  S_0 <- 0.7 * exp(-0.03*t^1.9) + 0.3*exp(-0.3*t^2.5)
  1/(1 + exp(log((1-S_0)/S_0)+beta_0*z))-u
}

for (i in 1:n){
  t[i] <- uniroot(f, c(0, 8), tol = .Machine$double.eps^0.5,
                  beta_0 = beta_0, u = u[i], z = z[i],
                  extendInt="yes" )$root
}

delta2 <- ifelse(t < c,1, 0)
u2     <- apply(cbind(t, c), 1, min)
dataSim$delta2 <- delta2
dataSim$u2 <- u2
dataSim$z <- z
eq1 <- u1 ~ s(u1, bs = "mpi") + z1 + s(z2)
eq2 <- u2 ~ s(u2, bs = "mpi") + z
eq3 <- ~ s(z2)

out <- gjrm(list(eq1, eq2), data = dataSim, surv = TRUE,
            margins = c("PH", "PO"),
            cens1 = delta1, cens2 = delta2, Model = "B")

PH margin fit can also be compared with cox.ph from mgcv.

conv.check(out)
res <- post.check(out)

Martingale residuals:

mr1 <- out$cens1 - res$qr1
mr2 <- out$cens2 - res$qr2

can be plotted against covariates obs index, survival time, rank order of surv times.

# To determine func form, one may use # res from null model against covariate.
# To test for PH, use:
# library(survival)
# fit <- coxph(Surv(u1, delta1) ~ z1 + z2, data = dataSim)
# temp <- cox.zph(fit)
# print(temp)
# plot(temp, resid = FALSE)

summary(out)
AIC(out)
BIC(out)

plot(out, eq = 1, scale = 0, pages = 1)
plot(out, eq = 2, scale = 0, pages = 1)

hazsurv.plot(out, eq = 1, newdata = data.frame(z1 = 0, z2 = 0),
             shade = TRUE, n.sim = 1000)

hazsurv.plot(out, eq = 1, newdata = data.frame(z1 = 0, z2 = 0),
             shade = TRUE, n.sim = 1000, type = "hazard")
hazsurv.plot(out, eq = 2, newdata = data.frame(z = 0),
             shade = TRUE, n.sim = 1000)
hazsurv.plot(out, eq = 2, newdata = data.frame(z = 0),
             shade = TRUE, n.sim = 1000, type = "hazard")

jc.probs(out, type = "joint", intervals = TRUE)[1:5,]
newd0 <- newd1 <- data.frame(z = 0, z1 = mean(dataSim$z1),
                             z2 = mean(dataSim$z2),
                             u1 = mean(dataSim$u1) + 1,
                             u2 = mean(dataSim$u2) + 1)
newd1$z <- 1

jc.probs(out, type = "joint", newdata = newd0, intervals = TRUE)
jc.probs(out, type = "joint", newdata = newd1, intervals = TRUE)

out1 <- gjrm(list(eq1, eq2, eq3), data = dataSim, surv = TRUE,
             margins = c("PH", "PO"),
             cens1 = delta1, cens2 = delta2, gamlssfit = TRUE,
             Model = "B")


# eq1 <- u1 ~ z1 + s(z2)
# eq2 <- u2 ~ z
# eq3 <- ~ s(z2)
# note that Weibull is implemented as AFT model (test case)
# out2 <- gjrm(list(eq1, eq2, ~ 1, ~ 1, eq3), data = dataSim, surv = TRUE,
# margins = c("WEI", "WEI"),
# cens1 = delta1, cens2 = delta2,
# Model = "B")

### Joint continuous and survival outcomes
# work in progress
#
# eq1 <- z2 ~ z1
# eq2 <- u2 ~ s(u2, bs = "mpi") + z
# eq3 <- ~ s(z2)
# eq4 <- ~ s(z2)
#
# f.l <- list(eq1, eq2, eq3, eq4)
#
# out3 <- gjrm(f.l, data = dataSim, surv = TRUE,
# margins = c("N", "PO"),
# cens1 = NULL, cens2 = delta2,
# gamlssfit = TRUE, Model = "B")
#
# conv.check(out3)
# post.check(out3)
# summary(out3)
# AIC(out3); BIC(out3)
# plot(out3, eq = 2, scale = 0, pages = 1)
# plot(out3, eq = 3, scale = 0, pages = 1)
# plot(out3, eq = 4, scale = 0, pages = 1)
#
# newd <- newd1 <- data.frame(z = 0, z1 = mean(dataSim$z1),
# z2 = mean(dataSim$z2),
# u2 = mean(dataSim$u2) + 1)
#
# jc.probs(out3, y1 = 0.6, type = "joint", newdata = newd, intervals = TRUE)