# NOT RUN {
# }
# NOT RUN {
# A joint model for the AIDS dataset:
# First we fit the linear mixed model for the longitudinal measurements of
# sqrt CD4 cell counts
lmeFit.aids <- lme(CD4 ~ obstime * drug, random = ~ obstime | patient, data = aids)
# next we fit the Cox model for the time to death (including the 'x = TRUE' argument)
survFit.aids <- coxph(Surv(Time, death) ~ drug, data = aids.id, x = TRUE)
# the corresponding joint model is fitted by (the default is to assume
# the current value parameterization)
jointFit.aids <- jointModelBayes(lmeFit.aids, survFit.aids, timeVar = "obstime")
summary(jointFit.aids)
# A joint model for the PBC dataset:
# We first fit the linear mixed and Cox models. In the first we include
# splines to model flexibly the subject-specific longitudinal trajectories
lmeFit.pbc <- lme(log(serBilir) ~ ns(year, 2),
random = list(id = pdDiag(form = ~ ns(year, 2))), data = pbc2)
survFit.pbc <- coxph(Surv(years, status2) ~ 1, data = pbc2.id, x = TRUE)
# the corresponding joint model is fitted by:
jointFit.pbc <- jointModelBayes(lmeFit.pbc, survFit.pbc, timeVar = "year",
baseHaz = "regression-splines")
summary(jointFit.pbc)
# we update the joint model fitted for the PBC dataset by including
# the time-dependent slopes term. To achieve this we need to define
# the 'extraForm' argument, in which we use function dns() to numerically
# compute the derivative of the natural cubic spline. In addition, we increase
# the number of MCMC iterations to 35000
dform = list(fixed = ~ 0 + dns(year, 2), random = ~ 0 + dns(year, 2),
indFixed = 2:3, indRandom = 2:3)
jointFit.pbc2 <- update(jointFit.pbc, param = "td-both", extraForm = dform,
n.iter = 35000)
summary(jointFit.pbc2)
# we fit the same model with the shared random effects formulation
jointFit.pbc3 <- update(jointFit.pbc, param = "shared-betasRE")
summary(jointFit.pbc3)
# a joint model for left censored longitudinal data
# we create artificial left censoring in serum bilirubin
pbc2$CensInd <- as.numeric(pbc2$serBilir <= 0.8)
pbc2$serBilir2 <- pbc2$serBilir
pbc2$serBilir2[pbc2$CensInd == 1] <- 0.8
censdLong <- function (y, eta.y, scale, log = FALSE, data) {
log.f <- dnorm(x = y, mean = eta.y, sd = scale, log = TRUE)
log.F <- pnorm(q = y, mean = eta.y, sd = scale, log.p = TRUE)
ind <- data$CensInd
if (log) {
(1 - ind) * log.f + ind * log.F
} else {
exp((1 - ind) * log.f + ind * log.F)
}
}
lmeFit.pbc2 <- lme(log(serBilir2) ~ ns(year, 2), data = pbc2,
random = ~ ns(year, 2) | id, method = "ML")
jointFit.pbc4 <- jointModelBayes(lmeFit.pbc2, survFit.pbc, timeVar = "year",
densLong = censdLong)
summary(jointFit.pbc4)
# a joint model for a binary outcome
pbc2$serBilirD <- 1 * (pbc2$serBilir > 1.8)
lmeFit.pbc3 <- glmmPQL(serBilirD ~ year, random = ~ year | id,
family = binomial, data = pbc2)
dLongBin <- function (y, eta.y, scale, log = FALSE, data) {
dbinom(x = y, size = 1, prob = plogis(eta.y), log = log)
}
jointFit.pbc5 <- jointModelBayes(lmeFit.pbc3, survFit.pbc, timeVar = "year",
densLong = dLongBin)
summary(jointFit.pbc5)
# create start-stop counting process variables
pbc <- pbc2[c("id", "serBilir", "drug", "year", "years",
"status2", "spiders")]
pbc$start <- pbc$year
splitID <- split(pbc[c("start", "years")], pbc$id)
pbc$stop <- unlist(lapply(splitID,
function (d) c(d$start[-1], d$years[1]) ))
pbc$event <- with(pbc, ave(status2, id,
FUN = function (x) c(rep(0, length(x)-1), x[1])))
pbc <- pbc[!is.na(pbc$spiders), ]
# left-truncation
pbc <- pbc[pbc$start != 0, ]
lmeFit.pbc <- lme(log(serBilir) ~ drug * ns(year, 2),
random = ~ ns(year, 2) | id, data = pbc)
tdCox.pbc <- coxph(Surv(start, stop, event) ~ drug * spiders + cluster(id),
data = pbc, x = TRUE, model = TRUE)
jointFit.pbc6 <- jointModelBayes(lmeFit.pbc, tdCox.pbc, timeVar = "year")
summary(jointFit.pbc6)
# }
Run the code above in your browser using DataLab