# NOT RUN {
# }
# NOT RUN {
pbc2.id$Time <- pbc2.id$years
pbc2.id$event <- as.numeric(pbc2.id$status != "alive")
##########################################################################################
##############################################
# Univariate joint model for serum bilirubin #
##############################################
# [1] Fit the mixed model using mvglmer(). The main arguments of this function are:
# 'formulas' a list of lme4-like formulas (a formular per outcome),
# 'data' a data.frame that contains all the variables specified in 'formulas' (NAs allowed),
# 'families' a list of family objects specifying the type of each outcome (currently only
# gaussian, binomial and poisson are allowed).
MixedModelFit1 <- mvglmer(list(log(serBilir) ~ year + (year | id)), data = pbc2,
families = list(gaussian))
# [2] Fit a Cox model, specifying the baseline covariates to be included in the joint
# model; you need to set argument 'model' to TRUE.
CoxFit1 <- coxph(Surv(Time, event) ~ drug + age, data = pbc2.id, model = TRUE)
# [3] The basic joint model is fitted using a call to mvJointModelBayes(), which is very
# similar to JointModelBayes(), i.e.,
JMFit1 <- mvJointModelBayes(MixedModelFit1, CoxFit1, timeVar = "year")
summary(JMFit1)
plot(JMFit1)
##########################################################################################
#########################################################
# Bivariate joint model for serum bilirubin and spiders #
#########################################################
MixedModelFit2 <- mvglmer(list(log(serBilir) ~ year + (year | id),
spiders ~ year + (1 | id)), data = pbc2,
families = list(gaussian, binomial))
CoxFit2 <- coxph(Surv(Time, event) ~ drug + age, data = pbc2.id, model = TRUE)
JMFit2 <- mvJointModelBayes(MixedModelFit2, CoxFit2, timeVar = "year")
summary(JMFit2)
plot(JMFit2)
##########################################################################################
#######################
# slopes & area terms #
#######################
# We extend model 'JMFit2' by including the value and slope term for bilirubin, and
# the area term for spiders (in the log-odds scale). To include these terms into the model
# we specify the 'Formulas' argument. This is specified in a similar manner as the
# 'derivForms' argument of jointModelBayes(). In particular, it should be a list of lists.
# Each component of the outer list should have as name the name of the corresponding
# outcome variable. Then in the inner list we need to specify four components, namely,
# 'fixed' & 'random' R formulas specifying the fixed and random effects part of the term
# to be included, and 'indFixed' & 'indRandom' integer indicices specifying which of the
# original fixed and random effects are involved in the claculation of the new term. In
# the inner list you can also optionally specify a name for the term you want to include.
# Notes: (1) For terms not specified in the 'Formulas' list, the default value functional
# form is used. (2) If for a particular outcome you want to include both the value
# functional form and an extra term, then you need to specify that in the 'Formulas'
# using two entries. To include the value functional form you only need to set the
# corresponding to 'value', and for the second term to specify the inner list. See
# example below on how to include the value and slope for serum bilirubin (for example,
# if the list below the entry '"log(serBilir)" = "value"' was not give, then only the
# slope term would have been included in the survival submodel).
Forms <- list("log(serBilir)" = "value",
"log(serBilir)" = list(fixed = ~ 1, random = ~ 1,
indFixed = 2, indRandom = 2, name = "slope"),
"spiders" = list(fixed = ~ 0 + year + I(year^2/2), random = ~ 0 + year,
indFixed = 1:2, indRandom = 1, name = "area"))
JMFit3 <- update(JMFit2, Formulas = Forms)
summary(JMFit3)
plot(JMFit3)
##########################################################################################
#####################
# Interaction terms #
#####################
# We further extend the previous model by including the interactions terms between the
# terms specified in 'Formulas' and 'drug'. The names specified in the list that defined
# the interaction factors should match with the names of the output from 'JMFit3'; the
# only exception is with regard to the 'value' functional form. See specification below
# (to include the interaction of the value term of 'log(serBilir)' with 'drug', in the
# list we can either specify as name of the corresponding formula 'log(serBilir)_value'
# or just 'log(serBilir)'):
Ints <- list("log(serBilir)" = ~ drug, "log(serBilir)_slope" = ~ drug,
"spiders_area" = ~ drug)
# because of the many association parameters we have, we place a shrinkage prior on the
# alpha coefficients. In particular, if we have K association parameters, we assume that
# alpha_k ~ N(0, tau * phi_k), k = 1, ..., K. The precision parameters tau and phi_k are
# given Gamma priors. Precision tau is global shrinkage parameter, and phi_k a specific
# per alpha coefficient.
JMFit4 <- update(JMFit3, Interactions = Ints, priors = list(shrink_alphas = TRUE))
summary(JMFit4)
plot(JMFit4)
##########################################################################################
########################
# Time-varying effects #
########################
# We allow the association parameter to vary with time; the time-varying coefficients are
# approximated with P-splines
JMFit_tveffect <- mvJointModelBayes(MixedModelFit1, CoxFit1, timeVar = "year",
Interactions = list("log(serBilir)_value" = ~ 0 + tve(Time, df = 8)))
plot(JMFit_tveffect, "tv_effect")
##########################################################################################
############################
# Interval censoring terms #
############################
# create artificial interval censoring in the PBC data set
pbc2$status2 <- as.numeric(pbc2$status != "alive")
pbc2.id$status2 <- as.numeric(pbc2.id$status != "alive")
pbc2$status3 <- as.character(pbc2$status)
ff <- function (x) {
out <- if (x[1L] %in% c('dead', 'transplanted')) 'dead' else
switch(sample(1:3, 1), '1' = "alive", '2' = "left", '3' = "interval")
rep(out, length(x))
}
pbc2$status3 <- unlist(with(pbc2, lapply(split(status3, id), ff)), use.names = FALSE)
pbc2$status3 <- unname(with(pbc2, sapply(status3, function (x)
switch(x, 'dead' = 1, 'alive' = 0, 'left' = 2, 'interval' = 3))))
pbc2$yearsU <- as.numeric(NA)
pbc2$yearsU[pbc2$status3 == 3] <- pbc2$years[pbc2$status3 == 3] +
runif(sum(pbc2$status3 == 3), 0, 4)
pbc2.id <- pbc2[!duplicated(pbc2$id), ]
# next we fit a weibull model for interval censored data
survFit <- survreg(Surv(years, yearsU, status3, type = "interval") ~ drug + age,
data = pbc2.id, model = TRUE)
# next we fit the joint model (we use 'MixedModelFit1' from above)
JMFit_intcens <- mvJointModelBayes(MixedModelFit1, survFit, timeVar = "year")
summary(JMFit_intcens)
# }
Run the code above in your browser using DataLab