In drizopoulos/JMbayes2: Extended Joint Models for Longitudinal and Time-to-Event Data
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
library("JMbayes2")
Non-Gaussian Joint Models with JMbayes2
Taking advantage of the versatility of the GLMMadaptive package, JMbayes2 can fit joint models with several different types of mixed-effects models. The following examples illustrate these capabilities. All examples have the same structure, namely, first, a short motivation for each mixed-model is given, followed by a piece of R code simulating data from a joint model with the respective mixed-effects sub-model, closing by the syntax to fit the joint model. In this last part, the main difference per example is the call to mixed_model().
Beta mixed models
With very few exceptions, continuous outcomes that we wish to analyze have some natural bounds. For example, the levels of blood biomarkers for a set of patients. However, most observations are often located far away from these natural bounds, and an assumption of a normal distribution for the outcome can be safely made. In some settings, though, we can have outcomes for which a substantial percentage of the observations are located near the boundaries leading to skewed or U-shaped distributions. A linear mixed model with normal error terms often does not provide a good fit for such longitudinal outcomes. A natural alternative is to select a distribution that respects the bounded nature of the outcome. The most well-known distribution for such outcomes is the Beta distribution defined in the $(0, 1)$ interval (note: a bounded outcome $Y^$ in the $(a, b)$ interval can be transformed to the $Y = (Y^ - a) / (b - a)$ in the $(0, 1)$ interval).
The following code illustrates how to simulate data from a joint model with a Beta mixed effects model. The default functional form is assumed, i.e., that the linear predictor $\eta(t)$ of the mixed model is associated with the hazard of an event at time $t$. The linear predictor is related to the mean $\mu(t)$ of the Beta distribution under the logit link function, i.e., $\log[\mu(t) / {1 - \mu(t)}] = \eta(t)$.
set.seed(1234)
n <- 200 # number of subjects
K <- 8 # number of measurements per subject
t_max <- 10 # maximum follow-up time
# we construct a data frame with the design:
# everyone has a baseline measurement, and then measurements at random
# follow-up times up to t_max
DF <- data.frame(id = rep(seq_len(n), each = K),
time = c(replicate(n, c(0, sort(runif(K - 1, 0, t_max))))),
sex = rep(gl(2, n/2, labels = c("male", "female")), each = K))
# design matrices for the fixed and random effects
X <- model.matrix(~ sex * time, data = DF)
Z <- model.matrix(~ time, data = DF)
betas <- c(-2.2, -0.25, 0.24, -0.05) # fixed effects coefficients
phi <- 5 # precision parameter of the Beta distribution
D11 <- 1.0 # variance of random intercepts
D22 <- 0.5 # variance of random slopes
# we simulate random effects
b <- cbind(rnorm(n, sd = sqrt(D11)), rnorm(n, sd = sqrt(D22)))
# linear predictor
eta_y <- as.vector(X %*% betas + rowSums(Z * b[DF$id, ]))
# mean of the Beta distribution
mu_y <- plogis(eta_y) # plogis(eta_y) = exp(eta_y) / (1 + exp(eta_y))
# we simulate Beta longitudinal data
DF$y <- rbeta(n * K, shape1 = mu_y * phi, shape2 = phi * (1 - mu_y))
# we transform to (0, 1)
DF$y <- (DF$y * (nrow(DF) - 1) + 0.5) / nrow(DF)
upp_Cens <- 15 # fixed Type I censoring time
shape_wb <- 5 # shape Weibull
alpha <- 0.8 # association coefficients
gammas <- c("(Intercept)" = -9, "sex" = 0.5)
W <- model.matrix(~ sex, data = DF[!duplicated(DF$id), ])
# linear predictor for the survival model
eta_t <- as.vector(W %*% gammas)
# to simulate event times we use inverse transform sampling
# (https://en.wikipedia.org/wiki/Inverse_transform_sampling). Namely, we want
# to find t, such that S(t) = u, where S(.) is the survival function, and u a
# number from the Unif(0, 1) distribution. The function below calculates
# log(u) - log(S(t)), and for a given u, we want to find t for which it equals
# zero. We do that below using the uniroot() function
invS <- function (t, i) {
# i denotes the subject
sex_i <- W[i, 2L]
# h() is the hazard function and we assume a Weibull baseline hazard
h <- function (s) {
X_at_s <- cbind(1, sex_i, s, sex_i * s)
Z_at_s <- cbind(1, s)
# the linear predictor from the mixed model evaluated at time s
f <- as.vector(X_at_s %*% betas +
rowSums(Z_at_s * b[rep(i, nrow(Z_at_s)), ]))
exp(log(shape_wb) + (shape_wb - 1) * log(s) + eta_t[i] + f * alpha)
}
# -log(S(t)) = H(t), where H(t) is the cumulative hazard function
integrate(h, lower = 0, upper = t)$value + log(u[i])
}
# we simulate the event times
u <- runif(n)
trueTimes <- numeric(n)
for (i in seq_len(n)) {
Up <- 100
Root <- try(uniroot(invS, interval = c(1e-05, Up), i = i)$root, TRUE)
trueTimes[i] <- if (!inherits(Root, "try-error")) Root else 150
}
# we use fixed Type I right censoring denoting the end of the trial.
Ctimes <- upp_Cens
Time <- pmin(trueTimes, Ctimes)
event <- as.numeric(trueTimes <= Ctimes) # event indicator
# we keep the longitudinal measurements before the event times
DF$Time <- Time[DF$id]
DF$event <- event[DF$id]
DF <- DF[DF$time <= DF$Time, ]
To fit the corresponding joint model, we fit first a Beta mixed model using the beta.fam() family object into the call of mixed_model():
DF_id <- DF[!duplicated(DF$id), ]
Cox_fit <- coxph(Surv(Time, event) ~ sex, data = DF_id)
Beta_MixMod <- mixed_model(y ~ sex * time, random = ~ time | id, data = DF,
family = beta.fam())
jointFit <- jm(Cox_fit, Beta_MixMod, time_var = "time")
summary(jointFit)
Censored linear mixed models
Some continuous longitudinal outcomes may have a censored nature. A typical example of such outcomes is when we have limit of detection issues. That is, the values of the outcome cannot be detected below a specified threshold having to do with the (laboratory) equipment used to determine the measurements. In these settings, and even if the complete data follows a normal distribution the observed censored data cannot be analyzed with a standard mixed model. The mixed_model() function can accommodate such outcomes using the censored.normal() family object.
The following piece of code simulates data from a joint model with a linear mixed model for the longitudinal outcomes but applies censoring in the realized longitudinal observations.
set.seed(1234)
n <- 200 # number of subjects
K <- 12 # number of measurements per subject
t_max <- 14 # maximum follow-up time
# we construct a data frame with the design:
# everyone has a baseline measurement, and then measurements at random
# follow-up times up to t_max
DF <- data.frame(id = rep(seq_len(n), each = K),
time = c(replicate(n, c(0, sort(runif(K - 1, 0, t_max))))),
sex = rep(gl(2, n/2, labels = c("male", "female")), each = K))
# design matrices for the fixed and random effects
X <- model.matrix(~ sex * time, data = DF)
Z <- model.matrix(~ time, data = DF)
betas <- c(-2.2, -0.25, 0.24, -0.05) # fixed effects coefficients
sigma <- 0.5 # errors' standard deviation
D11 <- 1.0 # variance of random intercepts
D22 <- 0.5 # variance of random slopes
# we simulate random effects
b <- cbind(rnorm(n, sd = sqrt(D11)), rnorm(n, sd = sqrt(D22)))
# linear predictor
eta_y <- as.vector(X %*% betas + rowSums(Z * b[DF$id, ]))
# we simulate normal longitudinal data
DF$y <- rnorm(n * K, mean = eta_y, sd = sigma)
# we assume that values below -4 are not observed, and set equal to -4
DF$ind <- as.numeric(DF$y < -4)
DF$y <- pmax(DF$y, -4)
upp_Cens <- 15 # fixed Type I censoring time
shape_wb <- 5 # shape Weibull
alpha <- 0.8 # association coefficients
gammas <- c("(Intercept)" = -9, "sex" = 0.5)
W <- model.matrix(~ sex, data = DF[!duplicated(DF$id), ])
# linear predictor for the survival model
eta_t <- as.vector(W %*% gammas)
# to simulate event times we use inverse transform sampling
# (https://en.wikipedia.org/wiki/Inverse_transform_sampling). Namely, we want
# to find t, such that S(t) = u, where S(.) is the survival function, and u a
# number from the Unif(0, 1) distribution. The function below calculates
# log(u) - log(S(t)), and for a given u, we want to find t for which it equals
# zero. We do that below using the uniroot() function
invS <- function (t, i) {
# i denotes the subject
sex_i <- W[i, 2L]
# h() is the hazard function and we assume a Weibull baseline hazard
h <- function (s) {
X_at_s <- cbind(1, sex_i, s, sex_i * s)
Z_at_s <- cbind(1, s)
# the linear predictor from the mixed model evaluated at time s
f <- as.vector(X_at_s %*% betas +
rowSums(Z_at_s * b[rep(i, nrow(Z_at_s)), ]))
exp(log(shape_wb) + (shape_wb - 1) * log(s) + eta_t[i] + f * alpha)
}
# -log(S(t)) = H(t), where H(t) is the cumulative hazard function
integrate(h, lower = 0, upper = t)$value + log(u[i])
}
# we simulate the event times
u <- runif(n)
trueTimes <- numeric(n)
for (i in seq_len(n)) {
Up <- 100
Root <- try(uniroot(invS, interval = c(1e-05, Up), i = i)$root, TRUE)
trueTimes[i] <- if (!inherits(Root, "try-error")) Root else 150
}
# we use fixed Type I right censoring denoting the end of the trial.
Ctimes <- upp_Cens
Time <- pmin(trueTimes, Ctimes)
event <- as.numeric(trueTimes <= Ctimes) # event indicator
# we keep the longitudinal measurements before the event times
DF$Time <- Time[DF$id]
DF$event <- event[DF$id]
DF <- DF[DF$time <= DF$Time, ]
The corresponding joint model is fitted with the following syntax:
DF_id <- DF[!duplicated(DF$id), ]
Cox_fit <- coxph(Surv(Time, event) ~ sex, data = DF_id)
CensNorm_MixMod <-
mixed_model(cbind(y, ind) ~ sex * time, random = ~ time | id, data = DF,
family = censored.normal())
jointFit <- jm(Cox_fit, CensNorm_MixMod, time_var = "time")
summary(jointFit)
Students's-t mixed models
Outlying observations are a common issue in practice. Several methods have been proposed in the literature for identifying such observations in the context of longitudinal data. However, removing such values from the analysis is generally not recommended unless we also have external information as to why these values are outlying. Hence, we would need to fit mixed models to accommodate these observations in these settings. A well-known approach to achieve this is replacing the normal distribution for the error terms in the linear mixed model with a Student's-t distribution with heavier tails.
The following syntax simulates data from a joint model with a Student's-t mixed effects model:
set.seed(1234)
n <- 200 # number of subjects
K <- 12 # number of measurements per subject
t_max <- 14 # maximum follow-up time
# we construct a data frame with the design:
# everyone has a baseline measurement, and then measurements at random
# follow-up times up to t_max
DF <- data.frame(id = rep(seq_len(n), each = K),
time = c(replicate(n, c(0, sort(runif(K - 1, 0, t_max))))),
sex = rep(gl(2, n/2, labels = c("male", "female")), each = K))
# design matrices for the fixed and random effects
X <- model.matrix(~ sex * time, data = DF)
Z <- model.matrix(~ time, data = DF)
betas <- c(-2.2, -0.25, 0.24, -0.05) # fixed effects coefficients
sigma <- 0.5 # error standard deviation
D11 <- 1.0 # variance of random intercepts
D22 <- 0.5 # variance of random slopes
# we simulate random effects
b <- cbind(rnorm(n, sd = sqrt(D11)), rnorm(n, sd = sqrt(D22)))
# linear predictor
eta_y <- as.vector(X %*% betas + rowSums(Z * b[DF$id, ]))
# we simulate Student's-t longitudinal data
DF$y <- eta_y + sigma * rt(n * K, df = 4)
upp_Cens <- 15 # fixed Type I censoring time
shape_wb <- 5 # shape Weibull
alpha <- 0.8 # association coefficients
gammas <- c("(Intercept)" = -9, "sex" = 0.5)
W <- model.matrix(~ sex, data = DF[!duplicated(DF$id), ])
# linear predictor for the survival model
eta_t <- as.vector(W %*% gammas)
# to simulate event times we use inverse transform sampling
# (https://en.wikipedia.org/wiki/Inverse_transform_sampling). Namely, we want
# to find t, such that S(t) = u, where S(.) is the survival function, and u a
# number from the Unif(0, 1) distribution. The function below calculates
# log(u) - log(S(t)), and for a given u, we want to find t for which it equals
# zero. We do that below using the uniroot() function
invS <- function (t, i) {
# i denotes the subject
sex_i <- W[i, 2L]
# h() is the hazard function and we assume a Weibull baseline hazard
h <- function (s) {
X_at_s <- cbind(1, sex_i, s, sex_i * s)
Z_at_s <- cbind(1, s)
# the linear predictor from the mixed model evaluated at time s
f <- as.vector(X_at_s %*% betas +
rowSums(Z_at_s * b[rep(i, nrow(Z_at_s)), ]))
exp(log(shape_wb) + (shape_wb - 1) * log(s) + eta_t[i] + f * alpha)
}
# -log(S(t)) = H(t), where H(t) is the cumulative hazard function
integrate(h, lower = 0, upper = t)$value + log(u[i])
}
# we simulate the event times
u <- runif(n)
trueTimes <- numeric(n)
for (i in seq_len(n)) {
Up <- 100
Root <- try(uniroot(invS, interval = c(1e-05, Up), i = i)$root, TRUE)
trueTimes[i] <- if (!inherits(Root, "try-error")) Root else 150
}
# we use fixed Type I right censoring denoting the end of the trial.
Ctimes <- upp_Cens
Time <- pmin(trueTimes, Ctimes)
event <- as.numeric(trueTimes <= Ctimes) # event indicator
# we keep the longitudinal measurements before the event times
DF$Time <- Time[DF$id]
DF$event <- event[DF$id]
DF <- DF[DF$time <= DF$Time, ]
To fit the corresponding joint model we use the students.t() family object in the call to mixed_model():
DF_id <- DF[!duplicated(DF$id), ]
Cox_fit <- coxph(Surv(Time, event) ~ sex, data = DF_id)
Stdt_MixMod <-
mixed_model(y ~ sex * time, random = ~ time | id, data = DF,
family = students.t(df = 4))
jointFit <- jm(Cox_fit, Stdt_MixMod, time_var = "time")
summary(jointFit)
Negative binomial mixed models
Count longitudinal outcomes are typically modeled with the Poisson distribution. However, these outcomes often exhibit more variance than what is allowed from the Poisson distribution, leading to the well-known problem of over-dispersion. To accommodate this over-dispersion, typically, the negative binomial distribution is used.
The following piece of code simulates data from a joint model for count longitudinal data that follow the negative binomial distribution:
set.seed(1234)
n <- 500 # number of subjects
K <- 10 # number of measurements per subject
t_max <- 5 # maximum follow-up time
# we construct a data frame with the design:
# everyone has a baseline measurement, and then measurements at random
# follow-up times up to t_max
DF <- data.frame(id = rep(seq_len(n), each = K),
time = c(replicate(n, c(0, sort(runif(K - 1, 0, t_max))))),
sex = rep(gl(2, n/2, labels = c("male", "female")), each = K))
# design matrices for the fixed and random effects
X <- model.matrix(~ sex * time, data = DF)
Z <- model.matrix(~ time, data = DF)
betas <- c(0.8, -0.5, 0.8, -0.5) # fixed effects coefficients
shape <- 2 # shape/size parameter of the negative binomial distribution
D11 <- 1.0 # variance of random intercepts
D22 <- 0.3 # variance of random slopes
# we simulate random effects
b <- cbind(rnorm(n, sd = sqrt(D11)), rnorm(n, sd = sqrt(D22)))
# linear predictor
eta_y <- as.vector(X %*% betas + rowSums(Z * b[DF$id, ]))
# mean of the Beta distribution
mu_y <- plogis(eta_y) # plogis(eta_y) = exp(eta_y) / (1 + exp(eta_y))
# we simulate negative binomial longitudinal data
DF$y <- rnbinom(n * K, size = shape, mu = exp(eta_y))
# simulate event times
upp_Cens <- 5 # fixed Type I censoring time
shape_wb <- 5 # shape Weibull
alpha <- 0.8 # association coefficient
gammas <- c("(Intercept)" = -9, "sex" = 0.5)
W <- model.matrix(~ sex, data = DF[!duplicated(DF$id), ])
# linear predictor for the survival model
eta_t <- as.vector(W %*% gammas)
# to simulate event times we use inverse transform sampling
# (https://en.wikipedia.org/wiki/Inverse_transform_sampling). Namely, we want
# to find t, such that S(t) = u, where S(.) is the survival function, and u a
# number from the Unif(0, 1) distribution. The function below calculates
# log(u) - log(S(t)), and for a given u, we want to find t for which it equals
# zero. We do that below using the uniroot() function
invS <- function (t, i) {
# i denotes the subject
sex_i <- W[i, 2L]
# h() is the hazard function and we assume a Weibull baseline hazard
h <- function (s) {
X_at_s <- cbind(1, sex_i, s, sex_i * s)
Z_at_s <- cbind(1, s)
# the linear predictor from the mixed model evaluated at time s
f <- as.vector(X_at_s %*% betas +
rowSums(Z_at_s * b[rep(i, nrow(Z_at_s)), ]))
exp(log(shape_wb) + (shape_wb - 1) * log(s) + eta_t[i] + f * alpha)
}
# -log(S(t)) = H(t), where H(t) is the cumulative hazard function
integrate(h, lower = 0, upper = t)$value + log(u[i])
}
# we simulate the event times
u <- runif(n)
trueTimes <- numeric(n)
for (i in seq_len(n)) {
Up <- 100
Root <- try(uniroot(invS, interval = c(1e-05, Up), i = i)$root, TRUE)
trueTimes[i] <- if (!inherits(Root, "try-error")) Root else 150
}
# we use fixed Type I right censoring denoting the end of the trial.
Ctimes <- upp_Cens
Time <- pmin(trueTimes, Ctimes)
event <- as.numeric(trueTimes <= Ctimes) # event indicator
# we keep the longitudinal measurements before the event times
DF$Time <- Time[DF$id]
DF$event <- event[DF$id]
DF <- DF[DF$time <= DF$Time, ]
The corresponding joint model is the fitted using the following syntax:
DF_id <- DF[!duplicated(DF$id), ]
Cox_fit <- coxph(Surv(Time, event) ~ sex, data = DF_id)
NB_MixMod <- mixed_model(y ~ sex * time, random = ~ time | id, data = DF,
family = negative.binomial())
jointFit <- jm(Cox_fit, NB_MixMod, time_var = "time")
summary(jointFit)
Beta-binomial longitudinal outcomes
For count data and binomial data, we may have an over-dispersion problem. To accommodate this, we can change the standard binomial distribution to a beta-binomial one.
The following piece of code simulates data from a joint model for binomial longitudinal data that follow the beta-binomial distribution:
set.seed(1234)
n <- 500 # number of subjects
K <- 8 # number of measurements per subject
t_max <- 10 # maximum follow-up time
# we construct a data frame with the design:
# everyone has a baseline measurement, and then measurements at random
# follow-up times up to t_max
DF <- data.frame(id = rep(seq_len(n), each = K),
time = c(replicate(n, c(0, sort(runif(K - 1, 0, t_max))))),
sex = rep(gl(2, n/2, labels = c("male", "female")), each = K))
# design matrices for the fixed and random effects
X <- model.matrix(~ sex * time, data = DF)
Z <- model.matrix(~ time, data = DF)
betas <- c(-2.2, -0.25, 0.24, -0.05) # fixed effects coefficients
phi <- 5 # precision parameter of the Beta distribution
D11 <- 1.0 # variance of random intercepts
D22 <- 0.5 # variance of random slopes
# we simulate random effects
b <- cbind(rnorm(n, sd = sqrt(D11)), rnorm(n, sd = sqrt(D22)))
# linear predictor
eta_y <- as.vector(X %*% betas + rowSums(Z * b[DF$id, ]))
# mean of the Beta distribution
mu_y <- plogis(eta_y) # plogis(eta_y) = exp(eta_y) / (1 + exp(eta_y))
# we simulate probabilities from the Beta distribution
probs <- rbeta(n * K, shape1 = mu_y * phi, shape2 = phi * (1 - mu_y))
# we transform to (0, 1)
probs <- (probs * (nrow(DF) - 1) + 0.5) / nrow(DF)
# we simulate binomial data use the probs
DF$y <- rbinom(n * K, size = 20, prob = probs)
upp_Cens <- 15 # fixed Type I censoring time
shape_wb <- 5 # shape Weibull
alpha <- 0.8 # association coefficients
gammas <- c("(Intercept)" = -9, "sex" = 0.5)
W <- model.matrix(~ sex, data = DF[!duplicated(DF$id), ])
# linear predictor for the survival model
eta_t <- as.vector(W %*% gammas)
# to simulate event times we use inverse transform sampling
# (https://en.wikipedia.org/wiki/Inverse_transform_sampling). Namely, we want
# to find t, such that S(t) = u, where S(.) is the survival function, and u a
# number from the Unif(0, 1) distribution. The function below calculates
# log(u) - log(S(t)), and for a given u, we want to find t for which it equals
# zero. We do that below using the uniroot() function
invS <- function (t, i) {
# i denotes the subject
sex_i <- W[i, 2L]
# h() is the hazard function and we assume a Weibull baseline hazard
h <- function (s) {
X_at_s <- cbind(1, sex_i, s, sex_i * s)
Z_at_s <- cbind(1, s)
# the linear predictor from the mixed model evaluated at time s
f <- as.vector(X_at_s %*% betas +
rowSums(Z_at_s * b[rep(i, nrow(Z_at_s)), ]))
exp(log(shape_wb) + (shape_wb - 1) * log(s) + eta_t[i] + f * alpha)
}
# -log(S(t)) = H(t), where H(t) is the cumulative hazard function
integrate(h, lower = 0, upper = t)$value + log(u[i])
}
# we simulate the event times
u <- runif(n)
trueTimes <- numeric(n)
for (i in seq_len(n)) {
Up <- 100
Root <- try(uniroot(invS, interval = c(1e-05, Up), i = i)$root, TRUE)
trueTimes[i] <- if (!inherits(Root, "try-error")) Root else 150
}
# we use fixed Type I right censoring denoting the end of the trial.
Ctimes <- upp_Cens
Time <- pmin(trueTimes, Ctimes)
event <- as.numeric(trueTimes <= Ctimes) # event indicator
# we keep the longitudinal measurements before the event times
DF$Time <- Time[DF$id]
DF$event <- event[DF$id]
DF <- DF[DF$time <= DF$Time, ]
The corresponding joint model is then fitted with the syntax:
DF_id <- DF[!duplicated(DF$id), ]
Cox_fit <- coxph(Surv(Time, event) ~ sex, data = DF_id)
BetaBinom_MixMod <-
mixed_model(cbind(y, 20 - y) ~ sex * time, random = ~ time | id, data = DF,
family = beta.binomial())
jointFit <- jm(Cox_fit, BetaBinom_MixMod, time_var = "time")
summary(jointFit)
drizopoulos/JMbayes2 documentation built on July 15, 2024, 11:13 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) library("JMbayes2")
Non-Gaussian Joint Models with JMbayes2
Taking advantage of the versatility of the GLMMadaptive package, JMbayes2 can fit joint models with several different types of mixed-effects models. The following examples illustrate these capabilities. All examples have the same structure, namely, first, a short motivation for each mixed-model is given, followed by a piece of R code simulating data from a joint model with the respective mixed-effects sub-model, closing by the syntax to fit the joint model. In this last part, the main difference per example is the call to mixed_model().
Beta mixed models
With very few exceptions, continuous outcomes that we wish to analyze have some natural bounds. For example, the levels of blood biomarkers for a set of patients. However, most observations are often located far away from these natural bounds, and an assumption of a normal distribution for the outcome can be safely made. In some settings, though, we can have outcomes for which a substantial percentage of the observations are located near the boundaries leading to skewed or U-shaped distributions. A linear mixed model with normal error terms often does not provide a good fit for such longitudinal outcomes. A natural alternative is to select a distribution that respects the bounded nature of the outcome. The most well-known distribution for such outcomes is the Beta distribution defined in the $(0, 1)$ interval (note: a bounded outcome $Y^$ in the $(a, b)$ interval can be transformed to the $Y = (Y^ - a) / (b - a)$ in the $(0, 1)$ interval).
The following code illustrates how to simulate data from a joint model with a Beta mixed effects model. The default functional form is assumed, i.e., that the linear predictor $\eta(t)$ of the mixed model is associated with the hazard of an event at time $t$. The linear predictor is related to the mean $\mu(t)$ of the Beta distribution under the logit link function, i.e., $\log[\mu(t) / {1 - \mu(t)}] = \eta(t)$.
set.seed(1234) n <- 200 # number of subjects K <- 8 # number of measurements per subject t_max <- 10 # maximum follow-up time # we construct a data frame with the design: # everyone has a baseline measurement, and then measurements at random # follow-up times up to t_max DF <- data.frame(id = rep(seq_len(n), each = K), time = c(replicate(n, c(0, sort(runif(K - 1, 0, t_max))))), sex = rep(gl(2, n/2, labels = c("male", "female")), each = K)) # design matrices for the fixed and random effects X <- model.matrix(~ sex * time, data = DF) Z <- model.matrix(~ time, data = DF) betas <- c(-2.2, -0.25, 0.24, -0.05) # fixed effects coefficients phi <- 5 # precision parameter of the Beta distribution D11 <- 1.0 # variance of random intercepts D22 <- 0.5 # variance of random slopes # we simulate random effects b <- cbind(rnorm(n, sd = sqrt(D11)), rnorm(n, sd = sqrt(D22))) # linear predictor eta_y <- as.vector(X %*% betas + rowSums(Z * b[DF$id, ])) # mean of the Beta distribution mu_y <- plogis(eta_y) # plogis(eta_y) = exp(eta_y) / (1 + exp(eta_y)) # we simulate Beta longitudinal data DF$y <- rbeta(n * K, shape1 = mu_y * phi, shape2 = phi * (1 - mu_y)) # we transform to (0, 1) DF$y <- (DF$y * (nrow(DF) - 1) + 0.5) / nrow(DF) upp_Cens <- 15 # fixed Type I censoring time shape_wb <- 5 # shape Weibull alpha <- 0.8 # association coefficients gammas <- c("(Intercept)" = -9, "sex" = 0.5) W <- model.matrix(~ sex, data = DF[!duplicated(DF$id), ]) # linear predictor for the survival model eta_t <- as.vector(W %*% gammas) # to simulate event times we use inverse transform sampling # (https://en.wikipedia.org/wiki/Inverse_transform_sampling). Namely, we want # to find t, such that S(t) = u, where S(.) is the survival function, and u a # number from the Unif(0, 1) distribution. The function below calculates # log(u) - log(S(t)), and for a given u, we want to find t for which it equals # zero. We do that below using the uniroot() function invS <- function (t, i) { # i denotes the subject sex_i <- W[i, 2L] # h() is the hazard function and we assume a Weibull baseline hazard h <- function (s) { X_at_s <- cbind(1, sex_i, s, sex_i * s) Z_at_s <- cbind(1, s) # the linear predictor from the mixed model evaluated at time s f <- as.vector(X_at_s %*% betas + rowSums(Z_at_s * b[rep(i, nrow(Z_at_s)), ])) exp(log(shape_wb) + (shape_wb - 1) * log(s) + eta_t[i] + f * alpha) } # -log(S(t)) = H(t), where H(t) is the cumulative hazard function integrate(h, lower = 0, upper = t)$value + log(u[i]) } # we simulate the event times u <- runif(n) trueTimes <- numeric(n) for (i in seq_len(n)) { Up <- 100 Root <- try(uniroot(invS, interval = c(1e-05, Up), i = i)$root, TRUE) trueTimes[i] <- if (!inherits(Root, "try-error")) Root else 150 } # we use fixed Type I right censoring denoting the end of the trial. Ctimes <- upp_Cens Time <- pmin(trueTimes, Ctimes) event <- as.numeric(trueTimes <= Ctimes) # event indicator # we keep the longitudinal measurements before the event times DF$Time <- Time[DF$id] DF$event <- event[DF$id] DF <- DF[DF$time <= DF$Time, ]
To fit the corresponding joint model, we fit first a Beta mixed model using the beta.fam() family object into the call of mixed_model():
DF_id <- DF[!duplicated(DF$id), ] Cox_fit <- coxph(Surv(Time, event) ~ sex, data = DF_id) Beta_MixMod <- mixed_model(y ~ sex * time, random = ~ time | id, data = DF, family = beta.fam()) jointFit <- jm(Cox_fit, Beta_MixMod, time_var = "time") summary(jointFit)
Censored linear mixed models
Some continuous longitudinal outcomes may have a censored nature. A typical example of such outcomes is when we have limit of detection issues. That is, the values of the outcome cannot be detected below a specified threshold having to do with the (laboratory) equipment used to determine the measurements. In these settings, and even if the complete data follows a normal distribution the observed censored data cannot be analyzed with a standard mixed model. The mixed_model() function can accommodate such outcomes using the censored.normal() family object.
The following piece of code simulates data from a joint model with a linear mixed model for the longitudinal outcomes but applies censoring in the realized longitudinal observations.
set.seed(1234) n <- 200 # number of subjects K <- 12 # number of measurements per subject t_max <- 14 # maximum follow-up time # we construct a data frame with the design: # everyone has a baseline measurement, and then measurements at random # follow-up times up to t_max DF <- data.frame(id = rep(seq_len(n), each = K), time = c(replicate(n, c(0, sort(runif(K - 1, 0, t_max))))), sex = rep(gl(2, n/2, labels = c("male", "female")), each = K)) # design matrices for the fixed and random effects X <- model.matrix(~ sex * time, data = DF) Z <- model.matrix(~ time, data = DF) betas <- c(-2.2, -0.25, 0.24, -0.05) # fixed effects coefficients sigma <- 0.5 # errors' standard deviation D11 <- 1.0 # variance of random intercepts D22 <- 0.5 # variance of random slopes # we simulate random effects b <- cbind(rnorm(n, sd = sqrt(D11)), rnorm(n, sd = sqrt(D22))) # linear predictor eta_y <- as.vector(X %*% betas + rowSums(Z * b[DF$id, ])) # we simulate normal longitudinal data DF$y <- rnorm(n * K, mean = eta_y, sd = sigma) # we assume that values below -4 are not observed, and set equal to -4 DF$ind <- as.numeric(DF$y < -4) DF$y <- pmax(DF$y, -4) upp_Cens <- 15 # fixed Type I censoring time shape_wb <- 5 # shape Weibull alpha <- 0.8 # association coefficients gammas <- c("(Intercept)" = -9, "sex" = 0.5) W <- model.matrix(~ sex, data = DF[!duplicated(DF$id), ]) # linear predictor for the survival model eta_t <- as.vector(W %*% gammas) # to simulate event times we use inverse transform sampling # (https://en.wikipedia.org/wiki/Inverse_transform_sampling). Namely, we want # to find t, such that S(t) = u, where S(.) is the survival function, and u a # number from the Unif(0, 1) distribution. The function below calculates # log(u) - log(S(t)), and for a given u, we want to find t for which it equals # zero. We do that below using the uniroot() function invS <- function (t, i) { # i denotes the subject sex_i <- W[i, 2L] # h() is the hazard function and we assume a Weibull baseline hazard h <- function (s) { X_at_s <- cbind(1, sex_i, s, sex_i * s) Z_at_s <- cbind(1, s) # the linear predictor from the mixed model evaluated at time s f <- as.vector(X_at_s %*% betas + rowSums(Z_at_s * b[rep(i, nrow(Z_at_s)), ])) exp(log(shape_wb) + (shape_wb - 1) * log(s) + eta_t[i] + f * alpha) } # -log(S(t)) = H(t), where H(t) is the cumulative hazard function integrate(h, lower = 0, upper = t)$value + log(u[i]) } # we simulate the event times u <- runif(n) trueTimes <- numeric(n) for (i in seq_len(n)) { Up <- 100 Root <- try(uniroot(invS, interval = c(1e-05, Up), i = i)$root, TRUE) trueTimes[i] <- if (!inherits(Root, "try-error")) Root else 150 } # we use fixed Type I right censoring denoting the end of the trial. Ctimes <- upp_Cens Time <- pmin(trueTimes, Ctimes) event <- as.numeric(trueTimes <= Ctimes) # event indicator # we keep the longitudinal measurements before the event times DF$Time <- Time[DF$id] DF$event <- event[DF$id] DF <- DF[DF$time <= DF$Time, ]
The corresponding joint model is fitted with the following syntax:
DF_id <- DF[!duplicated(DF$id), ] Cox_fit <- coxph(Surv(Time, event) ~ sex, data = DF_id) CensNorm_MixMod <- mixed_model(cbind(y, ind) ~ sex * time, random = ~ time | id, data = DF, family = censored.normal()) jointFit <- jm(Cox_fit, CensNorm_MixMod, time_var = "time") summary(jointFit)
Students's-t mixed models
Outlying observations are a common issue in practice. Several methods have been proposed in the literature for identifying such observations in the context of longitudinal data. However, removing such values from the analysis is generally not recommended unless we also have external information as to why these values are outlying. Hence, we would need to fit mixed models to accommodate these observations in these settings. A well-known approach to achieve this is replacing the normal distribution for the error terms in the linear mixed model with a Student's-t distribution with heavier tails.
The following syntax simulates data from a joint model with a Student's-t mixed effects model:
set.seed(1234) n <- 200 # number of subjects K <- 12 # number of measurements per subject t_max <- 14 # maximum follow-up time # we construct a data frame with the design: # everyone has a baseline measurement, and then measurements at random # follow-up times up to t_max DF <- data.frame(id = rep(seq_len(n), each = K), time = c(replicate(n, c(0, sort(runif(K - 1, 0, t_max))))), sex = rep(gl(2, n/2, labels = c("male", "female")), each = K)) # design matrices for the fixed and random effects X <- model.matrix(~ sex * time, data = DF) Z <- model.matrix(~ time, data = DF) betas <- c(-2.2, -0.25, 0.24, -0.05) # fixed effects coefficients sigma <- 0.5 # error standard deviation D11 <- 1.0 # variance of random intercepts D22 <- 0.5 # variance of random slopes # we simulate random effects b <- cbind(rnorm(n, sd = sqrt(D11)), rnorm(n, sd = sqrt(D22))) # linear predictor eta_y <- as.vector(X %*% betas + rowSums(Z * b[DF$id, ])) # we simulate Student's-t longitudinal data DF$y <- eta_y + sigma * rt(n * K, df = 4) upp_Cens <- 15 # fixed Type I censoring time shape_wb <- 5 # shape Weibull alpha <- 0.8 # association coefficients gammas <- c("(Intercept)" = -9, "sex" = 0.5) W <- model.matrix(~ sex, data = DF[!duplicated(DF$id), ]) # linear predictor for the survival model eta_t <- as.vector(W %*% gammas) # to simulate event times we use inverse transform sampling # (https://en.wikipedia.org/wiki/Inverse_transform_sampling). Namely, we want # to find t, such that S(t) = u, where S(.) is the survival function, and u a # number from the Unif(0, 1) distribution. The function below calculates # log(u) - log(S(t)), and for a given u, we want to find t for which it equals # zero. We do that below using the uniroot() function invS <- function (t, i) { # i denotes the subject sex_i <- W[i, 2L] # h() is the hazard function and we assume a Weibull baseline hazard h <- function (s) { X_at_s <- cbind(1, sex_i, s, sex_i * s) Z_at_s <- cbind(1, s) # the linear predictor from the mixed model evaluated at time s f <- as.vector(X_at_s %*% betas + rowSums(Z_at_s * b[rep(i, nrow(Z_at_s)), ])) exp(log(shape_wb) + (shape_wb - 1) * log(s) + eta_t[i] + f * alpha) } # -log(S(t)) = H(t), where H(t) is the cumulative hazard function integrate(h, lower = 0, upper = t)$value + log(u[i]) } # we simulate the event times u <- runif(n) trueTimes <- numeric(n) for (i in seq_len(n)) { Up <- 100 Root <- try(uniroot(invS, interval = c(1e-05, Up), i = i)$root, TRUE) trueTimes[i] <- if (!inherits(Root, "try-error")) Root else 150 } # we use fixed Type I right censoring denoting the end of the trial. Ctimes <- upp_Cens Time <- pmin(trueTimes, Ctimes) event <- as.numeric(trueTimes <= Ctimes) # event indicator # we keep the longitudinal measurements before the event times DF$Time <- Time[DF$id] DF$event <- event[DF$id] DF <- DF[DF$time <= DF$Time, ]
To fit the corresponding joint model we use the students.t() family object in the call to mixed_model():
DF_id <- DF[!duplicated(DF$id), ] Cox_fit <- coxph(Surv(Time, event) ~ sex, data = DF_id) Stdt_MixMod <- mixed_model(y ~ sex * time, random = ~ time | id, data = DF, family = students.t(df = 4)) jointFit <- jm(Cox_fit, Stdt_MixMod, time_var = "time") summary(jointFit)
Negative binomial mixed models
Count longitudinal outcomes are typically modeled with the Poisson distribution. However, these outcomes often exhibit more variance than what is allowed from the Poisson distribution, leading to the well-known problem of over-dispersion. To accommodate this over-dispersion, typically, the negative binomial distribution is used.
The following piece of code simulates data from a joint model for count longitudinal data that follow the negative binomial distribution:
set.seed(1234) n <- 500 # number of subjects K <- 10 # number of measurements per subject t_max <- 5 # maximum follow-up time # we construct a data frame with the design: # everyone has a baseline measurement, and then measurements at random # follow-up times up to t_max DF <- data.frame(id = rep(seq_len(n), each = K), time = c(replicate(n, c(0, sort(runif(K - 1, 0, t_max))))), sex = rep(gl(2, n/2, labels = c("male", "female")), each = K)) # design matrices for the fixed and random effects X <- model.matrix(~ sex * time, data = DF) Z <- model.matrix(~ time, data = DF) betas <- c(0.8, -0.5, 0.8, -0.5) # fixed effects coefficients shape <- 2 # shape/size parameter of the negative binomial distribution D11 <- 1.0 # variance of random intercepts D22 <- 0.3 # variance of random slopes # we simulate random effects b <- cbind(rnorm(n, sd = sqrt(D11)), rnorm(n, sd = sqrt(D22))) # linear predictor eta_y <- as.vector(X %*% betas + rowSums(Z * b[DF$id, ])) # mean of the Beta distribution mu_y <- plogis(eta_y) # plogis(eta_y) = exp(eta_y) / (1 + exp(eta_y)) # we simulate negative binomial longitudinal data DF$y <- rnbinom(n * K, size = shape, mu = exp(eta_y)) # simulate event times upp_Cens <- 5 # fixed Type I censoring time shape_wb <- 5 # shape Weibull alpha <- 0.8 # association coefficient gammas <- c("(Intercept)" = -9, "sex" = 0.5) W <- model.matrix(~ sex, data = DF[!duplicated(DF$id), ]) # linear predictor for the survival model eta_t <- as.vector(W %*% gammas) # to simulate event times we use inverse transform sampling # (https://en.wikipedia.org/wiki/Inverse_transform_sampling). Namely, we want # to find t, such that S(t) = u, where S(.) is the survival function, and u a # number from the Unif(0, 1) distribution. The function below calculates # log(u) - log(S(t)), and for a given u, we want to find t for which it equals # zero. We do that below using the uniroot() function invS <- function (t, i) { # i denotes the subject sex_i <- W[i, 2L] # h() is the hazard function and we assume a Weibull baseline hazard h <- function (s) { X_at_s <- cbind(1, sex_i, s, sex_i * s) Z_at_s <- cbind(1, s) # the linear predictor from the mixed model evaluated at time s f <- as.vector(X_at_s %*% betas + rowSums(Z_at_s * b[rep(i, nrow(Z_at_s)), ])) exp(log(shape_wb) + (shape_wb - 1) * log(s) + eta_t[i] + f * alpha) } # -log(S(t)) = H(t), where H(t) is the cumulative hazard function integrate(h, lower = 0, upper = t)$value + log(u[i]) } # we simulate the event times u <- runif(n) trueTimes <- numeric(n) for (i in seq_len(n)) { Up <- 100 Root <- try(uniroot(invS, interval = c(1e-05, Up), i = i)$root, TRUE) trueTimes[i] <- if (!inherits(Root, "try-error")) Root else 150 } # we use fixed Type I right censoring denoting the end of the trial. Ctimes <- upp_Cens Time <- pmin(trueTimes, Ctimes) event <- as.numeric(trueTimes <= Ctimes) # event indicator # we keep the longitudinal measurements before the event times DF$Time <- Time[DF$id] DF$event <- event[DF$id] DF <- DF[DF$time <= DF$Time, ]
The corresponding joint model is the fitted using the following syntax:
DF_id <- DF[!duplicated(DF$id), ] Cox_fit <- coxph(Surv(Time, event) ~ sex, data = DF_id) NB_MixMod <- mixed_model(y ~ sex * time, random = ~ time | id, data = DF, family = negative.binomial()) jointFit <- jm(Cox_fit, NB_MixMod, time_var = "time") summary(jointFit)
Beta-binomial longitudinal outcomes
For count data and binomial data, we may have an over-dispersion problem. To accommodate this, we can change the standard binomial distribution to a beta-binomial one.
The following piece of code simulates data from a joint model for binomial longitudinal data that follow the beta-binomial distribution:
set.seed(1234) n <- 500 # number of subjects K <- 8 # number of measurements per subject t_max <- 10 # maximum follow-up time # we construct a data frame with the design: # everyone has a baseline measurement, and then measurements at random # follow-up times up to t_max DF <- data.frame(id = rep(seq_len(n), each = K), time = c(replicate(n, c(0, sort(runif(K - 1, 0, t_max))))), sex = rep(gl(2, n/2, labels = c("male", "female")), each = K)) # design matrices for the fixed and random effects X <- model.matrix(~ sex * time, data = DF) Z <- model.matrix(~ time, data = DF) betas <- c(-2.2, -0.25, 0.24, -0.05) # fixed effects coefficients phi <- 5 # precision parameter of the Beta distribution D11 <- 1.0 # variance of random intercepts D22 <- 0.5 # variance of random slopes # we simulate random effects b <- cbind(rnorm(n, sd = sqrt(D11)), rnorm(n, sd = sqrt(D22))) # linear predictor eta_y <- as.vector(X %*% betas + rowSums(Z * b[DF$id, ])) # mean of the Beta distribution mu_y <- plogis(eta_y) # plogis(eta_y) = exp(eta_y) / (1 + exp(eta_y)) # we simulate probabilities from the Beta distribution probs <- rbeta(n * K, shape1 = mu_y * phi, shape2 = phi * (1 - mu_y)) # we transform to (0, 1) probs <- (probs * (nrow(DF) - 1) + 0.5) / nrow(DF) # we simulate binomial data use the probs DF$y <- rbinom(n * K, size = 20, prob = probs) upp_Cens <- 15 # fixed Type I censoring time shape_wb <- 5 # shape Weibull alpha <- 0.8 # association coefficients gammas <- c("(Intercept)" = -9, "sex" = 0.5) W <- model.matrix(~ sex, data = DF[!duplicated(DF$id), ]) # linear predictor for the survival model eta_t <- as.vector(W %*% gammas) # to simulate event times we use inverse transform sampling # (https://en.wikipedia.org/wiki/Inverse_transform_sampling). Namely, we want # to find t, such that S(t) = u, where S(.) is the survival function, and u a # number from the Unif(0, 1) distribution. The function below calculates # log(u) - log(S(t)), and for a given u, we want to find t for which it equals # zero. We do that below using the uniroot() function invS <- function (t, i) { # i denotes the subject sex_i <- W[i, 2L] # h() is the hazard function and we assume a Weibull baseline hazard h <- function (s) { X_at_s <- cbind(1, sex_i, s, sex_i * s) Z_at_s <- cbind(1, s) # the linear predictor from the mixed model evaluated at time s f <- as.vector(X_at_s %*% betas + rowSums(Z_at_s * b[rep(i, nrow(Z_at_s)), ])) exp(log(shape_wb) + (shape_wb - 1) * log(s) + eta_t[i] + f * alpha) } # -log(S(t)) = H(t), where H(t) is the cumulative hazard function integrate(h, lower = 0, upper = t)$value + log(u[i]) } # we simulate the event times u <- runif(n) trueTimes <- numeric(n) for (i in seq_len(n)) { Up <- 100 Root <- try(uniroot(invS, interval = c(1e-05, Up), i = i)$root, TRUE) trueTimes[i] <- if (!inherits(Root, "try-error")) Root else 150 } # we use fixed Type I right censoring denoting the end of the trial. Ctimes <- upp_Cens Time <- pmin(trueTimes, Ctimes) event <- as.numeric(trueTimes <= Ctimes) # event indicator # we keep the longitudinal measurements before the event times DF$Time <- Time[DF$id] DF$event <- event[DF$id] DF <- DF[DF$time <= DF$Time, ]
The corresponding joint model is then fitted with the syntax:
DF_id <- DF[!duplicated(DF$id), ] Cox_fit <- coxph(Surv(Time, event) ~ sex, data = DF_id) BetaBinom_MixMod <- mixed_model(cbind(y, 20 - y) ~ sex * time, random = ~ time | id, data = DF, family = beta.binomial()) jointFit <- jm(Cox_fit, BetaBinom_MixMod, time_var = "time") summary(jointFit)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.