# NOT RUN {
set.seed(123)
## Simulate survival data
## with random intercepts/slopes and a linear effect of time,
## constant association alpha and no effect of the derivative
d <- simJM(nsub = 200, long_setting = "linear",
alpha_setting = "constant",
dalpha_setting = "zero", full = FALSE)
# remove 0s from the longitudinal time variable to allow sparse calculations
d$obstime <- d$obstime - sqrt(.Machine$double.eps)
## Formula of the according joint model
f <- list(
Surv2(survtime, event, obs = y) ~ s(survtime, bs = "ps"),
gamma ~ s(x1, bs = "ps"),
mu ~ obstime + s(id, bs = "re") +
s(id, obstime, bs = "re"),
sigma ~ 1,
alpha ~ 1,
dalpha ~ -1
)
## Joint model estimation
## jm_bamlss() sets the default optimizer and sampler function.
## First, posterior mode estimates are computed using function
## jm.mode(), afterwards the sampler jm.mcmc() is started.
b <- bamlss(f, data = d, family = "jm",
timevar = "obstime", idvar = "id")
## Plot estimated effects.
plot(b)
## Predict event probabilities for two individuals
## at 12 time units after their last longitudinal measurement.
## The event probability is conditional on their survival
## up to their last observed measurement.
nd <- subset(d, id %in% c(192, 127))
p <- predict(b, type = "probabilities", dt = 12, FUN = c95)
print(p)
## Plot of survival probabilities and
## corresponding longitudinal effects
## for individual id.
jm.survplot(b, id = 3)
jm.survplot(b, id = 30)
## Simulate survival data
## with functional random intercepts and a nonlinear effect
## of time, time-varying association alpha and no effect
## of the derivative.
## Note: This specification is the simJM default.
d <- simJM(nsub = 200, full = FALSE)
## Formula of the according joint model
## number of knots for the smooth nonlinear effect of time
klong <- 8
f <- list(
Surv2(survtime, event, obs = y) ~ s(survtime, bs = "ps"),
gamma ~ s(x1, bs = "ps"),
mu ~ ti(id, bs = "re") +
ti(obstime, bs = "ps", k = klong) +
ti(id, obstime, bs = c("re", "ps"),
k = c(nlevels(d$id), klong)) +
s(x2, bs = "ps"),
sigma ~ 1,
alpha ~ s(survtime, bs = "ps"),
dalpha ~ -1
)
## Estimating posterior mode only using jm.mode()
b_mode <- bamlss(f, data = d, family = "jm",
timevar = "obstime", idvar = "id", sampler = FALSE)
## Estimating posterior means using jm.mcmc()
## with starting values generated from posterior mode
b_mean <- bamlss(f, data = d, family = "jm",
timevar = "obstime", idvar = "id", optimizer = FALSE,
start = parameters(b_mode), results = FALSE)
## Plot effects.
plot(b_mean, model = "alpha")
# }
Run the code above in your browser using DataLab