Joint Survival and Marker Models
This package provides means of estimating joint models for a mixture of different types of markers and different types of survival outcomes. The models are estimated with Gaussian variational approximations (GVAs). Description of the supported models and multiple examples are given in this document. A presentation of the package is available here which includes a simulation study to investigate bias of the parameter estimates and coverage of approximate likelihood ratio based confidence intervals.
A comparison with the JM, JMbayes, and joineRML is provided on Github at github.com/boennecd/VAJointSurv/tree/main/comparisons.
Installation
The package can be installed from Github by calling:
remotes::install_github("boennecd/VAJointSurv", build_vignettes = TRUE)The Model
We will start by covering the model, then cover some examples, and end with details about the implementation.
The Markers
We assume that there are \(L\) observed markers at each time point. We let \(\vec Y_{ij}\) denote the \(j\)th observed markers for individual \(i\) and let \(s_{ij}\) be the observation time. The model for the markers is
\[\begin{align*} \vec Y_{ij}&= \vec\mu_i(s_{ij}, \vec U_i) + \vec\epsilon_{ij} \\ \epsilon_{ij} &\sim N^{(L)}(\vec 0, \Sigma) \\ \mu_{i1}(s, \vec U_i) &= \vec x_{i1}^\top\vec\gamma_1 + \vec g_1(s)^\top\vec\beta_1 + \vec m_1(s)^\top\vec U_{i1} \\ \vdots &\hphantom{=}\vdots\\\ \mu_{iL}(s, \vec U_i) &= \vec x_{iL}^\top\vec\gamma_L + \vec g_L(s)^\top\vec\beta_L + \vec m_L(s)^\top\vec U_{iL} \\ \vec U_i &= \begin{pmatrix} \vec U_{i1} \\ \vdots \\ \vec U_{iL} \end{pmatrix}\sim N^{(R)}(\vec0, \Psi) \end{align*}\]
where \(\vec\epsilon_{ij}\) is an error term of the observation, \(\vec U_i\) is the \(R\) dimensional unobserved random effect for individual \(i\), and the \(\vec g_k\)s and \(\vec m_k\)s are basis expansions. We can write the model more compactly by letting
\[G(s) = \begin{pmatrix} \vec g_1(s)^\top & 0^\top & \cdots & \vec 0^\top \\ \vec 0^\top & \vec g_2(s)^\top & \ddots & \vdots \\ \vdots & \ddots & \ddots & \vec 0^\top \\ \vec 0^\top & \cdots & \vec 0^\top & \vec g_L(s)^\top \end{pmatrix}\]
and defining \(M(s)\) and \(X_i\) similarly. Furthermore, let \(\vec\gamma = (\vec\gamma_1^\top,\dots,\vec\gamma_L^\top)\) and define \(\vec\beta\) similarly. Then the conditional mean vector at time \(s_{ij}\) is
\[\vec\mu_i(s_{ij}, \vec U_i) = X_i\vec\gamma + G(s_{ij})\vec\beta + M(s_{ij})\vec U_i.\]
The Survival Outcomes
We assume that there are \(H\) different types of survival outcomes. The conditional hazards of the survival outcomes are
\[\begin{align*} h_{i1}(t\mid \vec U_i, \vec\xi_i) &= \exp\left( \vec z_{i1}^\top\vec \delta_1 + \omega_1^\top \vec b_1(t) + \vec\alpha_1^\top M(t)\vec U_i + \xi_{i1} \right) \\ \vdots &\hphantom{=}\vdots \\ h_{iH}(t\mid \vec U_i,\vec \xi_i) &= \exp\left( \vec z_{iL}^\top\vec \delta_H + \omega_H^\top\vec b_H(t) + \vec\alpha_H^\top M(t)\vec U_i + \xi_{iH} \right) \\ \vec\xi_i = \begin{pmatrix}\xi_{i1} \\ \vdots \\ \xi_{iH}\end{pmatrix} &\sim N^{(H)}(\vec 0, \Xi). \end{align*}\]
The \(\vec\alpha\)s are association parameters which makes the markers and the survival outcomes marginally dependent. The \(\exp\xi_{ih}\)s are frailty effects which makes the survival outcomes marginally dependent even if \(\vec\alpha = \vec 0\). By default, these are not included as it is assumed that events are often terminal in which case identifying the frailty may be hard and evidence of a frailty is both evidence of non-proportional hazards and heterogeneity.
The observation process, the \(s_{ij}\)s, can be modeled as one of \(H\) different types of survival process. This can be done by adding each \(s_{ij}\) as a left-truncated outcome of the given type of process except for \(s_{i1}\) which is not left-truncated.
More General Survival Sub-model
The conditional hazards shown above only allows the hazard to be associated with the population deviation for the mean marker through the current value
\[M(t)\vec U_i = \begin{pmatrix} \vec m_1(t)^\top\vec U_{i1} \\ \vdots \\ \vec m_L(t)^\top\vec U_{i1} \end{pmatrix}.\]
However, the researcher may be interested in using the cumulative value, \(\int^t \vec m_k(s)^\top ds\vec U_{ik}\), the derivative, \(\partial\vec m_k(t)^\top\vec U_{ik}/\partial t\), or a combination of these and the current value. This is supported in the package through the ders argument of surv_term. As an example, if we let \(\vec d_k(t) = \int^t \vec m_k(s) ds\), \(\vec r_k(t) = \partial\vec m_k(t)/\partial t\), and \(L = 3\) then it is possible to replace the \(\vec\alpha_1^\top M(t)\vec U_i\) in the hazard with
\[\vec\alpha_1^\top\begin{pmatrix} \vec d_1(t)^\top \vec U_{i1} \\ \vec m_1(t)^\top \vec U_{i1} \\ \vec m_2(t)^\top \vec U_{i2} \\ \vec d_3(t)^\top \vec U_{i3} \\ \vec r_3(t)^\top \vec U_{i3} \\ \end{pmatrix}\]
That is, a cumulative and current value for the first marker, only the current value for the second marker, and the cumulative and derivative for the third marker.
Model Support
The package can handle
- Mixed type of basis expansions in \(G\), \(M\), and \(b_1,\dots,b_H\).
- A maximum \(L = 31\) different markers, and an arbitrary number survival outcome types, \(H\) (though, identification of parameters quickly becomes an issue).
- Only observing a subset of the \(L\) markers at a given point in time \(s_{ij}\) (i.e. some may be missing).
- Any number of observed markers and survival outcomes for each individual.
- Left-truncation and right-censoring.
- Time-varying effects for covariates (both fixed and random effects) and non-proportional hazard effects for covariates.
The following is not supported
- Interval-censoring and left-censoring are not supported although this is not a complicated extension.
- There cannot be multiple observed markers of the same type at each point \(s_{ij}\) for each individual \(i\). This is not too complicated to implement but not yet implemented.
Examples
We show a few examples of using this package and illustrate its API in this section. The examples includes:
- Only one marker in the Univariate Marker Example section where we can compare with lme4. The section also includes some remarks about the implemented expansions.
- Two markers in the Two Markers section.
- A recurrent event in the Recurrent Event section. It is also shown how to change the quadrature rule that is used to approximate the approximate expected cumulative hazard.
- An example with two markers and a recurrent event is given in the Two Markers and a Recurrent Event section. The section also includes an example of how to construct approximate profile likelihood based confidence intervals and how to obtain Wald type tests and confidence intervals from the observed information matrix.
- A similar example is provided but without the frailty in the Two Markers and a Recurrent Event Without Frailty section. Notice that this is the default with
surv_term. It is also shown how one can get the variational parameters for each individual. - An example with two markers, the observation process, and a terminal event is given in the Two Markers, the Observation Time Process, and a Terminal Event section. The observation time process and the markers are not marginally independent in the true model. Thus, they need to be modeled together or the dependence accounted for in some other way. Examples and comments regarding caching of the expansions are also provided.
- A similar example is provided in the Two Markers, the Observation Time Process, and a Terminal Event with Delayed Entry section but where there is delayed entry.
- A similar example is provided in the Two Markers, the Observation Time Process, and a Terminal Event with Time-varying Effects section but where there time-varying fixed and random covariate effects along with non-proportional hazard effects.
- A similar example is provided in the Two Markers, the Observation Time Process, a Terminal Event and Mixed Dependencies section but where a mixture of the cumulative, the present value, and the derivative is used in the sub-survival models.
We start first with some general remarks.
General Remarks
There is not intercept by default in the splines. Usually, you almost always want one in the both the random effect and fixed effect part. The intercept is included by default in formula (unless you do . ~ ... - 1 and there is not factor variable). Thus, you do not want an intercept in the time_fixef in both marker_term and surv_term. On the other hand, you do need yourself to include the intercept in the time_rng argument of marker_term. Thus, a typical call may look like
mark <- marker_term(
y ~ X1 + X2, # there is an intercept by default
# `intercept` is FALSE by default
time_fixef = poly_term(time, degree = 3),
# we need to add the intercept
time_rng = poly_term(time, degree = 1, intercept = TRUE))
s_term <- surv_term(
Surv(tstop, event) ~ X1, # there is an intercept by default
# `intercept` is FALSE by default
time_fixef = poly_term(tstop, degree = 1))You may run into computational issues if you have a very flexible random effect model. This is quite typical with mixed models in a frequentist setting also with other estimation methods. One outcome of a too flexible random effect model is slow convergence or a NaN output because some of the scale matrices do not have full rank.
Univariate Marker Example
We make a comparison below with one marker only with the lme4 package. We do this as we can check that
- We get almost the same (the lower bound can be equal to the log marginal likelihood in this case).
- The lower bound is less than or equal to the log marginal likelihood.
lme4 should be used in this case as the linear model is tractable and this is used by the lme4 package. The computation time comparison is not fair as we use four threads with the methods in this package.
# settings for the simulation
library(splines)
g_func <- function(x)
ns(x, knots = c(.5, 1.5), Boundary.knots = c(0, 2))
m_func <- function(x)
ns(x, knots = 1, Boundary.knots = c(0, 2), intercept = TRUE)
fixef_vary_marker <- c(1.4, -1.2, -2.1) # beta
fixef_marker <- c(-.5, 1) # gamma
# Psi
vcov_vary <- structure(
c(0.18, 0.05, -0.05, 0.05, 0.34, -0.25, -0.05, -0.25, 0.24),
.Dim = c(3L, 3L))
vcov_marker <- matrix(.6^2, 1) # Sigma
# plot the true mean curve along with 95% confidence pointwise quantiles
library(VAJointSurv)
#> Loading required package: survival
par(mar = c(5, 5, 1, 1))
plot_marker(
time_fixef = ns_term(
knots = c(.5, 1.5), Boundary.knots = c(0, 2)),
time_rng = ns_term(
knots = 1, Boundary.knots = c(0, 2), intercept = TRUE),
fixef_vary = fixef_vary_marker, x_range = c(0, 2), vcov_vary = vcov_vary)
library(mvtnorm)
# simulates from the model by sampling a given number of individuals
sim_dat <- function(n_ids){
# simulate the outcomes
sig_sqrt <- sqrt(vcov_marker)
dat <- lapply(1:n_ids, function(id){
# sample the number of outcomes and the fixed effect covariates
n_obs <- sample.int(8L, 1L)
obs_time <- sort(runif(n_obs, 0, 2))
X <- cbind(1, rnorm(n_obs))
colnames(X) <- paste0("X", 1:NCOL(X) - 1L)
# sample the outcomes
U <- drop(rmvnorm(1, sigma = vcov_vary))
eta <- X %*% fixef_marker + drop(g_func(obs_time)) %*% fixef_vary_marker +
drop(m_func(obs_time)) %*% U
y <- drop(eta) + rnorm(n_obs, sd = sig_sqrt)
cbind(y = y, X = X[, -1, drop = FALSE], time = obs_time, id = id)
})
# combine the data and return
out <- as.data.frame(do.call(rbind, dat))
out$id <- as.integer(out$id)
out[sample.int(NROW(out)), ] # the order does not matter
}
# sample a moderately large data set
set.seed(2)
dat <- sim_dat(2000L)
# the number of observed outcomes
nrow(dat)
#> [1] 8992
# example of one individual's outcomes
subset(dat, id == 1)
#> y X1 time id
#> 5 -3.6206 -0.23970 1.8877 1
#> 2 0.1270 -0.08025 1.1467 1
#> 3 -0.3545 0.13242 1.4047 1
#> 1 -1.1089 -1.13038 0.3361 1
#> 4 -0.2127 0.70795 1.8869 1
# fit the model with lme4
library(lme4)
#> Loading required package: Matrix
system.time(
fit <- lmer(y ~ X1 + g_func(time) + (m_func(time) - 1| id), dat,
control = lmerControl(optimizer = "bobyqa"),
# to compare with the lower bound from this package
REML = FALSE))
#> user system elapsed
#> 0.815 0.008 0.837
# the maximum log likelihood
print(logLik(fit), digits = 8)
#> 'log Lik.' -9161.7723 (df=12)
# estimate the model with this package. Get the object we need for the
# optimization
system.time(comp_obj <- joint_ms_ptr(
markers = marker_term(
y ~ X1, id = id, dat,
time_fixef = ns_term(time, knots = c(.5, 1.5), Boundary.knots = c(0, 2)),
time_rng = ns_term(time, knots = 1, Boundary.knots = c(0, 2),
intercept = TRUE)),
max_threads = 4L))
#> user system elapsed
#> 0.020 0.000 0.021
# get the starting values
system.time(start_val <- joint_ms_start_val(comp_obj))
#> user system elapsed
#> 2.241 0.012 0.948
# lower bound at the starting values
print(-attr(start_val, "value"), digits = 8)
#> [1] -9161.7781
# check that the gradient is correct
f <- function(x){
start_val[seq_along(x)] <- x
joint_ms_lb(comp_obj, start_val)
}
all.equal(numDeriv::grad(f, head(start_val, 12 + 2 * 9)),
head(joint_ms_lb_gr(comp_obj, start_val), 12 + 2 * 9))
#> [1] "Mean relative difference: 1.029e-05"
# find the maximum lower bound estimate
system.time(opt_out <- joint_ms_opt(comp_obj, par = start_val, max_it = 1000L,
cg_tol = .2, c2 = .1))
#> user system elapsed
#> 0.505 0.000 0.128
# check the gradient norm. We may need to reduce the convergence tolerance if
# this is not small. In can also be a sign of convergence issues
sqrt(sum(joint_ms_lb_gr(comp_obj, opt_out$par)^2))
#> [1] 0.7437
opt_out$info # convergence code (0 == 'OK')
#> [1] 0
print(-opt_out$value, digits = 8) # maximum lower bound value
#> [1] -9161.7762
print(logLik(fit), digits = 8) # maximum likelihood from lme4
#> 'log Lik.' -9161.7723 (df=12)
# compare the estimates with the actual values. Start with the fixed effects
fmt_ests <- joint_ms_format(comp_obj, opt_out$par)
rbind(lme4 = fixef(fit),
`This package` = do.call(c, fmt_ests$markers[[1]]),
Actual = c(fixef_marker, fixef_vary_marker))
#> (Intercept) X1 g_func(time)1 g_func(time)2 g_func(time)3
#> lme4 -0.5206 1.007 1.372 -1.104 -2.089
#> This package -0.5207 1.007 1.372 -1.104 -2.089
#> Actual -0.5000 1.000 1.400 -1.200 -2.100
# then the random effect covariance matrix
vcov_est <- VarCorr(fit)[["id"]]
attributes(vcov_est)[c("stddev", "correlation")] <- NULL
vcov_est # lme4
#> m_func(time)1 m_func(time)2 m_func(time)3
#> m_func(time)1 0.16228 0.09634 -0.04534
#> m_func(time)2 0.09634 0.29559 -0.26821
#> m_func(time)3 -0.04534 -0.26821 0.29652
fmt_ests$vcov$vcov_vary # this package
#> [,1] [,2] [,3]
#> [1,] 0.16581 0.0955 -0.04641
#> [2,] 0.09550 0.2955 -0.26790
#> [3,] -0.04641 -0.2679 0.29643
vcov_vary # the actual values
#> [,1] [,2] [,3]
#> [1,] 0.18 0.05 -0.05
#> [2,] 0.05 0.34 -0.25
#> [3,] -0.05 -0.25 0.24
# the error term variance
c(lme4 = attr(VarCorr(fit), "sc")^2,
`This package` = fmt_ests$vcov$vcov_marker,
Actual = vcov_marker)
#> lme4 This package Actual
#> 0.3667 0.3665 0.3600Note on Basis Expansions
This is a technical section which you may skip. Special basis expansions for e.g. poly and ns are implemented in R as poly_term and ns_term. The reason is that the integration shown in the Survival Outcomes section has to be done many times. Thus, we have implemented all the expansions in C++ to be used in C++ to reduce the cost of the evaluations.
The _term functions provide the R interface to the C++ functions. Each return a list with elements like their R counterparts and an eval element. The latter element can be used to evaluate the basis, the derivatives of the basis if the der argument is greater than zero, or the integral
\[\vec d(u) = \int_l^u \vec m(s) ds\]
if the der argument is equal to minus one. The lower limit is passed with the lower_limit argument to the eval function. We illustrate this below first with a raw polynomial. We also show that the derivatives and integral are correct.
# raw polynomial with an without an intercept
pl <- poly_term(degree = 3, raw = TRUE, intercept = TRUE)
pl_no_inter <- poly_term(degree = 3, raw = TRUE, intercept = FALSE)
t(pl $eval(1:2))
#> [,1] [,2] [,3] [,4]
#> [1,] 1 1 1 1
#> [2,] 1 2 4 8
t(pl_no_inter$eval(1:2))
#> [,1] [,2] [,3]
#> [1,] 1 1 1
#> [2,] 2 4 8
outer(1:2, 0:3, `^`) # R version
#> [,1] [,2] [,3] [,4]
#> [1,] 1 1 1 1
#> [2,] 1 2 4 8
# derivatives
library(numDeriv)
t(pl$eval(2:3, der = 1))
#> [,1] [,2] [,3] [,4]
#> [1,] 0 1 4 12
#> [2,] 0 1 6 27
t(pl_no_inter$eval(2:3, der = 1))
#> [,1] [,2] [,3]
#> [1,] 1 4 12
#> [2,] 1 6 27
t(pl$eval(2:3, der = 2)) # second derivative
#> [,1] [,2] [,3] [,4]
#> [1,] 0 0 2 12
#> [2,] 0 0 2 18
t(pl_no_inter$eval(2:3, der = 2)) # second derivative
#> [,1] [,2] [,3]
#> [1,] 0 2 12
#> [2,] 0 2 18
# trivial to do compute by hand but we do it with numDeriv anyway
t(sapply(2:3, function(x) jacobian(function(z) outer(z, 0:3, `^`), x)))
#> [,1] [,2] [,3] [,4]
#> [1,] 0 1 4 12
#> [2,] 0 1 6 27
# integral
t(pl$eval(3:4, der = -1, lower_limit = 0))
#> [,1] [,2] [,3] [,4]
#> [1,] 3 4.5 9.00 20.25
#> [2,] 4 8.0 21.33 64.00
t(pl_no_inter$eval(3:4, der = -1, lower_limit = 0))
#> [,1] [,2] [,3]
#> [1,] 4.5 9.00 20.25
#> [2,] 8.0 21.33 64.00
# we could do this by hand but we write a function which we can reuse
do_int <- function(xs, f, lower_limit){
n_terms <- length(f(xs[1]))
g <- function(x, i) f(x)[, i]
outer(xs, 1:n_terms, Vectorize(
function(x, i) integrate(g, lower_limit, x, i = i)$value))
}
do_int(3:4, function(x) outer(x, 0:3, `^`), 0)
#> [,1] [,2] [,3] [,4]
#> [1,] 3 4.5 9.00 20.25
#> [2,] 4 8.0 21.33 64.00The functionality also works for orthogonal polynomials.
# equally works with orthogonal polynomials
xx <- seq(0, 4, length.out = 10)
pl <- poly_term(xx, degree = 3, raw = FALSE, intercept = TRUE)
t(pl$eval(2:3))
#> [,1] [,2] [,3] [,4]
#> [1,] 1 5.500e-17 -0.3590 2.391e-16
#> [2,] 1 2.477e-01 -0.1387 -3.881e-01
predict(poly(xx, degree = 3), 2:3)
#> 1 2 3
#> [1,] 5.500e-17 -0.3590 2.391e-16
#> [2,] 2.477e-01 -0.1387 -3.881e-01
# derivatives
t(pl$eval(2:3, der = 1))
#> [,1] [,2] [,3] [,4]
#> [1,] 0 0.2477 1.110e-16 -0.59310
#> [2,] 0 0.2477 4.406e-01 0.02176
t(sapply(2:3, function(x)
jacobian(function(z) predict(poly(xx, degree = 3), z), x)))
#> [,1] [,2] [,3]
#> [1,] 0.2477 0.0000 -0.59310
#> [2,] 0.2477 0.4406 0.02176
# integral
t(pl$eval(3:4, der = -1, lower_limit = 0))
#> [,1] [,2] [,3] [,4]
#> [1,] 3 -3.716e-01 -0.4162 1.211e-01
#> [2,] 4 2.220e-16 -0.2611 1.776e-15
do_int(3:4, function(x) predict(poly(xx, degree = 3), x), 0)
#> [,1] [,2] [,3]
#> [1,] -3.716e-01 -0.4162 1.211e-01
#> [2,] 2.248e-16 -0.2611 5.231e-16The derivatives and integral are also available for B-splines and natural cubic splines through the ns_term and bs_term functions.
# B-splines
library(splines)
xx <- 0:10
b <- bs_term(xx, df = 4, Boundary.knots = range(xx))
t(b$eval(4.33))
#> [,1] [,2] [,3] [,4]
#> [1,] 0.3598 0.4755 0.1624 0
c(predict(bs(xx, df = 4, Boundary.knots = range(xx)), 4.33))
#> [1] 0.3598 0.4755 0.1624 0.0000
# derivatives
t(b$eval(4.33, der = 1))
#> [,1] [,2] [,3] [,4]
#> [1,] -0.1713 0.06963 0.1125 0
f <- function(z) suppressWarnings(
predict(bs(xx, df = 4, Boundary.knots = range(xx)), z))
t(jacobian(f, 4.33))
#> [,1] [,2] [,3] [,4]
#> [1,] -0.1713 0.06963 0.1125 0
# integrals
all.equal( t(b$eval(12, der = -1, lower_limit = 11)), do_int(12, f, 11), 1e-6)
#> [1] TRUE
all.equal( t(b$eval(11, der = -1, lower_limit = 0)), do_int(11, f, 0), 1e-6)
#> [1] TRUE
all.equal(-t(b$eval( 0, der = -1, lower_limit = 11)), do_int(11, f, 0), 1e-6)
#> [1] TRUE
all.equal( t(b$eval( 5, der = -1, lower_limit = 1)), do_int( 5, f, 1), 1e-6)
#> [1] TRUE
all.equal( t(b$eval(-1, der = -1, lower_limit = 2)), do_int(-1, f, 2), 1e-6)
#> [1] TRUE
# natural cubic spline
xx <- 0:10
n <- ns_term(xx, df = 3, Boundary.knots = range(xx))
t(n$eval(4.33))
#> [,1] [,2] [,3]
#> [1,] 0.1723 0.5256 -0.346
c(predict(ns(xx, df = 3, Boundary.knots = range(xx)), 4.33))
#> [1] 0.1723 0.5256 -0.3460
# derivatives
t(n$eval(4.33, der = 1))
#> [,1] [,2] [,3]
#> [1,] 0.2186 -0.05737 0.05165
f <- function(z) predict(ns(xx, df = 3, Boundary.knots = range(xx)), z)
t(jacobian(f, 4.33))
#> [,1] [,2] [,3]
#> [1,] 0.2186 -0.05737 0.05165
# integrals
all.equal( t(n$eval(12, der = -1, lower_limit = 11)), do_int(12, f, 11), 1e-6)
#> [1] TRUE
all.equal( t(n$eval(11, der = -1, lower_limit = 0)), do_int(11, f, 0), 1e-6)
#> [1] TRUE
all.equal(-t(n$eval( 0, der = -1, lower_limit = 11)), do_int(11, f, 0), 1e-6)
#> [1] TRUE
all.equal( t(n$eval( 7, der = -1, lower_limit = -2)), do_int( 7, f, -2), 1e-6)
#> [1] TRUE
all.equal( t(n$eval(-1, der = -1, lower_limit = -3)), do_int(-1, f, -3), 1e-6)
#> [1] TRUECaution: the B-splines are only tested with degree = 3!
Log Transformations
The expansions can also be used on the log scale. This is mainly implemented to allow the log of the baseline hazard to be parameterized in terms of an expansions in log time to have the Weibull model as a special case. A few examples are given below.
# raw log polynomial with an without an intercept
pl <- poly_term(degree = 3, raw = TRUE, intercept = TRUE, use_log = TRUE)
pl_no_inter <- poly_term(degree = 3, raw = TRUE, intercept = FALSE,
use_log = TRUE)
t(pl $eval(1:2))
#> [,1] [,2] [,3] [,4]
#> [1,] 1 0.0000 0.0000 0.000
#> [2,] 1 0.6931 0.4805 0.333
t(pl_no_inter$eval(1:2))
#> [,1] [,2] [,3]
#> [1,] 0.0000 0.0000 0.000
#> [2,] 0.6931 0.4805 0.333
outer(log(1:2), 0:3, `^`) # R version
#> [,1] [,2] [,3] [,4]
#> [1,] 1 0.0000 0.0000 0.000
#> [2,] 1 0.6931 0.4805 0.333
# derivatives
library(numDeriv)
t(pl$eval(2:3, der = 1))
#> [,1] [,2] [,3] [,4]
#> [1,] 0 0.5000 0.6931 0.7207
#> [2,] 0 0.3333 0.7324 1.2069
t(pl_no_inter$eval(2:3, der = 1))
#> [,1] [,2] [,3]
#> [1,] 0.5000 0.6931 0.7207
#> [2,] 0.3333 0.7324 1.2069
# trivial to do compute by hand but we do it with numDeriv anyway
t(sapply(2:3, function(x) jacobian(function(z) outer(log(z), 0:3, `^`), x)))
#> [,1] [,2] [,3] [,4]
#> [1,] 0 0.5000 0.6931 0.7207
#> [2,] 0 0.3333 0.7324 1.2069# B-splines
library(splines)
xx <- 1:10
b <- bs_term(xx, df = 4, use_log = TRUE)
t(b$eval(4.33))
#> [,1] [,2] [,3] [,4]
#> [1,] 0.1736 0.4746 0.3491 0
c(predict(bs(log(xx), df = 4), log(4.33)))
#> [1] 0.1736 0.4746 0.3491 0.0000
# derivatives
t(b$eval(4.33, der = 1))
#> [,1] [,2] [,3] [,4]
#> [1,] -0.1223 -0.03495 0.165 0
f <- function(z) suppressWarnings(predict(bs(log(xx), df = 4), log(z)))
t(jacobian(f, 4.33))
#> [,1] [,2] [,3] [,4]
#> [1,] -0.1223 -0.03495 0.165 0
# natural cubic spline
xx <- 1:10
n <- ns_term(xx, df = 3, use_log = TRUE)
t(n$eval(4.33))
#> [,1] [,2] [,3]
#> [1,] 0.3855 0.4276 -0.307
c(predict(ns(log(xx), df = 3), log(4.33)))
#> [1] 0.3855 0.4276 -0.3070
# derivatives
t(n$eval(4.33, der = 1))
#> [,1] [,2] [,3]
#> [1,] 0.2536 -0.09559 0.07548
f <- function(z) predict(ns(log(xx), df = 3), log(z))
t(jacobian(f, 4.33))
#> [,1] [,2] [,3]
#> [1,] 0.2536 -0.09559 0.07548Time-varying Effects
We may also be interested in
- Time-varying fixed effects for a covariate.
- Time-varying random effects for a covariate.
- Non-proportional hazard effects for a covariate.
The way that we allow for this is with the weighted_term and stacked_term functions. The weighted_term function takes a basis expansion which we can denote as \(\vec g(t)\) and a symbol (weight) for a covariate which we denote by \(x\) and computes \(x\vec g(t)\). Two examples of using the function are given below. Note that the symbol is not evaluated till you call the eval function on the returned object.
# define two bases
term_bmi <- weighted_term(
poly_term(degree = 2, intercept = TRUE, raw = TRUE),
weight = bmi)
term_sex <- weighted_term(
poly_term(degree = 1, intercept = TRUE, raw = TRUE),
weight = is_male)
# use the two bases in an example
ti <- 2:3
df <- data.frame(bmi = c(23, 25), is_male = c(0, 1))
term_bmi$eval(x = ti, newdata = df)
#> [,1] [,2]
#> [1,] 23 25
#> [2,] 46 75
#> [3,] 92 225
term_sex$eval(x = ti, newdata = df)
#> [,1] [,2]
#> [1,] 0 1
#> [2,] 0 3The stacked_term function takes a number of expansions and stack them on top of each other. This illustrated below.
# add an "intercept" to the two other terms combined
term_comb <- stacked_term(poly_term(degree = 2, intercept = TRUE, raw = TRUE),
term_bmi, term_sex)
term_comb$eval(ti, newdata = df)
#> [,1] [,2]
#> [1,] 1 1
#> [2,] 2 3
#> [3,] 4 9
#> [4,] 23 25
#> [5,] 46 75
#> [6,] 92 225
#> [7,] 0 1
#> [8,] 0 3
# it is a bit obscure but we can also weight the combined term
term_comb_sex <- weighted_term(term_comb, is_male)
term_comb_sex$eval(ti, newdata = df)
#> [,1] [,2]
#> [1,] 0 1
#> [2,] 0 3
#> [3,] 0 9
#> [4,] 0 25
#> [5,] 0 75
#> [6,] 0 225
#> [7,] 0 1
#> [8,] 0 3Two Markers
We provide an example with two markers instead of one. We observe one of the two markers or both at each observation time. There is also a closed form solution in this case. Thus, this is not the optimal way of doing this and it is only shown as an example and as a sanity check.
# settings for the simulation
library(splines)
g_funcs <- list(
function(x)
ns(x, knots = c(.5, 1.5), Boundary.knots = c(0, 2)),
function(x)
# a raw polynomial
outer(x, 1:2, `^`))
m_funcs <- list(
function(x)
ns(x, knots = numeric(), Boundary.knots = c(0, 2), intercept = TRUE),
function(x)
# a raw polynomial
outer(x, 0:1, `^`))
fixef_vary_marker <- list(c(1.4, -1.2, -2.1), c(.5, .67)) # beta
fixef_marker <- list(c(-.5, 1), .25) # gamma
# Psi
vcov_vary <- structure(
c(0.35, 0.02, -0.05, 0.01, 0.02, 0.12, -0.06, -0.01, -0.05, -0.06, 0.32, 0.09, 0.01, -0.01, 0.09, 0.12),
.Dim = c(4L, 4L))
vcov_marker <- matrix(c(.6^2, .1, .1, .4^2), 2) # Sigma
# plot the markers' mean curve
library(VAJointSurv)
par(mar = c(5, 5, 1, 1))
plot_marker(
time_fixef = ns_term(
knots = c(.5, 1.5), Boundary.knots = c(0, 2)),
time_rng = ns_term(
knots = numeric(), Boundary.knots = c(0, 2), intercept = TRUE),
fixef_vary = fixef_vary_marker[[1]], x_range = c(0, 2),
vcov_vary = vcov_vary[1:2, 1:2], ylab = "Marker 1")
plot_marker(
time_fixef = poly_term(degree = 2, raw = TRUE),
time_rng = poly_term(degree = 1, raw = TRUE, intercept = TRUE),
fixef_vary = fixef_vary_marker[[2]], x_range = c(0, 2),
vcov_vary = vcov_vary[3:4, 3:4], ylab = "Marker 2")
library(mvtnorm)
# simulates from the model by sampling a given number of individuals
sim_dat <- function(n_ids){
# simulate the outcomes
dat <- lapply(1:n_ids, function(id){
# sample the number of outcomes and the fixed effect covariates
n_obs <- sample.int(8L, 1L)
obs_time <- sort(runif(n_obs, 0, 2))
X1 <- cbind(1, rnorm(n_obs))
X2 <- matrix(1, n_obs)
colnames(X1) <- paste0("X1_", 1:NCOL(X1) - 1L)
colnames(X2) <- paste0("X2_", 1:NCOL(X2) - 1L)
X <- list(X1, X2)
# sample the outcomes
U <- drop(rmvnorm(1, sigma = vcov_vary))
eta <- sapply(1:2, function(i)
X[[i]] %*% fixef_marker[[i]] +
drop(g_funcs[[i]](obs_time)) %*% fixef_vary_marker[[i]] +
drop(m_funcs[[i]](obs_time)) %*% U[1:2 + (i == 2) * 2])
ys <- eta + rmvnorm(n_obs, sigma = vcov_marker)
colnames(ys) <- paste0("Y", 1:2)
# mask some observations
do_mask <- sample.int(3L, n_obs, replace = TRUE)
ys[do_mask == 2, 1] <- NA
ys[do_mask == 3, 2] <- NA
X <- do.call(cbind, lapply(X, function(x) x[, -1, drop = FALSE]))
cbind(ys, X, time = obs_time, id = id)
})
# combine the data and return
out <- as.data.frame(do.call(rbind, dat))
out$id <- as.integer(out$id)
out[sample.int(NROW(out)), ] # the order does not matter
}
# sample a moderately large data set
set.seed(2)
dat <- sim_dat(1000L)
# example of the data for one individual
subset(dat, id == 1)
#> Y1 Y2 X1_1 time id
#> 3 NA 3.1915 0.13242 1.4047 1
#> 5 NA 4.3373 -0.23970 1.8877 1
#> 4 NA 4.8558 0.70795 1.8869 1
#> 2 1.399 1.7030 -0.08025 1.1467 1
#> 1 -1.286 0.4379 -1.13038 0.3361 1
# estimate the model with this package. Get the object we need for the
# optimization
marker_1 <- marker_term(
Y1 ~ X1_1, id = id, subset(dat, !is.na(Y1)),
time_fixef = ns_term(time, knots = c(.5, 1.5), Boundary.knots = c(0, 2)),
time_rng = ns_term(time, knots = numeric(), Boundary.knots = c(0, 2),
intercept = TRUE))
marker_2 <- marker_term(
Y2 ~ 1, id = id, subset(dat, !is.na(Y2)),
time_fixef = poly_term(time, degree = 2, raw = TRUE),
time_rng = poly_term(time, degree = 1, raw = TRUE, intercept = TRUE))
comp_obj <- joint_ms_ptr(markers = list(marker_1, marker_2),
max_threads = 4L)
rm(marker_1, marker_2)
# get the starting values
system.time(start_val <- joint_ms_start_val(comp_obj))
#> user system elapsed
#> 2.086 0.004 0.615
# lower bound at the starting values
print(-attr(start_val, "value"), digits = 8)
#> [1] -5729.7904
# check that the gradient is correct
f <- function(x){
start_val[seq_along(x)] <- x
joint_ms_lb(comp_obj, start_val)
}
all.equal(numDeriv::grad(f, head(start_val, 22 + 2 * 14)),
head(joint_ms_lb_gr(comp_obj, start_val), 22 + 2 * 14))
#> [1] "Mean relative difference: 0.0001464"
# find the maximum lower bound estimate
system.time(opt_out <- joint_ms_opt(comp_obj, par = start_val, max_it = 1000L,
pre_method = 3L, cg_tol = .2, c2 = .1,
gr_tol = .1))
#> user system elapsed
#> 0.342 0.000 0.086
# we set gr_tol in the call so this is the convergence criterion for the
# gradient
sqrt(sum(joint_ms_lb_gr(comp_obj, opt_out$par)^2))
#> [1] 0.09769
opt_out$info # convergence code (0 == 'OK')
#> [1] 0
print(-opt_out$value, digits = 8) # maximum lower bound value
#> [1] -5729.7904
opt_out$counts
#> function gradient n_cg
#> 42 30 41
# find the maximum lower bound with lbfgs
library(lbfgsb3c)
system.time(lbfgs_res <- lbfgsb3c(
start_val, function(x) joint_ms_lb(comp_obj, x),
function(x) joint_ms_lb_gr(comp_obj, x),
control = list(factr = 1e-8, maxit = 10000L)))
#> user system elapsed
#> 36.863 0.016 9.284
lbfgs_res$convergence # convergence code (0 == 'OK')
#> [1] 0
print(-lbfgs_res$value, digits = 8) # maximum lower bound value
#> [1] -5729.7903
lbfgs_res$counts # may have hit maxit!
#> [1] 4108 4108
# compare the estimates with the actual values. Start with the fixed effects
fmt_ests <- joint_ms_format(comp_obj, opt_out$par)
fmt_ests_lbfgs <- joint_ms_format(comp_obj, lbfgs_res$par)
# the parameters for the first marker
fmt_ests$markers[[1]]
#> $fixef
#> [1] -0.4967 1.0249
#>
#> $fixef_vary
#> [1] 1.544 -1.237 -2.114
fmt_ests_lbfgs$markers[[1]] # with lbfgs
#> $fixef
#> [1] -0.4966 1.0249
#>
#> $fixef_vary
#> [1] 1.544 -1.238 -2.114
fixef_marker[[1]] # true values
#> [1] -0.5 1.0
fixef_vary_marker[[1]] # true values
#> [1] 1.4 -1.2 -2.1
# the parameters for the second marker
fmt_ests$markers[[2]]
#> $fixef
#> [1] 0.2552
#>
#> $fixef_vary
#> [1] 0.5381 0.6592
fmt_ests_lbfgs$markers[[2]] # with lbfgs
#> $fixef
#> [1] 0.2552
#>
#> $fixef_vary
#> [1] 0.5382 0.6591
fixef_marker[[2]] # true values
#> [1] 0.25
fixef_vary_marker[[2]] # true values
#> [1] 0.50 0.67
# the parameters for covariance matrix of the random effects
fmt_ests$vcov$vcov_vary
#> [,1] [,2] [,3] [,4]
#> [1,] 0.4155040 0.02347 -0.006174 0.0004588
#> [2,] 0.0234666 0.19181 -0.087433 0.0167084
#> [3,] -0.0061743 -0.08743 0.331594 0.0688998
#> [4,] 0.0004588 0.01671 0.068900 0.0948096
fmt_ests_lbfgs$vcov$vcov_vary # with lbfgs
#> [,1] [,2] [,3] [,4]
#> [1,] 0.4155511 0.02345 -0.006134 0.0004161
#> [2,] 0.0234507 0.19154 -0.087515 0.0167512
#> [3,] -0.0061338 -0.08752 0.331641 0.0688570
#> [4,] 0.0004161 0.01675 0.068857 0.0948671
vcov_vary # the true values
#> [,1] [,2] [,3] [,4]
#> [1,] 0.35 0.02 -0.05 0.01
#> [2,] 0.02 0.12 -0.06 -0.01
#> [3,] -0.05 -0.06 0.32 0.09
#> [4,] 0.01 -0.01 0.09 0.12
# the parameters for the error term variance
fmt_ests$vcov$vcov_marker
#> [,1] [,2]
#> [1,] 0.3575 0.1003
#> [2,] 0.1003 0.1601
fmt_ests_lbfgs$vcov$vcov_marker # with lbfgs
#> [,1] [,2]
#> [1,] 0.3575 0.1003
#> [2,] 0.1003 0.1601
vcov_marker
#> [,1] [,2]
#> [1,] 0.36 0.10
#> [2,] 0.10 0.16Recurrent Event
We simulate from a model with a recurrent event in this section. This is a more relevant example as there is not a closed form solution of the marginal likelihood.
# the survival parameters
fixef_surv <- c(-.5, .25) # delta
fixef_vary_surv <- c(.5, .1, -.2, .11) # omega
vcov_surv <- matrix(.2^2, 1) # Xi
library(splines)
b_func <- function(x)
bs(x, knots = 1, Boundary.knots = c(0, 2))
# sample a few survival curves and plot them
library(SimSurvNMarker)
# time points where we evaluate the conditional survival functions
tis <- seq(0, 2, length.out = 50)
set.seed(1)
# compute the conditional survival functions disregarding the fixed effect
surv_vals <- replicate(100, {
err <- rnorm(1, sd = sqrt(vcov_surv))
eval_surv_base_fun(
ti = tis, omega = fixef_vary_surv, b_func = b_func,
gl_dat = get_gl_rule(100), delta = fixef_surv[1] + err)
})
par(mar = c(5, 5, 1, 1))
matplot(tis, surv_vals, type = "l", col = gray(0, .1), bty = "l", lty = 1,
ylim = c(0, 1), xaxs = "i", yaxs = "i", xlab = "Time",
ylab = "Conditional survival function")
# compute the average to get an estimate of the marginal survival curve
plot(tis, rowMeans(surv_vals), type = "l", bty = "l", ylim = c(0, 1),
xaxs = "i", yaxs = "i", xlab = "Time",
ylab = "Marginal survival function")
# simulates from the model by sampling a given number of individuals
sim_dat <- function(n_ids){
# simulate the outcomes
gl_dat <- get_gl_rule(100L)
dat <- lapply(1:n_ids, function(id){
# simulate the survival outcome
rng_surv <- rnorm(1, sd = sqrt(vcov_surv))
Z <- c(1, runif(1, -1, 1))
log_haz_offset <- sum(Z * fixef_surv) + rng_surv
# the conditional survival function
surv_func <- function(ti)
eval_surv_base_fun(
ti = ti, omega = fixef_vary_surv, b_func = b_func, gl_dat = gl_dat,
delta = log_haz_offset)
# simulate the recurrent events one at a time
max_sample <- 10L
Z <- matrix(rep(Z, each = max_sample), max_sample)
event <- y <- lf_trunc <- rep(NA_real_, max_sample)
lf_trunc_i <- 0
rngs <- runif(max_sample)
left_trunc_surv <- 1
for(i in 1:max_sample){
# sample a random uniform variable and invert the survival function
rng_i <- rngs[i]
lf_trunc[i] <- lf_trunc_i
root_func <- function(x) rng_i - surv_func(x) / left_trunc_surv
if(root_func(2) < 0){
# the observation is right-censored and we can exit
y[i] <- 2
event[i] <- 0
break
}
# we need to invert the survival function to find the observation time
root <- uniroot(root_func, c(lf_trunc_i, 2), tol = 1e-6)
lf_trunc_i <- y[i] <- root$root
event[i] <- 1
left_trunc_surv <- surv_func(y[i])
}
# keep the needed data and return
colnames(Z) <- paste0("Z_", 1:NCOL(Z) - 1L)
surv_data <- cbind(lf_trunc = lf_trunc[1:i], y = y[1:i], event = event[1:i],
Z[1:i, -1, drop = FALSE], id = id)
list(surv_data = surv_data)
})
# combine the data and return
surv_data <- as.data.frame(do.call(
rbind, lapply(dat, `[[`, "surv_data")))
surv_data$id <- as.integer(surv_data$id)
# the order does not matter
surv_data <- surv_data[sample.int(NROW(surv_data)), ]
list(surv_data = surv_data)
}
# sample a moderately large data set
set.seed(3)
dat <- sim_dat(1000L)
# we show a few properties of the data below
NROW(dat$surv_data) # number of survival outcomes
#> [1] 2388
sum(dat$surv_data$event) # number of observed events
#> [1] 1388
# example of the data for three individuals
subset(dat$surv_data, id %in% 1:3)
#> lf_trunc y event Z_1 id
#> 5 1.6383 2.0000 0 0.1208 3
#> 3 0.0000 0.3526 1 0.1208 3
#> 4 0.3526 1.6383 1 0.1208 3
#> 1 0.0000 2.0000 0 -0.2301 1
#> 2 0.0000 2.0000 0 0.6594 2
# distribution of observed events per individual
proportions(table(table(dat$surv_data$id) - 1L))
#>
#> 0 1 2 3 4 5 6
#> 0.266 0.343 0.220 0.103 0.050 0.012 0.006
# estimate the model with this package. Get the object we need for the
# optimization
library(survival)
library(VAJointSurv)
surv_obj <- surv_term(
Surv(lf_trunc, y, event) ~ Z_1, id = id, dat$surv_data,
time_fixef = bs_term(y, knots = 1, Boundary.knots = c(0, 2)),
with_frailty = TRUE)
comp_obj <- joint_ms_ptr(markers = list(),
survival_terms = surv_obj, max_threads = 4L)
rm(surv_obj)
# get the starting values
system.time(start_val <- joint_ms_start_val(comp_obj))
#> user system elapsed
#> 0.244 0.000 0.160
# lower bound at the starting values
print(-attr(start_val, "value"), digits = 8)
#> [1] -1864.1515
# check that the gradient is correct
f <- function(x){
start_val <- comp_obj$start_val
start_val[seq_along(x)] <- x
joint_ms_lb(comp_obj, start_val)
}
all.equal(numDeriv::grad(f, head(comp_obj$start_val, 7 + 2 * 2)),
head(joint_ms_lb_gr(comp_obj, comp_obj$start_val), 7 + 2 * 2),
tolerance = 1e-6)
#> [1] TRUE
# find the maximum lower bound estimate
system.time(opt_out <- joint_ms_opt(comp_obj, par = start_val, max_it = 1000L,
pre_method = 3L, cg_tol = .2, c2 = .1))
#> user system elapsed
#> 0.813 0.000 0.206
# check the gradient norm. We may need to reduce the convergence tolerance if
# this is not small. In can also be a sign of convergence issues
sqrt(sum(joint_ms_lb_gr(comp_obj, opt_out$par)^2))
#> [1] 0.1184
opt_out$info # convergence code (0 == 'OK')
#> [1] 0
print(-opt_out$value, digits = 8) # maximum lower bound value
#> [1] -1863.2795
opt_out$counts
#> function gradient n_cg
#> 156 101 165
# check the estimates
fmt_est <- joint_ms_format(comp_obj, opt_out$par)
rbind(estimate = fmt_est$survival[[1]]$fixef, truth = fixef_surv)
#>
#> estimate -0.5058 0.2582
#> truth -0.5000 0.2500
rbind(estimate = fmt_est$survival[[1]]$fixef_vary, truth = fixef_vary_surv)
#>
#> estimate 0.3965 0.2539 -0.3217 0.07835
#> truth 0.5000 0.1000 -0.2000 0.11000
c(estimate = sqrt(fmt_est$vcov$vcov_surv), truth = sqrt(vcov_surv))
#> estimate truth
#> 0.228 0.200Note on Quadrature Rule
A quadrature rule is used to integrate the expected cumulative hazard. The user can pass any quadrature rule with \(n\) nodes and weights each denoted by \((n_i, w_i)\) such that
\[\int_0^1f(x) dx \approx \sum_{i = 1}^nw_i f(n_i).\]
By default Gauss–Legendre quadrature is used. A very simple alternative is to use the midpoint rule. This is illustrated below.
# computes the midpoint rule
mid_rule <- function(n)
list(node = seq(0, 1, length.out = n), weight = rep(1 / n, n))
# compare the midpoint with different number of nodes
for(n in 2^(3:12)){
res <- joint_ms_lb(comp_obj, opt_out$par, n_threads = 4L,
quad_rule = mid_rule(n))
cat(sprintf("# nodes, lower bound %5d %12.4f\n", n, res))
}
#> # nodes, lower bound 8 1859.8265
#> # nodes, lower bound 16 1861.5170
#> # nodes, lower bound 32 1862.3858
#> # nodes, lower bound 64 1862.8294
#> # nodes, lower bound 128 1863.0536
#> # nodes, lower bound 256 1863.1663
#> # nodes, lower bound 512 1863.2228
#> # nodes, lower bound 1024 1863.2511
#> # nodes, lower bound 2048 1863.2652
#> # nodes, lower bound 4096 1863.2723
# the maximum lower bound value we got with Gauss–Legendre quadrature
print(-opt_out$value, digits = 8)
#> [1] -1863.2795We can compare this with using Gauss–Legendre quadrature with different number of nodes.
# get the rules
library(gaussquad)
#> Loading required package: polynom
#> Loading required package: orthopolynom
rules <- legendre.quadrature.rules(2^8, normalized = TRUE)
# rescale the rule
rules <- lapply(rules, function(x)
list(node = .5 * (x[, "x"] + 1), weight = .5 * x[, "w"]))
# use different number of nodes
for(n in 2^(2:8)){
res <- joint_ms_lb(comp_obj, opt_out$par, n_threads = 4L,
quad_rule = rules[[n]])
cat(sprintf("# nodes, lower bound %5d %12.8f\n", n, res))
}
#> # nodes, lower bound 4 1863.11441973
#> # nodes, lower bound 8 1863.26694567
#> # nodes, lower bound 16 1863.27861375
#> # nodes, lower bound 32 1863.27932308
#> # nodes, lower bound 64 1863.27936973
#> # nodes, lower bound 128 1863.27937270
#> # nodes, lower bound 256 1863.27937288
# the result we got
print(-opt_out$value, digits = 14)
#> [1] -1863.2795232813
length(comp_obj$quad_rule$node) # the number of nodes we used
#> [1] 25Evaluating the approximate expected cumulative hazard is by far the most computationally expensive part of the likelihood. Therefore, it may be nice to use fewer nodes e.g. when finding starting values.
Two Markers and a Recurrent Event
We simulate from a model with two different types of markers and a recurrent event in this section.
A Weibull model is selected for the baseline hazard by taking a polynomial of degree one and using the log of time. As an illustration, we select a more flexible expansion in the baseline hazard in the model we estimate.
# settings for the simulation
library(splines)
g_funcs <- list(
function(x)
ns(x, knots = c(.5, 1.5), Boundary.knots = c(0, 2)),
function(x)
# a raw polynomial
outer(x, 1:2, `^`))
m_funcs <- list(
function(x)
ns(x, knots = numeric(), Boundary.knots = c(0, 2), intercept = TRUE),
function(x)
# a raw polynomial
outer(x, 0:1, `^`))
fixef_vary_marker <- list(c(1.4, -1.2, -2.1), c(.5, .67)) # beta
fixef_marker <- list(c(-.5, 1), .25) # gamma
# Psi
vcov_vary <- structure(
c(0.35, 0.02, -0.05, 0.01, 0.02, 0.12, -0.06, -0.01, -0.05, -0.06, 0.32, 0.09, 0.01, -0.01, 0.09, 0.12),
.Dim = c(4L, 4L))
vcov_marker <- matrix(c(.6^2, .1, .1, .4^2), 2) # Sigma
# plot the markers' mean curve
library(VAJointSurv)
par(mar = c(5, 5, 1, 1))
plot_marker(
time_fixef = ns_term(
knots = c(.5, 1.5), Boundary.knots = c(0, 2)),
time_rng = ns_term(
knots = numeric(), Boundary.knots = c(0, 2), intercept = TRUE),
fixef_vary = fixef_vary_marker[[1]], x_range = c(0, 2),
vcov_vary = vcov_vary[1:2, 1:2], ylab = "Marker 1")
plot_marker(
time_fixef = poly_term(degree = 2, raw = TRUE),
time_rng = poly_term(degree = 1, raw = TRUE, intercept = TRUE),
fixef_vary = fixef_vary_marker[[2]], x_range = c(0, 2),
vcov_vary = vcov_vary[3:4, 3:4], ylab = "Marker 2")
# the survival parameters
fixef_surv <- c(-.5, .25) # delta
associations <- c(-.8, .7) # alpha
fixef_vary_surv <- c(.5) # omega
vcov_surv <- matrix(.2^2, 1) # Xi
b_term <- poly_term(degree = 1L, use_log = TRUE, raw = TRUE)
b_func <- function(x)
t(b_term$eval(x))
# plot the hazard with pointwise quantiles
plot_surv(
time_fixef = b_term,
time_rng = list(
ns_term(knots = numeric(), Boundary.knots = c(0, 2), intercept = TRUE),
poly_term(degree = 1, raw = TRUE, intercept = TRUE)),
x_range = c(0, 2), vcov_vary = vcov_vary, frailty_var = vcov_surv,
ps = c(.25, .5, .75), log_hazard_shift = fixef_surv[1],
fixef_vary = fixef_vary_surv, associations = associations)
# without the markers
plot_surv(
time_fixef = b_term,
time_rng = list(
ns_term(knots = numeric(), Boundary.knots = c(0, 2), intercept = TRUE),
poly_term(degree = 1, raw = TRUE, intercept = TRUE)),
x_range = c(0, 2), vcov_vary = matrix(0, 4, 4), frailty_var = vcov_surv,
ps = c(.25, .5, .75), log_hazard_shift = fixef_surv[1],
fixef_vary = fixef_vary_surv, associations = associations)
# compute a few survival curves and plot them
library(SimSurvNMarker)
# time points where we evaluate the conditional survival functions
tis <- seq(0, 2, length.out = 50)
set.seed(1)
Us <- mvtnorm::rmvnorm(250, sigma = vcov_vary) # the random effects
# compute the conditional survival functions disregarding the fixed effect
surv_vals <- apply(Us, 1, function(U){
expansion <- function(x)
cbind(b_func(x), m_funcs[[1]](x) %*% U[1:2],
m_funcs[[2]](x) %*% U[3:4])
err <- rnorm(1, sd = sqrt(vcov_surv))
eval_surv_base_fun(
ti = tis, omega = c(fixef_vary_surv, associations), b_func = expansion,
gl_dat = get_gl_rule(100), delta = fixef_surv[1] + err)
})
matplot(tis, surv_vals, type = "l", col = gray(0, .1), bty = "l", lty = 1,
ylim = c(0, 1), xaxs = "i", yaxs = "i", xlab = "Time",
ylab = "Conditional survival function")
# compute the average to get an estimate of the marginal survival curve
plot(tis, rowMeans(surv_vals), type = "l", bty = "l", ylim = c(0, 1),
xaxs = "i", yaxs = "i", xlab = "Time",
ylab = "Marginal survival function")
library(mvtnorm)
# simulates from the model by sampling a given number of individuals
sim_dat <- function(n_ids){
# simulate the outcomes
gl_dat <- get_gl_rule(100L)
dat <- lapply(1:n_ids, function(id){
# sample the number of outcomes and the fixed effect covariates
n_obs <- sample.int(8L, 1L)
obs_time <- sort(runif(n_obs, 0, 2))
X1 <- cbind(1, rnorm(n_obs))
X2 <- matrix(1, n_obs)
colnames(X1) <- paste0("X1_", 1:NCOL(X1) - 1L)
colnames(X2) <- paste0("X2_", 1:NCOL(X2) - 1L)
X <- list(X1, X2)
# sample the outcomes
U <- drop(rmvnorm(1, sigma = vcov_vary))
eta <- sapply(1:2, function(i)
X[[i]] %*% fixef_marker[[i]] +
drop(g_funcs[[i]](obs_time)) %*% fixef_vary_marker[[i]] +
drop(m_funcs[[i]](obs_time)) %*% U[1:2 + (i == 2) * 2])
ys <- eta + rmvnorm(n_obs, sigma = vcov_marker)
colnames(ys) <- paste0("Y", 1:2)
# mask some observations
do_mask <- sample.int(3L, n_obs, replace = TRUE)
ys[do_mask == 2, 1] <- NA
ys[do_mask == 3, 2] <- NA
X <- do.call(cbind, lapply(X, function(x) x[, -1, drop = FALSE]))
marker_data <- cbind(ys, X, time = obs_time, id = id)
# clean up
rm(list = setdiff(ls(), c("U", "marker_data", "id")))
# simulate the survival outcome
rng_surv <- rnorm(1, sd = sqrt(vcov_surv))
Z <- c(1, runif(1, -1, 1))
log_haz_offset <- sum(Z * fixef_surv) + rng_surv
expansion <- function(x)
cbind(b_func(x), m_funcs[[1]](x) %*% U[1:2],
m_funcs[[2]](x) %*% U[3:4])
# the conditional survival function
surv_func <- function(ti)
eval_surv_base_fun(
ti = ti, omega = c(fixef_vary_surv, associations), b_func = expansion,
gl_dat = gl_dat, delta = log_haz_offset)
# simulate the recurrent events one at a time
max_sample <- 10L
Z <- matrix(rep(Z, each = max_sample), max_sample)
event <- y <- lf_trunc <- rep(NA_real_, max_sample)
lf_trunc_i <- 0
left_trunc_surv <- 1
for(i in 1:max_sample){
# sample a random uniform variable and invert the survival function
rng_i <- runif(1)
lf_trunc[i] <- lf_trunc_i
root_func <- function(x) rng_i - surv_func(x) / left_trunc_surv
if(root_func(2) < 0){
# the observation is right-censored and we can exit
y[i] <- 2
event[i] <- 0
break
}
# we need to invert the survival function to find the observation time
root <- uniroot(root_func, c(lf_trunc_i, 2), tol = 1e-6)
lf_trunc_i <- y[i] <- root$root
event[i] <- 1
left_trunc_surv <- surv_func(y[i])
}
# keep the needed data and return
colnames(Z) <- paste0("Z_", 1:NCOL(Z) - 1L)
surv_data <- cbind(lf_trunc = lf_trunc[1:i], y = y[1:i], event = event[1:i],
Z[1:i, -1, drop = FALSE], id = id)
list(marker_data = marker_data, surv_data = surv_data)
})
# combine the data and return
marker_data <- as.data.frame(do.call(
rbind, lapply(dat, `[[`, "marker_data")))
marker_data$id <- as.integer(marker_data$id)
# the order does not matter
marker_data <- marker_data[sample.int(NROW(marker_data)), ]
surv_data <- as.data.frame(do.call(
rbind, lapply(dat, `[[`, "surv_data")))
surv_data$id <- as.integer(surv_data$id)
# the order does not matter
surv_data <- surv_data[sample.int(NROW(surv_data)), ]
list(marker_data = marker_data, surv_data = surv_data)
}
# sample a moderately large data set
set.seed(2)
dat <- sim_dat(1000L)
# we show a few properties of the data below
NROW(dat$marker_data) # number of observed marker
#> [1] 4463
NROW(dat$surv_data) # number of survival outcomes
#> [1] 2495
sum(dat$surv_data$event) # number of observed events
#> [1] 1500
# the data for one individual
subset(dat$marker_data, id == 1)
#> Y1 Y2 X1_1 time id
#> 5 NA 4.3373 -0.23970 1.8877 1
#> 4 NA 4.8558 0.70795 1.8869 1
#> 3 NA 3.1915 0.13242 1.4047 1
#> 2 1.399 1.7030 -0.08025 1.1467 1
#> 1 -1.286 0.4379 -1.13038 0.3361 1
subset(dat$surv_data, id == 1)
#> lf_trunc y event Z_1 id
#> 2 1.091 1.205 1 -0.9706 1
#> 3 1.205 2.000 0 -0.9706 1
#> 1 0.000 1.091 1 -0.9706 1
# distribution of observed events per individual
proportions(table(table(dat$surv_data$id) - 1L))
#>
#> 0 1 2 3 4 5 6 7 8 9
#> 0.320 0.320 0.154 0.090 0.057 0.031 0.008 0.005 0.004 0.011
# estimate the model with this package. Get the object we need for the
# optimization
marker_1 <- marker_term(
Y1 ~ X1_1, id = id, subset(dat$marker_data, !is.na(Y1)),
time_fixef = ns_term(time, knots = c(.5, 1.5), Boundary.knots = c(0, 2)),
time_rng = ns_term(time, knots = numeric(), Boundary.knots = c(0, 2),
intercept = TRUE))
marker_2 <- marker_term(
Y2 ~ 1, id = id, subset(dat$marker_data, !is.na(Y2)),
time_fixef = poly_term(time, degree = 2, raw = TRUE),
time_rng = poly_term(time, degree = 1, raw = TRUE, intercept = TRUE))
surv_obj <- surv_term(
Surv(lf_trunc, y, event) ~ Z_1, id = id, dat$surv_data,
# we select a more flexible model for the baseline hazard
time_fixef = poly_term(y, degree = 3L, raw = TRUE, use_log = TRUE),
with_frailty = TRUE)
comp_obj <- joint_ms_ptr(markers = list(marker_1, marker_2),
survival_terms = surv_obj, max_threads = 4L)
rm(marker_1, marker_2, surv_obj)
# get the starting values
system.time(start_val <- joint_ms_start_val(comp_obj, gr_tol = .1))
#> user system elapsed
#> 3.960 0.004 1.272
# lower bound at the starting values
print(-attr(start_val, "value"), digits = 8)
#> [1] -7542.5634
# check that the gradient is correct
f <- function(x){
start_val[seq_along(x)] <- x
joint_ms_lb(comp_obj, start_val)
}
all.equal(numDeriv::grad(f, head(start_val, 29 + 2 * 20)),
head(joint_ms_lb_gr(comp_obj, start_val), 29 + 2 * 20))
#> [1] "Mean relative difference: 6.815e-08"
# find the maximum lower bound estimate
system.time(opt_out <- joint_ms_opt(comp_obj, par = start_val, max_it = 1000L,
pre_method = 3L, cg_tol = .2, c2 = .1,
gr_tol = .1))
#> user system elapsed
#> 16.794 0.000 4.201
# we set gr_tol in the call so this is the convergence criterion for the
# gradient
sqrt(sum(joint_ms_lb_gr(comp_obj, opt_out$par)^2))
#> [1] 0.09355
opt_out$info # convergence code (0 == 'OK')
#> [1] 0
print(-opt_out$value, digits = 8) # maximum lower bound value
#> [1] -7257.0219
opt_out$counts
#> function gradient n_cg
#> 552 398 1464
# find the maximum lower bound with lbfgs
library(lbfgsb3c)
system.time(lbfgs_res <- lbfgsb3c(
start_val, function(x) joint_ms_lb(comp_obj, x),
function(x) joint_ms_lb_gr(comp_obj, x),
control = list(factr = 1e-8, maxit = 2000L)))
#> user system elapsed
#> 64.972 0.011 16.249
lbfgs_res$convergence # convergence code (0 == 'OK')
#> [1] 1
print(-lbfgs_res$value, digits = 8) # maximum lower bound value
#> [1] -7257.4011
lbfgs_res$counts # may have hit maxit!
#> [1] 2000 2000
# compare the estimates with the actual values. Start with the fixed effects
fmt_ests <- joint_ms_format(comp_obj, opt_out$par)
# the parameters for the first marker
fmt_ests$markers[[1]]
#> $fixef
#>
#> -0.4924 1.0171
#>
#> $fixef_vary
#>
#> 1.403 -1.191 -2.031
fixef_marker[[1]] # true values
#> [1] -0.5 1.0
fixef_vary_marker[[1]] # true values
#> [1] 1.4 -1.2 -2.1
# the parameters for the second marker
fmt_ests$markers[[2]]
#> $fixef
#>
#> 0.244
#>
#> $fixef_vary
#>
#> 0.6132 0.6280
fixef_marker[[2]] # true values
#> [1] 0.25
fixef_vary_marker[[2]] # true values
#> [1] 0.50 0.67
# the fixed effects for the survival outcome and the association parameters
fmt_ests$survival[[1]]
#> $fixef
#>
#> -0.5069 0.2612
#>
#> $fixef_vary
#>
#> 0.40372 0.07559 0.03288
#>
#> $associations
#>
#> -0.7302 0.7050
fixef_surv
#> [1] -0.50 0.25
fixef_vary_surv
#> [1] 0.5
associations
#> [1] -0.8 0.7
# the parameters for covariance matrix of the random effects
fmt_ests$vcov$vcov_vary
#> [,1] [,2] [,3] [,4]
#> [1,] 0.316844 0.008855 -0.04177 0.012386
#> [2,] 0.008855 0.167055 -0.05167 -0.004318
#> [3,] -0.041771 -0.051666 0.30518 0.090893
#> [4,] 0.012386 -0.004318 0.09089 0.139183
vcov_vary # the true values
#> [,1] [,2] [,3] [,4]
#> [1,] 0.35 0.02 -0.05 0.01
#> [2,] 0.02 0.12 -0.06 -0.01
#> [3,] -0.05 -0.06 0.32 0.09
#> [4,] 0.01 -0.01 0.09 0.12
# the parameters for the error term covariance matrix
fmt_ests$vcov$vcov_marker
#> [,1] [,2]
#> [1,] 0.37499 0.09843
#> [2,] 0.09843 0.16336
vcov_marker
#> [,1] [,2]
#> [1,] 0.36 0.10
#> [2,] 0.10 0.16
# the parameters for the frailty covariance matrix
fmt_ests$vcov$vcov_surv
#> [,1]
#> [1,] 0.03771
vcov_surv
#> [,1]
#> [1,] 0.04Observed Information Matrix and Approximate Wald Intervals
The package supplies an approximation of the observed information matrix and the full Hessian of all the parameters. This is illustrated below.
# compute the Hessian
system.time(hess <- joint_ms_hess(comp_obj, par = opt_out$par))
#> user system elapsed
#> 7.693 0.032 7.725
dim(hess$hessian_all) # the full matrix!
#> [1] 20029 20029
# compare parts it with those from numerical differentiation from R
n_comp <- 150L
hess_num <- numDeriv::jacobian(
function(x){
par <- opt_out$par
par[1:n_comp] <- x
joint_ms_lb_gr(comp_obj, par)[1:n_comp]
}, head(opt_out$par, n_comp))
# did we get the same?
all.equal(hess_num, as.matrix(hess$hessian_all[1:n_comp, 1:n_comp]),
check.attributes = FALSE)
#> [1] TRUE
# compute the covariance matrix from the approximate observed information matrix
# the part only for the model parameters accounting for the variational
# parameters
dim(hess$hessian)
#> [1] 29 29
obs_mat <- -hess$hessian # the observed information matrix
SEs <- sqrt(diag(solve(-obs_mat))) # approximate standard errors
# Wald type confidence intervals for the association parameters
idx_assoc <- comp_obj$indices$survival[[1]]$associations
matrix(opt_out$par[idx_assoc] + outer(SEs[idx_assoc], c(-1, 1) * 1.96), 2,
dimnames = list(c("alpha_1", "alpha_2"),
sprintf("%.2f pct", c(2.5, 97.5))))
#> 2.50 pct 97.50 pct
#> alpha_1 -1.0617 -0.3986
#> alpha_2 0.6314 0.7785
# we can do this for all parameters. We illustrate this by showing the
# estimates along with standard errors
do_comb <- function(est, se){
if(is.list(est))
setNames(mapply(do_comb, est, se, SIMPLIFY = FALSE), names(est))
else
rbind(Estimates = est, SE = se)
}
par_fmt <- joint_ms_format(comp_obj, opt_out$par)
par_fmt_SE <- joint_ms_format(comp_obj, SEs)
est_n_se <-
do_comb(par_fmt[c("markers", "survival")], par_fmt_SE[c("markers", "survival")])
# show the result
est_n_se$markers[[1]]
#> $fixef
#>
#> Estimates -0.49242 1.01709
#> SE 0.04692 0.01205
#>
#> $fixef_vary
#>
#> Estimates 1.40273 -1.1905 -2.0313
#> SE 0.05532 0.1179 0.0451
fixef_marker[[1]]
#> [1] -0.5 1.0
fixef_vary_marker[[1]]
#> [1] 1.4 -1.2 -2.1
est_n_se$markers[[2]]
#> $fixef
#>
#> Estimates 0.24400
#> SE 0.03147
#>
#> $fixef_vary
#>
#> Estimates 0.6132 0.62801
#> SE 0.0616 0.02978
fixef_marker[[2]]
#> [1] 0.25
fixef_vary_marker[[2]]
#> [1] 0.50 0.67
est_n_se$survival[[1]]
#> $fixef
#>
#> Estimates -0.50691 0.26119
#> SE 0.04753 0.05065
#>
#> $fixef_vary
#>
#> Estimates 0.40372 0.07559 0.03288
#> SE 0.05953 0.08393 0.02596
#>
#> $associations
#>
#> Estimates -0.7302 0.70498
#> SE 0.1692 0.03752
fixef_surv
#> [1] -0.50 0.25
fixef_vary_surv
#> [1] 0.5
associations
#> [1] -0.8 0.7We can compute standard error estimates for the covariance matrices’ parameters but this requires an application of the delta method because of the parameterization of the covariance matrices.
Approximate Likelihood Ratio based Confidence Intervals
Assuming that the lower bound is tight, we can construct approximate likelihood ratio based confidence intervals using the lower bound. We show how to do this below with the mask argument of joint_ms_opt.
# fixed input
level <- .95 # confidence level
which_fix <- 14L
# see the indices element for which element is fixed. It is an association
# parameter in this case
comp_obj$indices$survival[[1]]$associations # confidence interval for
#> [1] 14 15
# assumed plausible values of the parameter. The joint_ms_profile function
# shown later finds these automatically
params <- opt_out$par[which_fix]
params <- seq(params - .5, params + .5, length.out = 15)
# find the maximum lower bound values
lbs_max <- sapply(params, function(x){
par <- opt_out$par
par[which_fix] <- x
res <- joint_ms_opt(comp_obj, par = par, max_it = 1000L,
pre_method = 3L, cg_tol = .2, c2 = .1,
# -1 needed in the psqn package (zero-based indices)
mask = which_fix - 1L)
# return the maximum lower bound if the method converged
if(!res$convergence)
stop("did not converge")
-res$value
})# find the critical value and the approximate confidence interval
z_vals <- sqrt(pmax(-opt_out$value - lbs_max, 0) * 2)
z_vals <- ifelse(params < opt_out$par[which_fix], -1, 1) * z_vals
alpha <- 1 - level
pvs <- c(alpha / 2, 1 - alpha/2)
conf_int <- setNames(approx(z_vals, params, xout = qnorm(pvs))$y,
sprintf("%.2f pct.", 100 * pvs))
conf_int # the approximate confidence interval
#> 2.50 pct. 97.50 pct.
#> -1.0785 -0.4066
# plot the approximate log profile likelihood and highlight the critical value
par(mar = c(5, 5, 1, 1))
plot(params, lbs_max, pch = 16, bty = "l", xlab = "Association parameter",
ylab = "Approximate log profile likelihood")
grid()
smooth_est <- smooth.spline(params, lbs_max)
lines(predict(smooth_est, seq(min(params), max(params), length.out = 100)))
# mark the confidence interval
abline(h = -opt_out$value - qchisq(level, 1) / 2, lty = 2)
abline(v = conf_int, lty = 2)
The joint_ms_profile uses similar steps to the above to find approximate profile likelihood based confidence intervals. An example is shown below.
# construct the approximate likelihood ratio based confidence interval
system.time(
prof_conf <- joint_ms_profile(
comp_obj, opt_out = opt_out, which_prof = which_fix, delta = .25,
level = level,
max_it = 1000L, pre_method = 3L, cg_tol = .2, c2 = .1))
#>
#> Finding the upper limit of the approximate profile likelihood curve
#> LogLike: -7258.1624 at -0.480153
#> LogLike: -7261.6286 at -0.230153
#> LogLike: -7259.3813 at -0.371893. Lb, target, ub: -7259.3813, -7258.9426, -7258.1624
#> LogLike: -7258.8957 at -0.410649. Lb, target, ub: -7259.3813, -7258.9426, -7258.8957
#>
#> Finding the lower limit of the approximate profile likelihood curve
#> LogLike: -7258.0417 at -0.980153
#> LogLike: -7260.7667 at -1.230153
#> LogLike: -7259.3134 at -1.108822. Lb, target, ub: -7259.3134, -7258.9426, -7258.0417
#> LogLike: -7258.8517 at -1.067539. Lb, target, ub: -7259.3134, -7258.9426, -7258.8517
#> LogLike: -7257.0219 at -0.730153
#> user system elapsed
#> 89.602 0.016 22.409prof_conf$confs # the approximate confidence interval
#> 2.50 pct. 97.50 pct.
#> -1.0760 -0.4067
# plot the approximate log profile likelihood and highlight the critical value
par(mar = c(5, 5, 1, 1))
with(prof_conf, {
plot(xs, p_log_Lik, pch = 16, bty = "l",
xlab = "Association parameter", ylab = "Approximate log profile likelihood")
grid()
smooth_est <- smooth.spline(xs, p_log_Lik)
lines(predict(smooth_est, seq(min(xs), max(xs), length.out = 100)))
abline(v = confs, lty = 2)
})
The Hessian can be used to possibly decrease the computation time by passing it to the hess argument as illustrated below. Starting values along the profile likelihood curve is then found using an approximate normal distribution.
# construct the approximate likelihood ratio based confidence interval
system.time(
prof_conf_fast <- joint_ms_profile(
comp_obj, opt_out = opt_out, which_prof = which_fix,
# we can use the standard error
delta = 2.5 * sqrt(diag(solve(hess$hessian)))[which_fix],
level = level,
max_it = 1000L, pre_method = 3L, cg_tol = .2, c2 = .1,
hess = hess))
#> Computing the vector to find the starting values along the profile likelihood curve
#>
#> Finding the upper limit of the approximate profile likelihood curve
#> LogLike: -7260.3169 at -0.307229
#> LogLike: -7257.0219 at -0.730153
#> LogLike: -7257.8966 at -0.511013. Lb, target, ub: -7260.3169, -7258.9426, -7257.8966
#> LogLike: -7258.9272 at -0.407852. Lb, target, ub: -7260.3169, -7258.9426, -7258.9272
#> LogLike: -7259.0783 at -0.395511. Lb, target, ub: -7259.0783, -7258.9426, -7258.9272
#>
#> Finding the lower limit of the approximate profile likelihood curve
#> LogLike: -7259.8232 at -1.153077
#> LogLike: -7257.0219 at -0.730153
#> LogLike: -7258.1134 at -0.987911. Lb, target, ub: -7259.8232, -7258.9426, -7258.1134
#> LogLike: -7258.9025 at -1.072728. Lb, target, ub: -7259.8232, -7258.9426, -7258.9025
#> LogLike: -7259.0333 at -1.085092. Lb, target, ub: -7259.0333, -7258.9426, -7258.9025
#> LogLike: -7257.0219 at -0.730153
#> user system elapsed
#> 65.725 1.184 25.597# we got the same
prof_conf$confs
#> 2.50 pct. 97.50 pct.
#> -1.0760 -0.4067
prof_conf_fast$confs
#> 2.50 pct. 97.50 pct.
#> -1.0766 -0.4066Two Markers and a Recurrent Event Without Frailty
In this section, we estimate a model like before but without the frailty.
# set the frailty variance to ~0
vcov_surv <- matrix(1e-6^2, 1)
# plot the hazard with pointwise quantiles
par(mar = c(5, 5, 1, 1))
plot_surv(
time_fixef = b_term,
time_rng = list(
ns_term(knots = numeric(), Boundary.knots = c(0, 2), intercept = TRUE),
poly_term(degree = 1, raw = TRUE, intercept = TRUE)),
x_range = c(0, 2), vcov_vary = vcov_vary, frailty_var = vcov_surv,
ps = c(.25, .5, .75), log_hazard_shift = fixef_surv[1],
fixef_vary = fixef_vary_surv, associations = associations)
# sample a moderately large data set
set.seed(2)
dat <- sim_dat(1000L)
# we show a few properties of the data below
NROW(dat$marker_data) # number of observed marker
#> [1] 4480
NROW(dat$surv_data) # number of survival outcomes
#> [1] 2566
sum(dat$surv_data$event) # number of observed events
#> [1] 1571
# the data for one individual
subset(dat$marker_data, id == 1)
#> Y1 Y2 X1_1 time id
#> 5 NA 4.3373 -0.23970 1.8877 1
#> 4 NA 4.8558 0.70795 1.8869 1
#> 3 NA 3.1915 0.13242 1.4047 1
#> 2 1.399 1.7030 -0.08025 1.1467 1
#> 1 -1.286 0.4379 -1.13038 0.3361 1
subset(dat$surv_data, id == 1)
#> lf_trunc y event Z_1 id
#> 2 1.206 1.330 1 -0.9706 1
#> 3 1.330 2.000 0 -0.9706 1
#> 1 0.000 1.206 1 -0.9706 1
# distribution of observed events per individual
proportions(table(table(dat$surv_data$id) - 1L))
#>
#> 0 1 2 3 4 5 6 7 8 9
#> 0.317 0.279 0.184 0.102 0.057 0.029 0.010 0.007 0.004 0.011
# estimate the model with this package. Get the object we need for the
# optimization
marker_1 <- marker_term(
Y1 ~ X1_1, id = id, subset(dat$marker_data, !is.na(Y1)),
time_fixef = ns_term(time, knots = c(.5, 1.5), Boundary.knots = c(0, 2)),
time_rng = ns_term(time, knots = numeric(), Boundary.knots = c(0, 2),
intercept = TRUE))
marker_2 <- marker_term(
Y2 ~ 1, id = id, subset(dat$marker_data, !is.na(Y2)),
time_fixef = poly_term(time, degree = 2, raw = TRUE),
time_rng = poly_term(time, degree = 1, raw = TRUE, intercept = TRUE))
surv_obj <- surv_term(
Surv(lf_trunc, y, event) ~ Z_1, id = id, dat$surv_data,
# we select a more flexible model for the baseline hazard
time_fixef = poly_term(y, degree = 3L, raw = TRUE, use_log = TRUE))
comp_obj <- joint_ms_ptr(markers = list(marker_1, marker_2),
survival_terms = surv_obj, max_threads = 4L)
rm(marker_1, marker_2, surv_obj)
# get the starting values
system.time(start_val <- joint_ms_start_val(comp_obj))
#> user system elapsed
#> 2.596 0.000 0.750
# lower bound at the starting values
print(-attr(start_val, "value"), digits = 8)
#> [1] -7423.5688
# check that the gradient is correct
f <- function(x){
start_val[seq_along(x)] <- x
joint_ms_lb(comp_obj, start_val)
}
all.equal(numDeriv::grad(f, head(start_val, 28 + 2 * 14)),
head(joint_ms_lb_gr(comp_obj, start_val), 28 + 2 * 14))
#> [1] "Mean relative difference: 1.514e-07"
# find the maximum lower bound estimate
system.time(opt_out <- joint_ms_opt(comp_obj, par = start_val, max_it = 1000L,
pre_method = 3L, cg_tol = .2, c2 = .1,
rel_eps = 1e-10))
#> user system elapsed
#> 10.709 0.000 2.678
# check the gradient norm. We may need to reduce the convergence tolerance if
# this is not small. In can also be a sign of convergence issues
sqrt(sum(joint_ms_lb_gr(comp_obj, opt_out$par)^2))
#> [1] 0.01839
opt_out$info # convergence code (0 == 'OK')
#> [1] 0
print(-opt_out$value, digits = 8) # maximum lower bound value
#> [1] -7093.4009
opt_out$counts
#> function gradient n_cg
#> 412 308 685
# compare the estimates with the actual values. Start with the fixed effects
fmt_ests <- joint_ms_format(comp_obj, opt_out$par)
# the parameters for the first marker
fmt_ests$markers[[1]]
#> $fixef
#> [1] -0.5361 1.0058
#>
#> $fixef_vary
#> [1] 1.439 -1.134 -2.067
fixef_marker[[1]] # true values
#> [1] -0.5 1.0
fixef_vary_marker[[1]] # true values
#> [1] 1.4 -1.2 -2.1
# the parameters for the second marker
fmt_ests$markers[[2]]
#> $fixef
#> [1] 0.2043
#>
#> $fixef_vary
#> [1] 0.6436 0.6173
fixef_marker[[2]] # true values
#> [1] 0.25
fixef_vary_marker[[2]] # true values
#> [1] 0.50 0.67
# the fixed effects for the survival outcome and the association parameters
fmt_ests$survival[[1]]
#> $fixef
#> [1] -0.4367 0.2499
#>
#> $fixef_vary
#> [1] 0.426456 -0.006997 0.010090
#>
#> $associations
#> [1] -0.8852 0.7342
fixef_surv
#> [1] -0.50 0.25
fixef_vary_surv
#> [1] 0.5
associations
#> [1] -0.8 0.7
# the parameters for covariance matrix of the random effects
fmt_ests$vcov$vcov_vary
#> [,1] [,2] [,3] [,4]
#> [1,] 0.30519 0.01912 -0.06209 0.02126
#> [2,] 0.01912 0.08415 -0.03358 -0.02355
#> [3,] -0.06209 -0.03358 0.30457 0.08503
#> [4,] 0.02126 -0.02355 0.08503 0.14096
vcov_vary # the true values
#> [,1] [,2] [,3] [,4]
#> [1,] 0.35 0.02 -0.05 0.01
#> [2,] 0.02 0.12 -0.06 -0.01
#> [3,] -0.05 -0.06 0.32 0.09
#> [4,] 0.01 -0.01 0.09 0.12
# the parameters for the error term covariance matrix
fmt_ests$vcov$vcov_marker
#> [,1] [,2]
#> [1,] 0.36068 0.09824
#> [2,] 0.09824 0.15530
vcov_marker
#> [,1] [,2]
#> [1,] 0.36 0.10
#> [2,] 0.10 0.16
# the parameters for the frailty covariance matrix
fmt_ests$vcov$vcov_surv # NULL as there is not any
#> NULLUsing the Estimated Variational Parameters
The joint_ms_va_par function can be used to extract the estimated variational parameters for each individual. These can be used to get marginal estimates as the variational approximation for each individual is an approximation of the conditional distribution of the random effects given the observed data.
# get the variational parameters
va_par <- joint_ms_va_par(comp_obj, opt_out$par)
# it is a list with the estimated variational parameters for each individual
va_par[c("157", "345", "878")]
#> $`157`
#> $`157`$mean
#> [1] -0.47330 -0.02991 0.85787 0.37735
#>
#> $`157`$vcov
#> [,1] [,2] [,3] [,4]
#> [1,] 0.196318 -0.0094526 -0.011039 0.0255789
#> [2,] -0.009453 0.0673965 0.003971 0.0006195
#> [3,] -0.011039 0.0039709 0.098737 -0.0459746
#> [4,] 0.025579 0.0006195 -0.045975 0.0412363
#>
#>
#> $`345`
#> $`345`$mean
#> [1] -0.03225 -0.11438 0.92849 0.42211
#>
#> $`345`$vcov
#> [,1] [,2] [,3] [,4]
#> [1,] 0.184687 0.015121 -0.005131 0.032541
#> [2,] 0.015121 0.072400 -0.001714 0.001443
#> [3,] -0.005131 -0.001714 0.038312 -0.029648
#> [4,] 0.032541 0.001443 -0.029648 0.064868
#>
#>
#> $`878`
#> $`878`$mean
#> [1] 0.23917 0.03414 0.25500 0.14914
#>
#> $`878`$vcov
#> [,1] [,2] [,3] [,4]
#> [1,] 0.189760 2.005e-02 -0.004532 1.121e-02
#> [2,] 0.020050 7.357e-02 -0.004647 -5.375e-05
#> [3,] -0.004532 -4.647e-03 0.055990 -2.593e-02
#> [4,] 0.011207 -5.375e-05 -0.025934 4.042e-02
# they can be used to approximate marginal measures such as the mean curve for
# the individual along with 95% pointwise confidence intervals as shown below
# for where the true curve is. You can skip the code below and go the plots after
plot_mean_curve <- function(who){
# get the data for the individual
m_dat <- subset(dat$marker_data, id == who)
# plot the estimated mean population curve with pointwise confidence intervals
par(mfcol = c(1, 2), bty = "l")
fmt_ests <- joint_ms_format(comp_obj, opt_out$par)
plot_marker(
time_fixef = ns_term(
knots = c(.5, 1.5), Boundary.knots = c(0, 2)),
time_rng = ns_term(
knots = numeric(), Boundary.knots = c(0, 2), intercept = TRUE),
fixef_vary = fmt_ests$markers[[1]]$fixef_vary, x_range = c(0, 2),
vcov_vary = fmt_ests$vcov$vcov_vary[1:2, 1:2], ylab = "Marker 1")
# add the points for the individual
y_hat <- m_dat$Y1 - cbind(1, m_dat$X1_1) %*% fmt_ests$markers[[1]]$fixef
points(m_dat$time, y_hat, pch = 16)
# plot the estimated mean curve with pointwise confidence intervals
va_par_who <- va_par[[who]]
xs <- seq(0, 2, length.out = 100)
mea <- g_funcs[[1]](xs) %*% fmt_ests$markers[[1]]$fixef_vary +
m_funcs[[1]](xs) %*% va_par_who$mean[1:2]
lines(xs, mea, lty = 2)
se <- apply(m_funcs[[1]](xs), 1, function(x)
sqrt(x %*% va_par_who$vcov[1:2, 1:2] %*% x))
polygon(c(xs, rev(xs)), c(mea + 1.96 * se, rev(mea - 1.95 * se)), border = NA,
col = gray(0, .1))
# plot the estimated mean population curve with pointwise confidence intervals
plot_marker(
time_fixef = poly_term(degree = 2, raw = TRUE),
time_rng = poly_term(degree = 1, raw = TRUE, intercept = TRUE),
fixef_vary = fmt_ests$markers[[2]]$fixef_vary, x_range = c(0, 2),
vcov_vary = fmt_ests$vcov$vcov_vary[3:4, 3:4], ylab = "Marker 2")
# add the points for the individual
y_hat <- m_dat$Y2 - fmt_ests$markers[[2]]$fixef
points(m_dat$time, y_hat, pch = 16)
# plot the estimated mean curve with confidence intervals
xs <- seq(0, 2, length.out = 100)
mea <- g_funcs[[2]](xs) %*% fmt_ests$markers[[2]]$fixef_vary +
m_funcs[[2]](xs) %*% va_par_who$mean[2:3]
lines(xs, mea, lty = 2)
se <- apply(m_funcs[[2]](xs), 1, function(x)
sqrt(x %*% va_par_who$vcov[3:4, 3:4] %*% x))
polygon(c(xs, rev(xs)), c(mea + 1.96 * se, rev(mea - 1.95 * se)), border = NA,
col = gray(0, .1))
}
# plot the curves for three individuals
plot_mean_curve(157)
plot_mean_curve(235)
plot_mean_curve(878)
Two Markers, the Observation Time Process, and a Terminal Event
We simulate and fit a model in this section where we have two markers, a recurrent event which is the observation times, and a terminal event.
# settings for the simulation
library(splines)
g_funcs <- list(
function(x)
ns(x, knots = c(3.33, 6.67), Boundary.knots = c(0, 10)),
function(x)
# a raw polynomial
outer(x, 1:2, `^`))
m_funcs <- list(
function(x)
ns(x, knots = numeric(), Boundary.knots = c(0, 10), intercept = TRUE),
function(x)
# a raw polynomial
outer(x, 0:1, `^`))
fixef_vary_marker <- list(c(1.4, 1.2, -2.1), c(.5, -.02)) # beta
fixef_marker <- list(c(-.5, 2), 1) # gamma
# Psi
vcov_vary <- structure(c(0.35, 0.08, -0.05, 0.01, 0.08, 1.92, -0.24, -0.04,
-0.05, -0.24, 0.32, 0.09, 0.01, -0.04, 0.09, 0.12),
.Dim = c(4L, 4L))
vcov_marker <- matrix(c(.6^2, .1, .1, .4^2), 2)
# plot the markers' mean curve
library(VAJointSurv)
par(mar = c(5, 5, 1, 1))
plot_marker(
time_fixef = ns_term(
knots = c(3.33, 6.67), Boundary.knots = c(0, 10), intercept = FALSE),
time_rng = ns_term(
knots = numeric(), Boundary.knots = c(0, 10), intercept = TRUE),
fixef_vary = fixef_vary_marker[[1]], x_range = c(0, 10),
vcov_vary = vcov_vary[1:2, 1:2], ylab = "Marker 1")
plot_marker(
time_fixef = poly_term(degree = 2, raw = TRUE),
time_rng = poly_term(degree = 1, raw = TRUE, intercept = TRUE),
fixef_vary = fixef_vary_marker[[2]], x_range = c(0, 10),
vcov_vary = vcov_vary[3:4, 3:4], ylab = "Marker 2")
# the survival parameters
vcov_surv <- matrix(c(.2^2, .15^2, .15^2, .25^2), 2) # Xi
fixef_surv <- list(c(-1, .25), .2)
associations <- list(c(.6, -.4), c(-.7, .2))
fixef_vary_surv <- list(c(.5, .1, -.2, .11),
c(-1, -.25))
b_funcs <- list(
function(x) bs(x, knots = 5, Boundary.knots = c(0, 10)),
function(x) ns(x, knots = 5, Boundary.knots = c(0, 10)))
# plot the log hazard with the 25%, 50% and 75% quantiles
library(SimSurvNMarker)
plot_surv(
time_fixef = bs_term(knots = 5, Boundary.knots = c(0, 10)),
time_rng = list(
ns_term(knots = numeric(), Boundary.knots = c(0, 10), intercept = TRUE),
poly_term(degree = 1, raw = TRUE, intercept = TRUE)),
x_range = c(0, 10), fixef_vary = fixef_vary_surv[[1]],
vcov_vary = vcov_vary, frailty_var = vcov_surv[1, 1], ps = c(.25, .5, .75),
associations = associations[[1]], log_hazard_shift = fixef_surv[[1]][1],
ylab = "Terminal event")
plot_surv(
time_fixef = ns_term(knots = 5, Boundary.knots = c(0, 10)),
time_rng = list(
ns_term(knots = numeric(), Boundary.knots = c(0, 10), intercept = TRUE),
poly_term(degree = 1, raw = TRUE, intercept = TRUE)),
x_range = c(0, 10), fixef_vary = fixef_vary_surv[[2]],
vcov_vary = vcov_vary, frailty_var = vcov_surv[2, 2], ps = c(.25, .5, .75),
associations = associations[[2]], log_hazard_shift = fixef_surv[[2]][1],
ylab = "Observation process")
library(mvtnorm)
# simulates from the model by sampling a given number of individuals
sim_dat <- function(n_ids){
# simulate the outcomes
gl_dat <- get_gl_rule(100L)
dat <- lapply(1:n_ids, function(id){
# sample the terminal event time and the censoring time
cens <- min(rexp(1, rate = 1/10), 10)
U <- drop(rmvnorm(1, sigma = vcov_vary))
frailties <- drop(rmvnorm(1, sigma = vcov_surv))
Z1 <- c(1, runif(1, -1, 1))
log_haz_offset <- sum(Z1 * fixef_surv[[1]]) + frailties[1]
# assign the conditional survival function
expansion <- function(x, b_func)
cbind(b_func(x), m_funcs[[1]](x) %*% U[1:2],
m_funcs[[2]](x) %*% U[3:4])
surv_func <- function(ti, fixef_vary_surv, associations, b_func){
formals(expansion)$b_func <- b_func
eval_surv_base_fun(
ti = ti, omega = c(fixef_vary_surv, associations), b_func = expansion,
gl_dat = gl_dat, delta = log_haz_offset)
}
# sample the survival time
rng <- runif(1)
root_func <- function(x, rng)
rng - surv_func(x, fixef_vary_surv = fixef_vary_surv[[1]],
associations = associations[[1]], b_func = b_funcs[[1]])
if(root_func(cens, rng) < 0){
# the observation is censored
y_terminal <- cens
event <- 0
} else {
# find the event time
root <- uniroot(root_func, c(0, cens), tol = 1e-6, rng = rng)
y_terminal <- root$root
event <- 1
}
terminal_outcome <- cbind(y = y_terminal, event = event, Z1 = Z1[2],
id = id)
# clean up
rm(list = setdiff(ls(), c("y_terminal", "terminal_outcome", "expansion",
"surv_func", "frailties", "U", "id")))
# simulate the observation times
Z2 <- 1
log_haz_offset <- sum(Z2 * fixef_surv[[2]]) + frailties[2]
root_func <- function(x, left_trunc_surv, rng)
rng - surv_func(x, fixef_vary_surv = fixef_vary_surv[[2]],
associations = associations[[2]], b_func = b_funcs[[2]]) /
left_trunc_surv
max_sample <- 1000L
left_trunc_surv <- 1
Z2 <- matrix(rep(Z2, each = max_sample), max_sample)
event <- y <- lf_trunc <- rep(NA_real_, max_sample)
lf_trunc_i <- 0
for(i in 1:max_sample){
# sample a random uniform variable and invert the survival function
rng_i <- runif(1)
lf_trunc[i] <- lf_trunc_i
if(root_func(y_terminal, left_trunc_surv, rng_i) < 0){
# the observation is right-censored and we can exit
y[i] <- y_terminal
event[i] <- 0
break
}
# we need to invert the survival function to find the observation time
root <- uniroot(root_func, c(lf_trunc_i, y_terminal), tol = 1e-6,
left_trunc_surv = left_trunc_surv, rng = rng_i)
lf_trunc_i <- y[i] <- root$root
event[i] <- 1
left_trunc_surv <- surv_func(
y[i], fixef_vary_surv = fixef_vary_surv[[2]], associations = associations[[2]],
b_func = b_funcs[[2]])
}
colnames(Z2) <- paste0("Z", 1:NCOL(Z2) - 1L)
obs_process <- cbind(lf_trunc = lf_trunc[1:i], y = y[1:i],
event = event[1:i], Z2[1:i, -1, drop = FALSE],
id = id)
# clean up
rm(list = setdiff(ls(), c("terminal_outcome", "U", "id",
"obs_process")))
# sample the fixed effect covariates
obs_time <- c(0, obs_process[obs_process[, "event"] == 1, "y"])
n_obs <- length(obs_time)
X1 <- cbind(1, rnorm(n_obs))
X2 <- matrix(1, n_obs)
colnames(X1) <- paste0("X1_", 1:NCOL(X1) - 1L)
colnames(X2) <- paste0("X2_", 1:NCOL(X2) - 1L)
X <- list(X1, X2)
# sample the outcomes
eta <- sapply(1:2, function(i)
X[[i]] %*% fixef_marker[[i]] +
drop(g_funcs[[i]](obs_time)) %*% fixef_vary_marker[[i]] +
drop(m_funcs[[i]](obs_time)) %*% U[1:2 + (i == 2) * 2])
ys <- eta + rmvnorm(n_obs, sigma = vcov_marker)
colnames(ys) <- paste0("Y", 1:2)
# mask some observations
do_mask <- sample.int(3L, n_obs, replace = TRUE)
ys[do_mask == 2, 1] <- NA
ys[do_mask == 3, 2] <- NA
X <- do.call(cbind, lapply(X, function(x) x[, -1, drop = FALSE]))
marker_data <- cbind(ys, X, time = obs_time, id = id)
return(list(marker_data = marker_data, obs_process = obs_process,
terminal_outcome = terminal_outcome))
})
# combine the data and return
marker_data <- as.data.frame(do.call(
rbind, lapply(dat, `[[`, "marker_data")))
marker_data$id <- as.integer(marker_data$id)
# the order does not matter
marker_data <- marker_data[sample.int(NROW(marker_data)), ]
obs_process <- as.data.frame(do.call(
rbind, lapply(dat, `[[`, "obs_process")))
obs_process$id <- as.integer(obs_process$id)
# the order does not matter
obs_process <- obs_process[sample.int(NROW(obs_process)), ]
terminal_outcome <- as.data.frame(do.call(
rbind, lapply(dat, `[[`, "terminal_outcome")))
terminal_outcome$id <- as.integer(terminal_outcome$id)
# the order does not matter
terminal_outcome <- terminal_outcome[sample.int(NROW(terminal_outcome)), ]
list(marker_data = marker_data, obs_process = obs_process,
terminal_outcome = terminal_outcome)
}
# sample a moderately large data set
set.seed(2)
dat <- sim_dat(1000L)
# we show a few properties of the data below
mean(dat$terminal_outcome$event) # mean event rate
#> [1] 0.782
sum(dat$obs_process$event) # number of observed markers less the individuals
#> [1] 2421
NROW(dat$marker_data) # number of observed markers less the individuals
#> [1] 3421
# distribution of observed marker per individual
proportions(table(table(dat$obs_process$id)))
#>
#> 1 2 3 4 5 6 7 8 9 10 11 12 13
#> 0.353 0.194 0.135 0.076 0.069 0.045 0.030 0.024 0.023 0.010 0.005 0.006 0.008
#> 14 15 16 17 18 19 20 21 22 23 24 25 26
#> 0.003 0.001 0.004 0.001 0.003 0.001 0.001 0.001 0.003 0.001 0.001 0.001 0.001
# show data for one individual
subset(dat$marker_data, id == 1)
#> Y1 Y2 X1_1 time id
#> X.3 NA 4.041 -0.02815 2.4032 1
#> X.2 NA 2.164 -1.03429 0.2988 1
#> X.1 -0.3812 2.646 -0.28571 0.2603 1
#> X -2.6308 1.584 -1.43968 0.0000 1
subset(dat$obs_process, id == 1)
#> lf_trunc y event id
#> 2 0.2603 0.2988 1 1
#> 1 0.0000 0.2603 1 1
#> 4 2.4032 2.8900 0 1
#> 3 0.2988 2.4032 1 1
subset(dat$terminal_outcome, id == 1)
#> y event Z1 id
#> 1 2.89 1 -0.6384 1
# estimate the model with this package. Get the object we need for the
# optimization
marker_1 <- marker_term(
Y1 ~ X1_1, id = id, subset(dat$marker_data, !is.na(Y1)),
time_fixef = ns_term(time, knots = c(3.33, 6.67), Boundary.knots = c(0, 10)),
time_rng = ns_term(time, knots = numeric(), Boundary.knots = c(0, 10),
intercept = TRUE))
marker_2 <- marker_term(
Y2 ~ 1, id = id, subset(dat$marker_data, !is.na(Y2)),
time_fixef = poly_term(time, degree = 2, raw = TRUE),
time_rng = poly_term(time, degree = 1, raw = TRUE, intercept = TRUE))
library(survival)
surv_terminal <- surv_term(
Surv(y, event) ~ Z1, id = id, dat$terminal_outcome,
time_fixef = bs_term(y, knots = 5, Boundary.knots = c(0, 10)),
with_frailty = TRUE)
surv_obs <- surv_term(
Surv(lf_trunc, y, event) ~ 1, id = id, dat$obs_process,
time_fixef = ns_term(y, knots = 5, Boundary.knots = c(0, 10)),
with_frailty = TRUE)
comp_obj <- joint_ms_ptr(markers = list(marker_1, marker_2),
survival_terms = list(surv_terminal, surv_obs),
max_threads = 4L)
rm(marker_1, marker_2, surv_terminal, surv_obs)
# get the starting values
system.time(start_val <- joint_ms_start_val(comp_obj, gr_tol = .1))
#> user system elapsed
#> 31.820 0.027 8.329
# lower bound at the starting values
print(-attr(start_val, "value"), digits = 8)
#> [1] -7885.9054
# check that the gradient is correct
f <- function(x){
start_val[seq_along(x)] <- x
joint_ms_lb(comp_obj, start_val)
}
all.equal(numDeriv::grad(f, head(start_val, 37 + 2 * 27)),
head(joint_ms_lb_gr(comp_obj, start_val), 37 + 2 * 27),
tolerance = 1e-6)
#> [1] TRUE
# find the maximum lower bound estimate
system.time(opt_out <- joint_ms_opt(comp_obj, par = start_val, max_it = 2000L,
pre_method = 3L, cg_tol = .2, c2 = .1,
gr_tol = .1))
#> user system elapsed
#> 241.003 0.112 60.287
# we set gr_tol in the call so this is the convergence criterion for the
# gradient
sqrt(sum(joint_ms_lb_gr(comp_obj, opt_out$par)^2))
#> [1] 0.09475
opt_out$info # convergence code (0 == 'OK')
#> [1] 0
print(-opt_out$value, digits = 8) # maximum lower bound value
#> [1] -7600.1724
opt_out$counts
#> function gradient n_cg
#> 3748 2665 19654
# compare the estimates with the actual values. Start with the fixed effects
fmt_ests <- joint_ms_format(comp_obj, opt_out$par)
# the parameters for the first marker
fmt_ests$markers[[1]]
#> $fixef
#>
#> -0.5136 1.9734
#>
#> $fixef_vary
#>
#> 1.400 1.561 -1.606
fixef_marker[[1]] # true values
#> [1] -0.5 2.0
fixef_vary_marker[[1]] # true values
#> [1] 1.4 1.2 -2.1
# the parameters for the second marker
fmt_ests$markers[[2]]
#> $fixef
#>
#> 1.036
#>
#> $fixef_vary
#>
#> 0.48410 -0.02007
fixef_marker[[2]] # true values
#> [1] 1
fixef_vary_marker[[2]] # true values
#> [1] 0.50 -0.02
# the fixed effects for the survival outcome and the association parameters
# for the terminal event
fmt_ests$survival[[1]]
#> $fixef
#>
#> -0.8571 0.2568
#>
#> $fixef_vary
#>
#> -0.07556 1.17483 -1.33168 0.31447
#>
#> $associations
#>
#> 0.5676 -0.3621
fixef_surv[[1]]
#> [1] -1.00 0.25
fixef_vary_surv[[1]]
#> [1] 0.50 0.10 -0.20 0.11
associations[[1]]
#> [1] 0.6 -0.4
# same for the observation process
fmt_ests$survival[[2]]
#> $fixef
#>
#> 0.2476
#>
#> $fixef_vary
#>
#> -1.3520 -0.3493
#>
#> $associations
#>
#> -0.7369 0.2078
fixef_surv[[2]]
#> [1] 0.2
fixef_vary_surv[[2]]
#> [1] -1.00 -0.25
associations[[2]]
#> [1] -0.7 0.2
# the parameters for covariance matrix of the random effects
fmt_ests$vcov$vcov_vary
#> [,1] [,2] [,3] [,4]
#> [1,] 0.294322 -0.145392 0.05288 0.001029
#> [2,] -0.145392 1.481999 -0.12440 0.007732
#> [3,] 0.052880 -0.124399 0.27871 0.091861
#> [4,] 0.001029 0.007732 0.09186 0.116766
vcov_vary # the true values
#> [,1] [,2] [,3] [,4]
#> [1,] 0.35 0.08 -0.05 0.01
#> [2,] 0.08 1.92 -0.24 -0.04
#> [3,] -0.05 -0.24 0.32 0.09
#> [4,] 0.01 -0.04 0.09 0.12
# the parameters for the error term covariance matrix
fmt_ests$vcov$vcov_marker
#> [,1] [,2]
#> [1,] 0.3890 0.1104
#> [2,] 0.1104 0.1688
vcov_marker
#> [,1] [,2]
#> [1,] 0.36 0.10
#> [2,] 0.10 0.16
# the parameters for the frailty covariance matrix
fmt_ests$vcov$vcov_surv
#> [,1] [,2]
#> [1,] 0.00693 -0.01237
#> [2,] -0.01237 0.02767
vcov_surv
#> [,1] [,2]
#> [1,] 0.0400 0.0225
#> [2,] 0.0225 0.0625Caching Expansions
Some basis expansions like bs_term and ns_term take relatively long time to evaluate in the approximation of the approximate expected cumulative hazard. Thus, it may be advantageous to save the expansions if the same quadrature rule is used. This is done by setting the cache_expansions argument to true. The pros of doing this is that the expensive basis expansions are only evaluated once which may decrees the computation time. The cons are
- We no longer use the CPU cache efficiently with present hardware. Thus, you may not see great advantages of using many threads and the performance
may even be worse. This is more likely to be an issue with larger data set. - It requires a lot more memory which may be an issue larger problems.
We illustrate this by showing the computation time where we change the number of threads we use.
# with caching
w_caching <- bench::mark(
`w/ caching 1 thread` =
joint_ms_lb(comp_obj, opt_out$par, n_threads = 1, cache_expansions = TRUE),
`w/ caching 2 thread` =
joint_ms_lb(comp_obj, opt_out$par, n_threads = 2, cache_expansions = TRUE),
`w/ caching 3 thread` =
joint_ms_lb(comp_obj, opt_out$par, n_threads = 3, cache_expansions = TRUE),
`w/ caching 4 thread` =
joint_ms_lb(comp_obj, opt_out$par, n_threads = 4, cache_expansions = TRUE))
w_caching[, c("expression", "median")]
#> # A tibble: 4 × 2
#> expression median
#> <bch:expr> <bch:tm>
#> 1 w/ caching 1 thread 7.29ms
#> 2 w/ caching 2 thread 3.93ms
#> 3 w/ caching 3 thread 2.66ms
#> 4 w/ caching 4 thread 2.17ms
# difference between one and four threads
with(w_caching, median[4] / median[1])
#> [1] 297ms
# w/o caching
wo_caching <- bench::mark(
`w/o caching 1 thread` =
joint_ms_lb(comp_obj, opt_out$par, n_threads = 1, cache_expansions = FALSE),
`w/o caching 2 thread` =
joint_ms_lb(comp_obj, opt_out$par, n_threads = 2, cache_expansions = FALSE),
`w/o caching 3 thread` =
joint_ms_lb(comp_obj, opt_out$par, n_threads = 3, cache_expansions = FALSE),
`w/o caching 4 thread` =
joint_ms_lb(comp_obj, opt_out$par, n_threads = 4, cache_expansions = FALSE))
wo_caching[, c("expression", "median")]
#> # A tibble: 4 × 2
#> expression median
#> <bch:expr> <bch:tm>
#> 1 w/o caching 1 thread 53.6ms
#> 2 w/o caching 2 thread 30.8ms
#> 3 w/o caching 3 thread 21.9ms
#> 4 w/o caching 4 thread 17.7ms
# difference between one and four threads
with(wo_caching, median[4] / median[1])
#> [1] 331msTwo Markers, the Observation Time Process, and a Terminal Event with Delayed Entry
We alter the previous example in this section by adding delayed entries. In this case, the likelihood has to be altered to account for the delayed entry as described by Crowther et al. (2016) and Berg and Drepper (2016).
To show the likelihood we use, assume that there is only one type of survival outcome, \(H = 1\). The log marginal likelihood of individual \(i\) where we do not properly account for the delayed entry is
\[\log E\left( f_i(\vec y_i\mid \vec U_i)h_{i1}(t_{i1}\mid \vec U_i)^{d_{i1}} \frac{S_{i1}(t_{i1}\mid \vec U_i)}{S_{i1}(v_{i1}\mid \vec U_i)}\right)\]
where \(\vec y_i\) is observed markers of individual \(i\), \(f_i\) is the conditional density for the markers of individual \(i\) given the random effects \(\vec U_i\), \(t_{i1}\) is the observed time, \(d_{i1}\) is an event indicator, \(h_{i1}\) is the conditional hazard, \(S_{i1}\) is the conditional survival function, and \(v_{i1}\in [0,\infty)\) is the left-truncation time. We approximate the log marginal likelihood term with a lower bound given by
\[\log E\left( f_i(\vec y_i\mid \vec U_i)h_{i1}(t_{i1}\mid \vec U_i)^{d_{i1}} \frac{S_{i1}(t_{i1}\mid \vec U_i)}{S_{i1}(v_{i1}\mid \vec U_i)}\right) \geq E_{Q_{i1}}\left(\log\left( f_i(\vec y_i\mid \vec U_i)h_{i1}(t_{i1}\mid \vec U_i)^{d_{i1}} \frac{S_{i1}(t_{i1}\mid \vec U_i)v(\vec U_i)} {S_{i1}(v_{i1}\mid \vec U_i)q_{\vec\theta_{i1}}(\vec U_i)} \right)\right)\]
where the expectation on the right-hand side is using the density \(q_{\vec\theta_{i1}}\) and \(v\) is the unconditional density of the random effects. Since this is a lower bound, we can jointly optimize the lower bound over the \(\vec\theta_{i1}\)s and the model parameters.
The correct likelihood with delayed entry is
\[\log E\left( f_i(\vec y_i\mid \vec U_i)h_{i1}(t_{i1}\mid \vec U_i)^{d_{i1}} S_{i1}(t_{i1}\mid \vec U_i)\right) - \log E\left(S_{i1}(v_{i1}\mid \vec U_i)\right).\]
Of course, we can replace this with an approximation in the form of
\[E_{Q_{i1}}\left(\log\left( f_i(\vec y_i\mid \vec U_i)h_{i1}(t_{i1}\mid \vec U_i)^{d_{i1}} S_{i1}(t_{i1}\mid \vec U_i)\frac{v(\vec U_i)}{q_{\theta_{i1}}(\vec U_i)} \right)\right) - E_{Q_{i2}}\left(\log\left(S_{i1}(v_{i1}\mid \vec U_i) \frac{v(\vec U_i)}{q_{\theta_{i2}}(\vec U_i)}\right)\right).\]
This has some implications:
- We still need to maximize over the \(\vec\theta_{i1}\)s and the model parameters. However, we have to minimize over the \(\vec\theta_{i2}\)s as they provide an upper bound on the log marginal likelihood. Thus, we have a maxmin problem rather than a pure maximization problem.
- The approximation is no longer guaranteed to be a lower bound on the log marginal likelihood because we are adding a lower bound on one term with an upper bound on the other term.
- The approximation of \(\log E\left(S_{i1}(v_{i1}\mid \vec U_i)\right)\) may be poor as we do not have the multivariate normal density from the markers.
Currently, we compute the terms in the marginal likelihood due to delayed entry with adaptive Gauss-Hermite quadrature. This can be quite slow if there are many random effects that needs to marginalized out. However, it seems that quite few quadrature nodes can be used to approximate the terms from delayed entry when the marginal survival probabilities are close to one (and the variance of the random effects are not extremely large). However, more quadrature nodes seems to be needed for the gradient of marginal likelihood terms due to delayed entry. Thus, you may experience odd behavior (very small improvements in the likelihood approximation at each step) of the optimizer if more quadrature nodes are needed.
We do believe that the maxmin problem mentioned in 1. can be solved efficiently but our preliminary “brute force” attempts were not fast.
Delayed entry is handled by the delayed argument of surv_term as illustrated in the example below.
library(mvtnorm)
# simulates from the model by sampling a given number of individuals
sim_dat <- function(n_ids){
# simulate the outcomes
gl_dat <- get_gl_rule(100L)
dat <- lapply(1:n_ids, function(id){
# sample the delayed entry time
delayed_entry <- pmax(runif(1, -.2, 1), 0)
# sample the censoring time
cens <- -Inf
while(cens < delayed_entry)
cens <- min(rexp(1, rate = 1/10), 10)
# sample the terminal event time and the random effects
Z1 <- c(1, runif(1, -1, 1))
y_terminal <- -Inf
while(y_terminal < delayed_entry){
U <- drop(rmvnorm(1, sigma = vcov_vary))
frailties <- drop(rmvnorm(1, sigma = vcov_surv))
log_haz_offset <- sum(Z1 * fixef_surv[[1]]) + frailties[1]
# assign the conditional survival function
expansion <- function(x, b_func)
cbind(b_func(x), m_funcs[[1]](x) %*% U[1:2],
m_funcs[[2]](x) %*% U[3:4])
surv_func <- function(ti, fixef_vary_surv, associations, b_func){
formals(expansion)$b_func <- b_func
eval_surv_base_fun(
ti = ti, omega = c(fixef_vary_surv, associations), b_func = expansion,
gl_dat = gl_dat, delta = log_haz_offset)
}
# sample the survival time
rng <- runif(1)
root_func <- function(x, rng)
rng - surv_func(x, fixef_vary_surv = fixef_vary_surv[[1]],
associations = associations[[1]], b_func = b_funcs[[1]])
if(root_func(cens, rng) < 0){
# the observation is censored
y_terminal <- cens
event <- 0
} else {
# find the event time
root <- uniroot(root_func, c(0, cens), tol = 1e-6, rng = rng)
y_terminal <- root$root
event <- 1
}
}
terminal_outcome <- cbind(y = y_terminal, event = event, Z1 = Z1[2],
id = id, delayed_entry = delayed_entry)
# clean up
rm(list = setdiff(ls(), c(
"y_terminal", "terminal_outcome", "expansion", "surv_func", "frailties",
"U", "id", "delayed_entry")))
# simulate the observation times
Z2 <- 1
log_haz_offset <- sum(Z2 * fixef_surv[[2]]) + frailties[2]
root_func <- function(x, left_trunc_surv, rng)
rng - surv_func(x, fixef_vary_surv = fixef_vary_surv[[2]],
associations = associations[[2]], b_func = b_funcs[[2]]) /
left_trunc_surv
max_sample <- 1000L
left_trunc_surv <- 1
Z2 <- matrix(rep(Z2, each = max_sample), max_sample)
event <- y <- lf_trunc <- rep(NA_real_, max_sample)
lf_trunc_i <- 0
for(i in 1:max_sample){
# sample a random uniform variable and invert the survival function
rng_i <- runif(1)
lf_trunc[i] <- lf_trunc_i
if(root_func(y_terminal, left_trunc_surv, rng_i) < 0){
# the observation is right-censored and we can exit
y[i] <- y_terminal
event[i] <- 0
break
}
# we need to invert the survival function to find the observation time
root <- uniroot(root_func, c(lf_trunc_i, y_terminal), tol = 1e-6,
left_trunc_surv = left_trunc_surv, rng = rng_i)
lf_trunc_i <- y[i] <- root$root
event[i] <- 1
left_trunc_surv <- surv_func(
y[i], fixef_vary_surv = fixef_vary_surv[[2]], associations = associations[[2]],
b_func = b_funcs[[2]])
}
colnames(Z2) <- paste0("Z", 1:NCOL(Z2) - 1L)
obs_process <- cbind(lf_trunc = lf_trunc[1:i], y = y[1:i],
event = event[1:i], Z2[1:i, -1, drop = FALSE],
id = id)
# account for the delayed entry
obs_process[, "lf_trunc"] <- pmax(delayed_entry, obs_process[, "lf_trunc"])
obs_process <- obs_process[
obs_process[, "y"] > delayed_entry, , drop = FALSE]
# clean up
rm(list = setdiff(ls(), c("terminal_outcome", "U", "id",
"obs_process", "delayed_entry")))
# sample the fixed effect covariates
obs_time <- c(delayed_entry, obs_process[obs_process[, "event"] == 1, "y"])
n_obs <- length(obs_time)
X1 <- cbind(1, rnorm(n_obs))
X2 <- matrix(1, n_obs)
colnames(X1) <- paste0("X1_", 1:NCOL(X1) - 1L)
colnames(X2) <- paste0("X2_", 1:NCOL(X2) - 1L)
X <- list(X1, X2)
# sample the outcomes
eta <- sapply(1:2, function(i)
X[[i]] %*% fixef_marker[[i]] +
drop(g_funcs[[i]](obs_time)) %*% fixef_vary_marker[[i]] +
drop(m_funcs[[i]](obs_time)) %*% U[1:2 + (i == 2) * 2])
ys <- eta + rmvnorm(n_obs, sigma = vcov_marker)
colnames(ys) <- paste0("Y", 1:2)
# mask some observations
do_mask <- sample.int(3L, n_obs, replace = TRUE)
ys[do_mask == 2, 1] <- NA
ys[do_mask == 3, 2] <- NA
X <- do.call(cbind, lapply(X, function(x) x[, -1, drop = FALSE]))
marker_data <- cbind(ys, X, time = obs_time, id = id)
return(list(marker_data = marker_data, obs_process = obs_process,
terminal_outcome = terminal_outcome))
})
# combine the data and return
marker_data <- as.data.frame(do.call(
rbind, lapply(dat, `[[`, "marker_data")))
marker_data$id <- as.integer(marker_data$id)
# the order does not matter
marker_data <- marker_data[sample.int(NROW(marker_data)), ]
obs_process <- as.data.frame(do.call(
rbind, lapply(dat, `[[`, "obs_process")))
obs_process$id <- as.integer(obs_process$id)
# the order does not matter
obs_process <- obs_process[sample.int(NROW(obs_process)), ]
terminal_outcome <- as.data.frame(do.call(
rbind, lapply(dat, `[[`, "terminal_outcome")))
terminal_outcome$id <- as.integer(terminal_outcome$id)
# the order does not matter
terminal_outcome <- terminal_outcome[sample.int(NROW(terminal_outcome)), ]
list(marker_data = marker_data, obs_process = obs_process,
terminal_outcome = terminal_outcome)
}
# sample a moderate sized data set
set.seed(2)
dat <- sim_dat(1000L)
# we show a few properties of the data below
mean(dat$terminal_outcome$event) # mean event rate
#> [1] 0.776
sum(dat$obs_process$event) # number of observed markers less the individuals
#> [1] 2620
NROW(dat$marker_data) # number of observed markers less the individuals
#> [1] 3620
# fraction of individuals with a delayed entry time
mean(dat$terminal_outcome$delayed_entry > 0)
#> [1] 0.821
# distribution of observed marker per individual
proportions(table(table(dat$obs_process$id)))
#>
#> 1 2 3 4 5 6 7 8 9 10 11 12 13
#> 0.345 0.214 0.130 0.074 0.052 0.050 0.025 0.021 0.019 0.019 0.010 0.008 0.004
#> 14 15 16 17 18 19 20 21 22 24 28 29 34
#> 0.002 0.003 0.003 0.004 0.003 0.001 0.003 0.002 0.002 0.001 0.001 0.001 0.001
#> 41 75
#> 0.001 0.001
# show data for one individual
subset(dat$marker_data, id == 1)
#> Y1 Y2 X1_1 time id
#> X NA 1.977 0.4177 0.02186 1
#> X.1 -0.9811 NA 0.9818 0.53440 1
#> X.2 -2.0826 2.277 -0.3927 0.86052 1
subset(dat$obs_process, id == 1)
#> lf_trunc y event id
#> 3 0.86052 1.1949 0 1
#> 1 0.02186 0.5344 1 1
#> 2 0.53440 0.8605 1 1
subset(dat$terminal_outcome, id == 1)
#> y event Z1 id delayed_entry
#> 1 1.195 1 0.1467 1 0.02186
# estimate the model with this package. Get the object we need for the
# optimization while NOT account and accounting for the delayed entry
marker_1 <- marker_term(
Y1 ~ X1_1, id = id, subset(dat$marker_data, !is.na(Y1)),
time_fixef = ns_term(time, knots = c(3.33, 6.67), Boundary.knots = c(0, 10)),
time_rng = ns_term(time, knots = numeric(), Boundary.knots = c(0, 10),
intercept = TRUE))
marker_2 <- marker_term(
Y2 ~ 1, id = id, subset(dat$marker_data, !is.na(Y2)),
time_fixef = poly_term(time, degree = 2, raw = TRUE),
time_rng = poly_term(time, degree = 1, raw = TRUE, intercept = TRUE))
# the wrong way
library(survival)
surv_terminal_wrong <- surv_term(
Surv(delayed_entry, y, event) ~ Z1, id = id, dat$terminal_outcome,
time_fixef = bs_term(y, knots = 5, Boundary.knots = c(0, 10)),
with_frailty = TRUE)
# the right way
surv_terminal <- surv_term(
Surv(delayed_entry, y, event) ~ Z1, id = id, dat$terminal_outcome,
time_fixef = bs_term(y, knots = 5, Boundary.knots = c(0, 10)),
with_frailty = TRUE,
# some have delayed entry
delayed = delayed_entry > 0)
surv_obs <- surv_term(
Surv(lf_trunc, y, event) ~ 1, id = id, dat$obs_process,
time_fixef = ns_term(y, knots = 5, Boundary.knots = c(0, 10)),
with_frailty = TRUE)
# the wrong way
comp_obj_wrong <- joint_ms_ptr(
markers = list(marker_1, marker_2),
survival_terms = list(surv_terminal_wrong, surv_obs),
max_threads = 4L)
# the right way
comp_obj <- joint_ms_ptr(
markers = list(marker_1, marker_2),
survival_terms = list(surv_terminal, surv_obs),
max_threads = 4L)
# the default number of Gauss-Hermite quadrature nodes we use
length(comp_obj$gh_quad_rule$node)
#> [1] 4
# we try with one fewer (has a big effect with this amount of random variables
# per individual)
ghq_fewer <- with(fastGHQuad::gaussHermiteData(3),
list(node = x, weight = w))
# get the starting values
system.time(start_val_wrong <- joint_ms_start_val(comp_obj_wrong, gr_tol = .1))
#> user system elapsed
#> 17.717 0.000 4.767
system.time(start_val <- joint_ms_start_val(comp_obj, gr_tol = .1))
#> user system elapsed
#> 25.646 0.004 6.745
system.time(start_val_few <- joint_ms_start_val(comp_obj, gr_tol = .1,
gh_quad_rule = ghq_fewer))
#> user system elapsed
#> 18.274 0.008 4.857
# lower bound at the starting values
print(-attr(start_val_wrong, "value"), digits = 8)
#> [1] -8370.427
print(-attr(start_val, "value"), digits = 8)
#> [1] -8370.4091
print(-attr(start_val_few, "value"), digits = 8)
#> [1] -8370.4091
# check that the gradient is correct
test_grad <- function(cmp, par, gh_quad_rule = cmp$gh_quad_rule){
f <- function(x){
par[seq_along(x)] <- x
joint_ms_lb(cmp, par, gh_quad_rule = gh_quad_rule)
}
all.equal(numDeriv::grad(f, head(par, 37 + 2 * 27)),
head(joint_ms_lb_gr(cmp, par, gh_quad_rule = gh_quad_rule),
37 + 2 * 27),
tolerance = 1e-6)
}
test_grad(comp_obj, start_val, comp_obj$gh_quad_rule)
#> [1] TRUE
test_grad(comp_obj, start_val_few, ghq_fewer) # less precise
#> [1] "Mean relative difference: 2.486e-06"
# find the maximum lower bound estimate
system.time(opt_out_wrong <- joint_ms_opt(
comp_obj_wrong, par = start_val_wrong, max_it = 2000L, pre_method = 3L,
cg_tol = .2, c2 = .1, gr_tol = .1))
#> user system elapsed
#> 95.22 0.02 23.81
# optimize in the right way with different number of Gauss-Hermite quadrature
# nodes
est_w_ghq <- function(par, gh_quad_rule = comp_obj$gh_quad_rule)
joint_ms_opt(
comp_obj, par = par, max_it = 2000L, pre_method = 3L,
# the function is more expensive. Thus, it makes more sense to get a better
# approximation in the conjugate gradient step
cg_tol = .1, c2 = .9,
# it takes much longer to get the same precision so we just stop when we do
# at a lower threshold
gr_tol = 1, gh_quad_rule = gh_quad_rule)
system.time(opt_out <- est_w_ghq(start_val))
#> user system elapsed
#> 1514.165 0.387 379.564
system.time(opt_out_fewer <- est_w_ghq(start_val_few, ghq_fewer))
#> user system elapsed
#> 345.31 0.14 86.44
# check the gradients again (expect to be somewhat off now)
test_grad(comp_obj, opt_out$par, comp_obj$gh_quad_rule)
#> [1] "Mean relative difference: 0.05749"
test_grad(comp_obj, opt_out_fewer$par, ghq_fewer) # less precise
#> [1] "Mean relative difference: 0.4548"
# we set gr_tol in some of the calls so this is the convergence criterion
# for the gradient in those cases
sqrt(sum(joint_ms_lb_gr(comp_obj_wrong, opt_out_wrong$par)^2))
#> [1] 0.08385
sqrt(sum(joint_ms_lb_gr(comp_obj, opt_out$par)^2))
#> [1] 0.791
sqrt(sum(joint_ms_lb_gr(comp_obj, opt_out_fewer$par,
gh_quad_rule = ghq_fewer)^2))
#> [1] 0.457
opt_out_wrong$info # convergence code (0 == 'OK')
#> [1] 0
opt_out$info # convergence code (0 == 'OK')
#> [1] 0
opt_out_fewer$info # convergence code (0 == 'OK')
#> [1] 0
print(-opt_out_wrong$value, digits = 8) # maximum lower bound value
#> [1] -7957.3085
print(-opt_out$value, digits = 8) # maximum lower bound value
#> [1] -7954.8461
print(-opt_out_fewer$value, digits = 8) # maximum lower bound value
#> [1] -7954.832
# the precision using fewer quadrature nodes may be a bit off. We evaluate the
# lower bound using more quadrature nodes to check this
print(-joint_ms_lb(comp_obj, opt_out_fewer$par), digits = 8)
#> [1] -7954.8565
opt_out_wrong$counts
#> function gradient n_cg
#> 1604 1166 5052
opt_out$counts
#> function gradient n_cg
#> 1299 1272 52601
opt_out_fewer$counts
#> function gradient n_cg
#> 742 715 25231
# compare the estimates with the actual values. Start with the fixed effects
fmt_ests_wrong <- joint_ms_format(comp_obj_wrong, opt_out_wrong$par)
fmt_ests <- joint_ms_format(comp_obj, opt_out$par)
fmt_ests_fewer <- joint_ms_format(comp_obj, opt_out_fewer$par)
# the parameters for the first marker
mapply(rbind, SIMPLIFY = FALSE,
wrong = fmt_ests_wrong$markers[[1]],
right = fmt_ests $markers[[1]],
`right (fewer)` = fmt_ests_fewer$markers[[1]])
#> $fixef
#>
#> wrong -0.5691 1.987
#> right -0.4965 1.987
#> right (fewer) -0.4966 1.987
#>
#> $fixef_vary
#>
#> wrong 1.675 1.732 -1.534
#> right 1.585 1.598 -1.652
#> right (fewer) 1.585 1.599 -1.653
fixef_marker[[1]] # true values
#> [1] -0.5 2.0
fixef_vary_marker[[1]] # true values
#> [1] 1.4 1.2 -2.1
# the parameters for the second marker
mapply(rbind, SIMPLIFY = FALSE,
wrong = fmt_ests_wrong$markers[[2]],
right = fmt_ests $markers[[2]],
`right (fewer)` = fmt_ests_fewer$markers[[2]])
#> $fixef
#>
#> wrong 1.034
#> right 1.023
#> right (fewer) 1.023
#>
#> $fixef_vary
#>
#> wrong 0.4476 -0.01808
#> right 0.4289 -0.01791
#> right (fewer) 0.4289 -0.01790
fixef_marker[[2]] # true values
#> [1] 1
fixef_vary_marker[[2]] # true values
#> [1] 0.50 -0.02
# the fixed effects for the survival outcome and the association parameters
# for the terminal event
mapply(rbind, SIMPLIFY = FALSE,
wrong = fmt_ests_wrong$survival[[1]],
right = fmt_ests $survival[[1]],
`right (fewer)` = fmt_ests_fewer$survival[[1]])
#> $fixef
#>
#> wrong -0.9226 0.1476
#> right -0.8757 0.1480
#> right (fewer) -0.8776 0.1480
#>
#> $fixef_vary
#>
#> wrong -0.052288 1.059 -0.8572 1.223
#> right -0.007732 1.024 -0.8176 1.228
#> right (fewer) -0.002374 1.021 -0.7997 1.227
#>
#> $associations
#>
#> wrong 0.6321 -0.4301
#> right 0.6251 -0.4346
#> right (fewer) 0.6271 -0.4348
fixef_surv[[1]]
#> [1] -1.00 0.25
fixef_vary_surv[[1]]
#> [1] 0.50 0.10 -0.20 0.11
associations[[1]]
#> [1] 0.6 -0.4
# same for the observation process
mapply(rbind, SIMPLIFY = FALSE,
wrong = fmt_ests_wrong$survival[[2]],
right = fmt_ests $survival[[2]],
`right (fewer)` = fmt_ests_fewer$survival[[2]])
#> $fixef
#>
#> wrong 0.2704
#> right 0.2241
#> right (fewer) 0.2238
#>
#> $fixef_vary
#>
#> wrong -1.499 -0.6836
#> right -1.430 -0.6212
#> right (fewer) -1.437 -0.6233
#>
#> $associations
#>
#> wrong -0.7097 0.2110
#> right -0.7062 0.2106
#> right (fewer) -0.7064 0.2109
fixef_surv[[2]]
#> [1] 0.2
fixef_vary_surv[[2]]
#> [1] -1.00 -0.25
associations[[2]]
#> [1] -0.7 0.2
# the parameters for covariance matrix of the random effects
fmt_ests_wrong$vcov$vcov_vary
#> [,1] [,2] [,3] [,4]
#> [1,] 0.4761 0.26833 -0.15198 -0.03440
#> [2,] 0.2683 2.24614 -0.37479 -0.04865
#> [3,] -0.1520 -0.37479 0.33351 0.09637
#> [4,] -0.0344 -0.04865 0.09637 0.14012
fmt_ests $vcov$vcov_vary
#> [,1] [,2] [,3] [,4]
#> [1,] 0.48951 0.24783 -0.15465 -0.03837
#> [2,] 0.24783 2.26102 -0.37253 -0.03968
#> [3,] -0.15465 -0.37253 0.33392 0.09793
#> [4,] -0.03837 -0.03968 0.09793 0.14206
fmt_ests_fewer$vcov$vcov_vary
#> [,1] [,2] [,3] [,4]
#> [1,] 0.48998 0.24214 -0.15456 -0.03844
#> [2,] 0.24214 2.24544 -0.37181 -0.03961
#> [3,] -0.15456 -0.37181 0.33366 0.09798
#> [4,] -0.03844 -0.03961 0.09798 0.14202
vcov_vary # the true values
#> [,1] [,2] [,3] [,4]
#> [1,] 0.35 0.08 -0.05 0.01
#> [2,] 0.08 1.92 -0.24 -0.04
#> [3,] -0.05 -0.24 0.32 0.09
#> [4,] 0.01 -0.04 0.09 0.12
norm(fmt_ests_wrong$vcov$vcov_vary - vcov_vary, "F")
#> [1] 0.5051
norm(fmt_ests $vcov$vcov_vary - vcov_vary, "F")
#> [1] 0.5046
norm(fmt_ests_fewer$vcov$vcov_vary - vcov_vary, "F")
#> [1] 0.4901
# the parameters for the error term covariance matrix
fmt_ests_wrong$vcov$vcov_marker
#> [,1] [,2]
#> [1,] 0.3728 0.1053
#> [2,] 0.1053 0.1648
fmt_ests $vcov$vcov_marker
#> [,1] [,2]
#> [1,] 0.3729 0.1051
#> [2,] 0.1051 0.1648
fmt_ests_fewer$vcov$vcov_marker
#> [,1] [,2]
#> [1,] 0.3729 0.1051
#> [2,] 0.1051 0.1648
vcov_marker
#> [,1] [,2]
#> [1,] 0.36 0.10
#> [2,] 0.10 0.16
norm(fmt_ests_wrong$vcov$vcov_marker - vcov_marker, "F")
#> [1] 0.01559
norm(fmt_ests $vcov$vcov_marker - vcov_marker, "F")
#> [1] 0.01551
norm(fmt_ests_fewer$vcov$vcov_marker - vcov_marker, "F")
#> [1] 0.01557
# the parameters for the frailty covariance matrix
fmt_ests_wrong$vcov$vcov_surv
#> [,1] [,2]
#> [1,] 0.05595 0.05547
#> [2,] 0.05547 0.06358
fmt_ests $vcov$vcov_surv
#> [,1] [,2]
#> [1,] 0.05308 0.05650
#> [2,] 0.05650 0.06311
fmt_ests_fewer$vcov$vcov_surv
#> [,1] [,2]
#> [1,] 0.05434 0.05447
#> [2,] 0.05447 0.06309
vcov_surv
#> [,1] [,2]
#> [1,] 0.0400 0.0225
#> [2,] 0.0225 0.0625
norm(fmt_ests_wrong$vcov$vcov_surv - vcov_surv, "F")
#> [1] 0.04929
norm(fmt_ests $vcov$vcov_surv - vcov_surv, "F")
#> [1] 0.04983
norm(fmt_ests_fewer$vcov$vcov_surv - vcov_surv, "F")
#> [1] 0.04743Two Markers, the Observation Time Process, and a Terminal Event with Time-varying Effects
We simulate from a model with both fixed and random time-varying covariate effects and non-proportional hazard effects of some of the covariates. This is possible within the package by using the weighted_term and stacked_term functions. See Note on Basis Expansions section for examples of how these functions work. We set up a simulation study below.
# settings for the simulation
library(splines)
g_funcs <- list(
function(x, data)
cbind(ns(x, knots = c(3.33, 6.67), Boundary.knots = c(0, 10)),
# a time-varying fixed effect
x * data$W1),
function(x, data)
cbind(outer(x, 1:2, `^`),
# a time-varying fixed effect
outer(x, 1:2, `^`) * data$W2))
m_funcs <- list(
function(x, data)
cbind(
# a random intercept
outer(x, 0, `^`),
# random slope
data$W1),
function(x, data)
cbind(outer(x, 0:1, `^`),
# a time-varying random effect
outer(x, 0:1, `^`) * data$W2))
fixef_vary_marker <- list(c(1.4, 1.2, -2.1, .1), c(.5, -.02, .25, -.1)) # beta
fixef_marker <- list(c(-.5, 2, 1), c(1, -1)) # gamma
# Psi
vcov_vary <-
structure(c(
0.622, -0.02, 0.13, -0.081, -0.103, 0.016, -0.02, 0.473, 0.121, 0.022, -0.032, 0.059, 0.13, 0.121, 0.801, -0.009, -0.193, 0, -0.081, 0.022, -0.009, 0.529, 0.04, 0.05, -0.103, -0.032, -0.193, 0.04, 0.636, -0.144, 0.016, 0.059, 0, 0.05, -0.144, 0.354), dim = c(6L, 6L))
vcov_marker <- matrix(c(.6^2, .1, .1, .4^2), 2)
# plot the markers' mean curve
library(VAJointSurv)
par(mar = c(5, 5, 1, 1))
rng_basis_one <- stacked_term(
poly_term(degree = 0, intercept = TRUE, raw = TRUE),
poly_term(degree = 0, intercept = TRUE, raw = TRUE) |>
weighted_term(W1))
show_marker_one <- \(newdata, ylab)
plot_marker(
time_fixef = stacked_term(
ns_term(
knots = c(3.33, 6.67), Boundary.knots = c(0, 10), intercept = FALSE),
poly_term(degree = 1, raw = TRUE) |>
weighted_term(W1)),
time_rng = rng_basis_one,
fixef_vary = fixef_vary_marker[[1]], x_range = c(0, 10),
vcov_vary = vcov_vary[1:2, 1:2], ylab = ylab, newdata = newdata)
show_marker_one(data.frame(W1 = 1), ylab = "Marker 1 (W1 = 1)")
show_marker_one(data.frame(W1 = -1), ylab = "Marker 1 (W1 = -1)")
rng_basis_two <- stacked_term(
poly_term(degree = 1, raw = TRUE, intercept = TRUE),
poly_term(degree = 1, raw = TRUE, intercept = TRUE) |>
weighted_term(W2))
show_marker_two <- \(newdata, ylab)
plot_marker(
time_fixef = stacked_term(
poly_term(degree = 2, raw = TRUE),
poly_term(degree = 2, raw = TRUE)|>
weighted_term(W2)),
time_rng = rng_basis_two,
fixef_vary = fixef_vary_marker[[2]], x_range = c(0, 10),
vcov_vary = vcov_vary[3:6, 3:6], ylab = ylab, newdata = newdata)
show_marker_two(data.frame(W2 = 1), ylab = "Marker 2 (W2 = 1)")
show_marker_two(data.frame(W2 = -1), ylab = "Marker 2 (W2 = -1)")
# the survival parameters
vcov_surv <- matrix(c(.2^2, .15^2, .15^2, .25^2), 2) # Xi
fixef_surv <- list(c(-1.25, .25), .75)
associations <- list(c(.4, -.3), c(-.5, .2))
fixef_vary_surv <- list(c(.5, .1, -.2, .11, -.2),
c(-1, -.25))
b_funcs <- list(
function(x, data)
cbind(bs(x, knots = 5, Boundary.knots = c(0, 10)),
data$Z1 * x),
function(x, data) ns(x, knots = 5, Boundary.knots = c(0, 10)))
# plot the log hazard with the 25%, 50% and 75% quantiles
library(SimSurvNMarker)
plot_terminal <- \(newdata, ylab)
plot_surv(
time_fixef = stacked_term(
bs_term(knots = 5, Boundary.knots = c(0, 10)),
poly_term(degree = 1, raw = TRUE) |>
weighted_term(Z1)),
time_rng = list(rng_basis_one, rng_basis_two),
x_range = c(0, 10), fixef_vary = fixef_vary_surv[[1]],
vcov_vary = vcov_vary, frailty_var = vcov_surv[1, 1], ps = c(.25, .5, .75),
associations = associations[[1]], log_hazard_shift = fixef_surv[[1]][1],
ylab = ylab, newdata = newdata)
plot_terminal(data.frame(Z1 = 1, W1 = 1, W2 = 1),
"Terminal event (Z1 = 1, W1 = 1, W2 = 1)")
plot_terminal(data.frame(Z1 = 0, W1 = 1, W2 = 1),
"Terminal event (Z1 = 0, W1 = 1, W2 = 1)")
plot_terminal(data.frame(Z1 = 0, W1 = 0, W2 = 0),
"Terminal event (Z1 = 0, W1 = 0, W2 = 0)")
plot_recurrent <- \(newdata, ylab)
plot_surv(
time_fixef = ns_term(knots = 5, Boundary.knots = c(0, 10)),
time_rng = list(rng_basis_one, rng_basis_two),
x_range = c(0, 10), fixef_vary = fixef_vary_surv[[2]],
vcov_vary = vcov_vary, frailty_var = vcov_surv[2, 2], ps = c(.25, .5, .75),
associations = associations[[2]], log_hazard_shift = fixef_surv[[2]][1],
ylab = ylab, newdata = newdata)
plot_recurrent(
data.frame(W1 = 1, W2 = 1), "Observation process (W1 = 1, W2 = 1)")
plot_recurrent(
data.frame(W1 = 0, W2 = 0), "Observation process (W1 = 0, W2 = 0)")
library(mvtnorm)
# simulates from the model by sampling a given number of individuals
sim_dat <- function(n_ids){
# simulate the outcomes
gl_dat <- get_gl_rule(100L)
dat <- lapply(1:n_ids, function(id){
# sample the terminal event time and the censoring time
cens <- min(rexp(1, rate = 1/10), 10)
U <- drop(rmvnorm(1, sigma = vcov_vary))
frailties <- drop(rmvnorm(1, sigma = vcov_surv))
Z1 <- c(1, runif(1, -1, 1))
W1 <- runif(1, -1, 1)
W2 <- runif(1, -1, 1)
data_pass <- data.frame(Z1 = Z1[2], W1 = W1, W2 = W2)
log_haz_offset <- sum(Z1 * fixef_surv[[1]]) + frailties[1]
# assign the conditional survival function
expansion <- function(x, b_func)
cbind(b_func(x, data_pass),
m_funcs[[1]](x, data_pass) %*% U[1:2],
m_funcs[[2]](x, data_pass) %*% U[3:6])
surv_func <- function(ti, fixef_vary_surv, associations, b_func){
formals(expansion)$b_func <- b_func
eval_surv_base_fun(
ti = ti, omega = c(fixef_vary_surv, associations), b_func = expansion,
gl_dat = gl_dat, delta = log_haz_offset)
}
# sample the survival time
rng <- runif(1)
root_func <- function(x, rng)
rng - surv_func(x, fixef_vary_surv = fixef_vary_surv[[1]],
associations = associations[[1]], b_func = b_funcs[[1]])
if(root_func(cens, rng) < 0){
# the observation is censored
y_terminal <- cens
event <- 0
} else {
# find the event time
root <- uniroot(root_func, c(0, cens), tol = 1e-6, rng = rng)
y_terminal <- root$root
event <- 1
}
terminal_outcome <- cbind(y = y_terminal, event = event, Z1 = Z1[2],
id = id, W1 = W1, W2 = W2)
# clean up
rm(list = setdiff(
ls(),
c("y_terminal", "terminal_outcome", "expansion", "surv_func", "frailties",
"U", "id", "data_pass", "W1", "W2")))
# simulate the observation times
Z2 <- 1
log_haz_offset <- sum(Z2 * fixef_surv[[2]]) + frailties[2]
root_func <- function(x, left_trunc_surv, rng)
rng - surv_func(x, fixef_vary_surv = fixef_vary_surv[[2]],
associations = associations[[2]], b_func = b_funcs[[2]]) /
left_trunc_surv
max_sample <- 1000L
left_trunc_surv <- 1
Z2 <- matrix(rep(Z2, each = max_sample), max_sample)
event <- y <- lf_trunc <- rep(NA_real_, max_sample)
lf_trunc_i <- 0
for(i in 1:max_sample){
# sample a random uniform variable and invert the survival function
rng_i <- runif(1)
lf_trunc[i] <- lf_trunc_i
if(root_func(y_terminal, left_trunc_surv, rng_i) < 0){
# the observation is right-censored and we can exit
y[i] <- y_terminal
event[i] <- 0
break
}
# we need to invert the survival function to find the observation time
root <- uniroot(root_func, c(lf_trunc_i, y_terminal), tol = 1e-6,
left_trunc_surv = left_trunc_surv, rng = rng_i)
lf_trunc_i <- y[i] <- root$root
event[i] <- 1
left_trunc_surv <- surv_func(
y[i], fixef_vary_surv = fixef_vary_surv[[2]], associations = associations[[2]],
b_func = b_funcs[[2]])
}
colnames(Z2) <- paste0("Z", 1:NCOL(Z2) - 1L)
obs_process <- cbind(lf_trunc = lf_trunc[1:i], y = y[1:i],
event = event[1:i], Z2[1:i, -1, drop = FALSE],
id = id, W1 = W1, W2 = W2)
# clean up
rm(list = setdiff(
ls(),
c("terminal_outcome", "U", "id", "obs_process", "data_pass", "W1", "W2")))
# sample the fixed effect covariates
obs_time <- c(0, obs_process[obs_process[, "event"] == 1, "y"])
n_obs <- length(obs_time)
X1 <- cbind(1, rnorm(n_obs), W1)
X2 <- cbind(matrix(1, n_obs), W2)
colnames(X1) <- paste0("X1_", 1:NCOL(X1) - 1L)
colnames(X1)[NCOL(X1)] <- "W1"
colnames(X2) <- paste0("X2_", 1:NCOL(X2) - 1L)
colnames(X2)[NCOL(X2)] <- "W2"
X <- list(X1, X2)
# sample the outcomes
eta <- sapply(1:2, function(i)
X[[i]] %*% fixef_marker[[i]] +
drop(g_funcs[[i]](obs_time, data_pass)) %*% fixef_vary_marker[[i]] +
drop(m_funcs[[i]](obs_time, data_pass)) %*% U[if(i == 1) 1:2 else 3:6])
ys <- eta + rmvnorm(n_obs, sigma = vcov_marker)
colnames(ys) <- paste0("Y", 1:2)
# mask some observations
do_mask <- sample.int(3L, n_obs, replace = TRUE)
ys[do_mask == 2, 1] <- NA
ys[do_mask == 3, 2] <- NA
X <- do.call(cbind, lapply(X, function(x) x[, -1, drop = FALSE]))
marker_data <- cbind(ys, X, time = obs_time, id = id)
return(list(marker_data = marker_data, obs_process = obs_process,
terminal_outcome = terminal_outcome))
})
# combine the data and return
marker_data <- as.data.frame(do.call(
rbind, lapply(dat, `[[`, "marker_data")))
marker_data$id <- as.integer(marker_data$id)
# the order does not matter
marker_data <- marker_data[sample.int(NROW(marker_data)), ]
obs_process <- as.data.frame(do.call(
rbind, lapply(dat, `[[`, "obs_process")))
obs_process$id <- as.integer(obs_process$id)
# the order does not matter
obs_process <- obs_process[sample.int(NROW(obs_process)), ]
terminal_outcome <- as.data.frame(do.call(
rbind, lapply(dat, `[[`, "terminal_outcome")))
terminal_outcome$id <- as.integer(terminal_outcome$id)
# the order does not matter
terminal_outcome <- terminal_outcome[sample.int(NROW(terminal_outcome)), ]
list(marker_data = marker_data, obs_process = obs_process,
terminal_outcome = terminal_outcome)
}
# sample a moderately large data set
set.seed(2)
dat <- sim_dat(1000L)
# we show a few properties of the data below
mean(dat$terminal_outcome$event) # mean event rate
#> [1] 0.713
sum(dat$obs_process$event) # number of observed markers less the individuals
#> [1] 6919
NROW(dat$marker_data) # number of observed markers less the individuals
#> [1] 7919
# distribution of observed marker per individual
proportions(table(table(dat$obs_process$id)))
#>
#> 1 2 3 4 5 6 7 8 9 10 11 12 13
#> 0.214 0.161 0.121 0.102 0.069 0.062 0.039 0.040 0.024 0.019 0.018 0.018 0.010
#> 14 15 16 17 18 19 20 21 22 23 24 25 27
#> 0.007 0.008 0.005 0.002 0.008 0.004 0.005 0.002 0.004 0.004 0.002 0.004 0.002
#> 28 29 30 31 32 33 38 39 40 41 43 45 46
#> 0.005 0.005 0.001 0.003 0.001 0.002 0.001 0.001 0.003 0.001 0.002 0.001 0.001
#> 49 54 58 63 66 72 77 81 82 84 98 119 135
#> 0.001 0.001 0.001 0.001 0.001 0.002 0.001 0.001 0.001 0.001 0.001 0.001 0.001
#> 172 220 247 288 318
#> 0.001 0.001 0.001 0.001 0.001
# show data for one individual
subset(dat$marker_data, id == 1)
#> Y1 Y2 X1_1 W1 W2 time id
#> X.1 NA 3.894 -1.1152 -0.1104 -0.85 0.3592 1
#> X -1.462 2.890 -0.3663 -0.1104 -0.85 0.0000 1
#> X.4 2.027 6.237 0.3244 -0.1104 -0.85 1.4799 1
#> X.3 NA 5.071 1.1210 -0.1104 -0.85 1.1119 1
#> X.2 -1.951 4.066 -0.9756 -0.1104 -0.85 0.4256 1
#> X.5 NA 6.153 -0.5684 -0.1104 -0.85 1.7229 1
subset(dat$obs_process, id == 1)
#> lf_trunc y event id W1 W2
#> 1 0.0000 0.3592 1 1 -0.1104 -0.85
#> 6 1.7229 2.2592 0 1 -0.1104 -0.85
#> 2 0.3592 0.4256 1 1 -0.1104 -0.85
#> 3 0.4256 1.1119 1 1 -0.1104 -0.85
#> 4 1.1119 1.4799 1 1 -0.1104 -0.85
#> 5 1.4799 1.7229 1 1 -0.1104 -0.85
subset(dat$terminal_outcome, id == 1)
#> y event Z1 id W1 W2
#> 1 2.259 1 -0.5483 1 -0.1104 -0.85
# estimate the model with this package. Get the object we need for the
# optimization
marker_1 <- marker_term(
Y1 ~ X1_1 + W1, id = id, subset(dat$marker_data, !is.na(Y1)),
time_fixef =
stacked_term(
ns_term(time, knots = c(3.33, 6.67), Boundary.knots = c(0, 10)),
poly_term(time, degree = 1, raw = TRUE) |>
weighted_term(W1)),
time_rng =
stacked_term(
poly_term(time, degree = 0, intercept = TRUE, raw = TRUE),
poly_term(time, degree = 0, intercept = TRUE, raw = TRUE) |>
weighted_term(W1)))
marker_2 <- marker_term(
Y2 ~ W2, id = id, subset(dat$marker_data, !is.na(Y2)),
time_fixef =
stacked_term(
poly_term(time, degree = 2, raw = TRUE),
poly_term(time, degree = 2, raw = TRUE) |>
weighted_term(W2)),
time_rng =
stacked_term(
poly_term(time, degree = 1, raw = TRUE, intercept = TRUE),
poly_term(time, degree = 1, raw = TRUE, intercept = TRUE) |>
weighted_term(W2)))
library(survival)
surv_terminal <- surv_term(
Surv(y, event) ~ Z1, id = id, dat$terminal_outcome,
time_fixef =
stacked_term(
bs_term(y, knots = 5, Boundary.knots = c(0, 10)),
poly_term(y, degree = 1, raw = TRUE) |>
weighted_term(Z1)),
with_frailty = TRUE)
surv_obs <- surv_term(
Surv(lf_trunc, y, event) ~ 1, id = id, dat$obs_process,
time_fixef = ns_term(y, knots = 5, Boundary.knots = c(0, 10)),
with_frailty = TRUE)
comp_obj <- joint_ms_ptr(markers = list(marker_1, marker_2),
survival_terms = list(surv_terminal, surv_obs),
max_threads = 4L)
rm(marker_1, marker_2, surv_terminal, surv_obs)
# get the starting values
system.time(start_val <- joint_ms_start_val(comp_obj, gr_tol = .1))
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
#> unable to evaluate scaled gradient
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
#> Model failed to converge: degenerate Hessian with 1 negative eigenvalues
#> user system elapsed
#> 180.653 0.144 46.950
# lower bound at the starting values
print(-attr(start_val, "value"), digits = 8)
#> [1] -8865.4428
# check that the gradient is correct
f <- function(x){
start_val[seq_along(x)] <- x
joint_ms_lb(comp_obj, start_val)
}
all.equal(numDeriv::grad(f, head(start_val, 53 + 2 * 44)),
head(joint_ms_lb_gr(comp_obj, start_val), 53 + 2 * 44),
tolerance = 1e-6)
#> [1] TRUE
# find the maximum lower bound estimate
system.time(opt_out <- joint_ms_opt(
comp_obj, par = start_val, max_it = 10000L, pre_method = 3L,
cg_tol = .1, c2 = .1, gr_tol = .1))
#> user system elapsed
#> 771.851 0.408 195.882
# we set gr_tol in the call so this is the convergence criterion for the
# gradient
sqrt(sum(joint_ms_lb_gr(comp_obj, opt_out$par)^2))
#> [1] 0.08981
opt_out$info # convergence code (0 == 'OK')
#> [1] 0
print(-opt_out$value, digits = 8) # maximum lower bound value
#> [1] -7333.7229
opt_out$counts
#> function gradient n_cg
#> 2508 1835 17924
# compare the estimates with the actual values. Start with the fixed effects
fmt_ests <- joint_ms_format(comp_obj, opt_out$par)
# the parameters for the first marker
fmt_ests$markers[[1]]
#> $fixef
#>
#> -0.5297 1.9959 1.0058
#>
#> $fixef_vary
#>
#> 1.3503 1.3160 -2.0364 0.1038
fixef_marker[[1]] # true values
#> [1] -0.5 2.0 1.0
fixef_vary_marker[[1]] # true values
#> [1] 1.4 1.2 -2.1 0.1
# the parameters for the second marker
fmt_ests$markers[[2]]
#> $fixef
#>
#> 1.004 -1.030
#>
#> $fixef_vary
#>
#> 0.57345 -0.02072 0.30081 -0.10075
fixef_marker[[2]] # true values
#> [1] 1 -1
fixef_vary_marker[[2]] # true values
#> [1] 0.50 -0.02 0.25 -0.10
# the fixed effects for the survival outcome and the association parameters
# for the terminal event
fmt_ests$survival[[1]]
#> $fixef
#>
#> -1.2353 0.2632
#>
#> $fixef_vary
#>
#> 0.8067 -0.3803 -0.8106 0.1701 -0.1939
#>
#> $associations
#>
#> 0.2831 -0.2842
fixef_surv[[1]]
#> [1] -1.25 0.25
fixef_vary_surv[[1]]
#> [1] 0.50 0.10 -0.20 0.11 -0.20
associations[[1]]
#> [1] 0.4 -0.3
# same for the observation process
fmt_ests$survival[[2]]
#> $fixef
#>
#> 0.7519
#>
#> $fixef_vary
#>
#> -0.8182 -0.1391
#>
#> $associations
#>
#> -0.5121 0.1978
fixef_surv[[2]]
#> [1] 0.75
fixef_vary_surv[[2]]
#> [1] -1.00 -0.25
associations[[2]]
#> [1] -0.5 0.2
# the parameters for covariance matrix of the random effects
fmt_ests$vcov$vcov_vary
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 0.71220 -0.08031 0.20699 -0.105911 -0.09079 -0.026239
#> [2,] -0.08031 0.45789 0.05584 -0.059417 0.03967 -0.012826
#> [3,] 0.20699 0.05584 0.89696 -0.142425 -0.09388 0.068093
#> [4,] -0.10591 -0.05942 -0.14242 0.448486 0.02705 0.006592
#> [5,] -0.09079 0.03967 -0.09388 0.027049 0.17862 0.070847
#> [6,] -0.02624 -0.01283 0.06809 0.006592 0.07085 0.480699
vcov_vary # the true values
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 0.622 -0.020 0.130 -0.081 -0.103 0.016
#> [2,] -0.020 0.473 0.121 0.022 -0.032 0.059
#> [3,] 0.130 0.121 0.801 -0.009 -0.193 0.000
#> [4,] -0.081 0.022 -0.009 0.529 0.040 0.050
#> [5,] -0.103 -0.032 -0.193 0.040 0.636 -0.144
#> [6,] 0.016 0.059 0.000 0.050 -0.144 0.354
# the parameters for the error term covariance matrix
fmt_ests$vcov$vcov_marker
#> [,1] [,2]
#> [1,] 0.3805 0.1128
#> [2,] 0.1128 0.1673
vcov_marker
#> [,1] [,2]
#> [1,] 0.36 0.10
#> [2,] 0.10 0.16
# the parameters for the frailty covariance matrix
fmt_ests$vcov$vcov_surv
#> [,1] [,2]
#> [1,] 0.02250 0.03248
#> [2,] 0.03248 0.05625
vcov_surv
#> [,1] [,2]
#> [1,] 0.0400 0.0225
#> [2,] 0.0225 0.0625Two Markers, the Observation Time Process, a Terminal Event, and Mixed Dependencies
We simulate and fit a model in this section where we have two markers, a recurrent event which is the observation times, and a terminal event like in the previous section. Unlike in the previous section, we let the two time-to-event outcomes depend on the integral or derivative of the deviation of the mean marker and the current value.
# settings for the simulation
library(VAJointSurv)
g1 <- ns_term(knots = c(3.33, 6.67), Boundary.knots = c(0, 10))
g2 <- poly_term(degree = 2, raw = TRUE)
g_funcs <- list(g1$eval, g2$eval)
m1 <- ns_term(knots = numeric(), Boundary.knots = c(0, 10), intercept = TRUE)
m2 <- poly_term(degree = 1, intercept = TRUE, raw = TRUE)
m_funcs <- list(m1$eval, m2$eval)
fixef_vary_marker <- list(c(1.4, 1.2, -2.1), c(.5, -.02)) # beta
fixef_marker <- list(c(-.5, 2), 1) # gamma
# Psi
vcov_vary <- structure(c(0.35, 0.08, -0.05, 0.01, 0.08, 1.92, -0.24, -0.04,
-0.05, -0.24, 0.32, 0.09, 0.01, -0.04, 0.09, 0.12),
.Dim = c(4L, 4L))
vcov_marker <- matrix(c(.6^2, .1, .1, .4^2), 2)
# plot the markers' mean curve
par(mar = c(5, 5, 1, 1))
plot_marker(
time_fixef = g1, time_rng = m1,
fixef_vary = fixef_vary_marker[[1]], x_range = c(0, 10),
vcov_vary = vcov_vary[1:2, 1:2], ylab = "Marker 1")
plot_marker(
time_fixef = g2, time_rng = m2,
fixef_vary = fixef_vary_marker[[2]], x_range = c(0, 10),
vcov_vary = vcov_vary[3:4, 3:4], ylab = "Marker 2")
# the survival parameters
vcov_surv <- matrix(c(.2^2, .15^2, .15^2, .25^2), 2) # Xi
fixef_surv <- list(c(-1, .25), .2)
# needs more parameters now
associations <- list(c(-.4, .6, -.4), c(-.7, .2, -1))
fixef_vary_surv <- list(c(.5, .1, -.2, .11),
c(-1, -.25))
# specify the dependence with the random effect from the markers
ders <- list(
list(c(-1L, 0L), # cumulative and present value
0L), # present value
list(0L, # present value
c(0L, 1L))) # present value and derivative
b_funcs <- list(
function(x) bs(x, knots = 5, Boundary.knots = c(0, 10)),
function(x) ns(x, knots = 5, Boundary.knots = c(0, 10)))
# plot the log hazard with the 25%, 50% and 75% quantiles
library(SimSurvNMarker)
plot_surv(
time_fixef = bs_term(knots = 5, Boundary.knots = c(0, 10)),
time_rng = list(m1, m2), ders = ders[[1]],
x_range = c(0, 10), fixef_vary = fixef_vary_surv[[1]],
vcov_vary = vcov_vary, frailty_var = vcov_surv[1, 1], ps = c(.25, .5, .75),
associations = associations[[1]], log_hazard_shift = fixef_surv[[1]][1],
ylab = "Terminal event")
plot_surv(
time_fixef = ns_term(knots = 5, Boundary.knots = c(0, 10)),
time_rng = list(m1, m2), ders = ders[[2]],
x_range = c(0, 10), fixef_vary = fixef_vary_surv[[2]],
vcov_vary = vcov_vary, frailty_var = vcov_surv[2, 2], ps = c(.25, .5, .75),
associations = associations[[2]], log_hazard_shift = fixef_surv[[2]][1],
ylab = "Observation process")
library(mvtnorm)
# simulates from the model by sampling a given number of individuals
sim_dat <- function(n_ids){
# simulate the outcomes
gl_dat <- get_gl_rule(100L)
dat <- lapply(1:n_ids, function(id){
# sample the terminal event time and the censoring time
cens <- min(rexp(1, rate = 1/10), 10)
U <- drop(rmvnorm(1, sigma = vcov_vary))
frailties <- drop(rmvnorm(1, sigma = vcov_surv))
Z1 <- c(1, runif(1, -1, 1))
log_haz_offset <- sum(Z1 * fixef_surv[[1]]) + frailties[1]
# assign the conditional survival function
expansion <- function(x, b_func)
cbind(b_func(x),
t(U[1:2] %*% m_funcs[[1]](x, ders[[1]][[1]][1])),
t(U[1:2] %*% m_funcs[[1]](x, ders[[1]][[1]][2])),
t(U[3:4] %*% m_funcs[[2]](x, ders[[1]][[2]][1])))
surv_func <- function(ti, fixef_vary_surv, associations, b_func){
formals(expansion)$b_func <- b_func
eval_surv_base_fun(
ti = ti, omega = c(fixef_vary_surv, associations), b_func = expansion,
gl_dat = gl_dat, delta = log_haz_offset)
}
# sample the survival time
rng <- runif(1)
root_func <- function(x, rng)
rng - surv_func(x, fixef_vary_surv = fixef_vary_surv[[1]],
associations = associations[[1]], b_func = b_funcs[[1]])
if(root_func(cens, rng) < 0){
# the observation is censored
y_terminal <- cens
event <- 0
} else {
# find the event time
root <- uniroot(root_func, c(0, cens), tol = 1e-6, rng = rng)
y_terminal <- root$root
event <- 1
}
terminal_outcome <- cbind(y = y_terminal, event = event, Z1 = Z1[2],
id = id)
# clean up
rm(list = setdiff(ls(), c("y_terminal", "terminal_outcome", "expansion",
"surv_func", "frailties", "U", "id")))
# simulate the observation times
Z2 <- 1
log_haz_offset <- sum(Z2 * fixef_surv[[2]]) + frailties[2]
expansion <- function(x, b_func)
cbind(b_func(x),
t(U[1:2] %*% m_funcs[[1]](x, ders[[2]][[1]][1])),
t(U[3:4] %*% m_funcs[[2]](x, ders[[2]][[2]][1])),
t(U[3:4] %*% m_funcs[[2]](x, ders[[2]][[2]][2])))
root_func <- function(x, left_trunc_surv, rng)
rng - surv_func(x, fixef_vary_surv = fixef_vary_surv[[2]],
associations = associations[[2]], b_func = b_funcs[[2]]) /
left_trunc_surv
max_sample <- 1000L
left_trunc_surv <- 1
Z2 <- matrix(rep(Z2, each = max_sample), max_sample)
event <- y <- lf_trunc <- rep(NA_real_, max_sample)
lf_trunc_i <- 0
for(i in 1:max_sample){
# sample a random uniform variable and invert the survival function
rng_i <- runif(1)
lf_trunc[i] <- lf_trunc_i
if(root_func(y_terminal, left_trunc_surv, rng_i) < 0){
# the observation is right-censored and we can exit
y[i] <- y_terminal
event[i] <- 0
break
}
# we need to invert the survival function to find the observation time
root <- uniroot(root_func, c(lf_trunc_i, y_terminal), tol = 1e-6,
left_trunc_surv = left_trunc_surv, rng = rng_i)
lf_trunc_i <- y[i] <- root$root
event[i] <- 1
left_trunc_surv <- surv_func(
y[i], fixef_vary_surv = fixef_vary_surv[[2]], associations = associations[[2]],
b_func = b_funcs[[2]])
}
colnames(Z2) <- paste0("Z", 1:NCOL(Z2) - 1L)
obs_process <- cbind(lf_trunc = lf_trunc[1:i], y = y[1:i],
event = event[1:i], Z2[1:i, -1, drop = FALSE],
id = id)
# clean up
rm(list = setdiff(ls(), c("terminal_outcome", "U", "id",
"obs_process")))
# sample the fixed effect covariates
obs_time <- c(0, obs_process[obs_process[, "event"] == 1, "y"])
n_obs <- length(obs_time)
X1 <- cbind(1, rnorm(n_obs))
X2 <- matrix(1, n_obs)
colnames(X1) <- paste0("X1_", 1:NCOL(X1) - 1L)
colnames(X2) <- paste0("X2_", 1:NCOL(X2) - 1L)
X <- list(X1, X2)
# sample the outcomes
eta <- sapply(1:2, function(i)
X[[i]] %*% fixef_marker[[i]] +
drop(t(g_funcs[[i]](obs_time))) %*% fixef_vary_marker[[i]] +
drop(t(m_funcs[[i]](obs_time))) %*% U[1:2 + (i == 2) * 2])
ys <- eta + rmvnorm(n_obs, sigma = vcov_marker)
colnames(ys) <- paste0("Y", 1:2)
# mask some observations
do_mask <- sample.int(3L, n_obs, replace = TRUE)
ys[do_mask == 2, 1] <- NA
ys[do_mask == 3, 2] <- NA
X <- do.call(cbind, lapply(X, function(x) x[, -1, drop = FALSE]))
marker_data <- cbind(ys, X, time = obs_time, id = id)
return(list(marker_data = marker_data, obs_process = obs_process,
terminal_outcome = terminal_outcome))
})
# combine the data and return
marker_data <- as.data.frame(do.call(
rbind, lapply(dat, `[[`, "marker_data")))
marker_data$id <- as.integer(marker_data$id)
# the order does not matter
marker_data <- marker_data[sample.int(NROW(marker_data)), ]
obs_process <- as.data.frame(do.call(
rbind, lapply(dat, `[[`, "obs_process")))
obs_process$id <- as.integer(obs_process$id)
# the order does not matter
obs_process <- obs_process[sample.int(NROW(obs_process)), ]
terminal_outcome <- as.data.frame(do.call(
rbind, lapply(dat, `[[`, "terminal_outcome")))
terminal_outcome$id <- as.integer(terminal_outcome$id)
# the order does not matter
terminal_outcome <- terminal_outcome[sample.int(NROW(terminal_outcome)), ]
list(marker_data = marker_data, obs_process = obs_process,
terminal_outcome = terminal_outcome)
}
# sample a moderately large data set
set.seed(2)
dat <- sim_dat(1000L)
# we show a few properties of the data below
mean(dat$terminal_outcome$event) # mean event rate
#> [1] 0.77
sum(dat$obs_process$event) # number of observed markers less the individuals
#> [1] 2425
NROW(dat$marker_data) # number of observed markers less the individuals
#> [1] 3425
# distribution of observed marker per individual
proportions(table(table(dat$obs_process$id)))
#>
#> 1 2 3 4 5 6 7 8 9 10 11 12 13
#> 0.319 0.190 0.141 0.109 0.078 0.041 0.033 0.027 0.013 0.012 0.015 0.003 0.004
#> 14 16 17 27 29 31 38
#> 0.005 0.003 0.003 0.001 0.001 0.001 0.001
# show data for one individual
subset(dat$marker_data, id == 1)
#> Y1 Y2 X1_1 time id
#> X.6 NA 8.382 -0.56843 7.8724 1
#> X.3 0.7605 NA -0.97556 3.6495 1
#> X.7 -3.5564 9.325 -1.03442 8.4949 1
#> X.2 -0.5523 NA -1.11525 0.4897 1
#> X.9 -3.9194 9.798 -0.98139 9.5827 1
#> X.1 NA 1.726 -0.36632 0.4271 1
#> X NA 1.987 -0.02815 0.0000 1
#> X.4 4.1065 6.374 1.12102 4.9005 1
#> X.8 -1.5241 9.405 -0.53302 8.6015 1
#> X.5 0.3933 NA 0.32444 7.5675 1
subset(dat$obs_process, id == 1)
#> lf_trunc y event id
#> 8 8.4949 8.6015 1 1
#> 9 8.6015 9.5827 1 1
#> 2 0.4271 0.4897 1 1
#> 1 0.0000 0.4271 1 1
#> 7 7.8724 8.4949 1 1
#> 10 9.5827 10.0000 0 1
#> 5 4.9005 7.5675 1 1
#> 6 7.5675 7.8724 1 1
#> 4 3.6495 4.9005 1 1
#> 3 0.4897 3.6495 1 1
subset(dat$terminal_outcome, id == 1)
#> y event Z1 id
#> 1 10 0 -0.6384 1
# estimate the model with this package. Get the object we need for the
# optimization
marker_1 <- marker_term(
Y1 ~ X1_1, id = id, subset(dat$marker_data, !is.na(Y1)),
time_fixef = ns_term(time, knots = c(3.33, 6.67), Boundary.knots = c(0, 10)),
time_rng = ns_term(time, knots = numeric(), Boundary.knots = c(0, 10),
intercept = TRUE))
marker_2 <- marker_term(
Y2 ~ 1, id = id, subset(dat$marker_data, !is.na(Y2)),
time_fixef = poly_term(time, degree = 2, raw = TRUE),
time_rng = poly_term(time, degree = 1, raw = TRUE, intercept = TRUE))
library(survival)
surv_terminal <- surv_term(
Surv(y, event) ~ Z1, id = id, dat$terminal_outcome,
time_fixef = bs_term(y, knots = 5, Boundary.knots = c(0, 10)),
with_frailty = TRUE)
surv_obs <- surv_term(
Surv(lf_trunc, y, event) ~ 1, id = id, dat$obs_process,
time_fixef = ns_term(y, knots = 5, Boundary.knots = c(0, 10)),
with_frailty = TRUE)
comp_obj <- joint_ms_ptr(markers = list(marker_1, marker_2),
survival_terms = list(surv_terminal, surv_obs),
max_threads = 4L, ders = ders)
rm(marker_1, marker_2, surv_terminal, surv_obs)
# get the starting values
system.time(start_val <- joint_ms_start_val(comp_obj, gr_tol = .1))
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
#> Model failed to converge with max|grad| = 0.00249371 (tol = 0.002, component 1)
#> user system elapsed
#> 25.047 0.012 6.636
# lower bound at the starting values
print(-attr(start_val, "value"), digits = 8)
#> [1] -7950.3559
# check that the gradient is correct
f <- function(x){
start_val[seq_along(x)] <- x
joint_ms_lb(comp_obj, start_val)
}
all.equal(numDeriv::grad(f, head(start_val, 39 + 2 * 27)),
head(joint_ms_lb_gr(comp_obj, start_val), 39 + 2 * 27),
tolerance = 1e-6)
#> [1] TRUE
# find the maximum lower bound estimate
system.time(opt_out <- joint_ms_opt(comp_obj, par = start_val, max_it = 1000L,
pre_method = 3L, cg_tol = .2, c2 = .1,
gr_tol = .1))
#> user system elapsed
#> 242.687 0.152 60.773
# we set gr_tol in the call so this is the convergence criterion for the
# gradient
sqrt(sum(joint_ms_lb_gr(comp_obj, opt_out$par)^2))
#> [1] 0.09699
opt_out$info # convergence code (0 == 'OK')
#> [1] 0
print(-opt_out$value, digits = 8) # maximum lower bound value
#> [1] -7596.2613
opt_out$counts
#> function gradient n_cg
#> 2851 2059 14027
# compare the estimates with the actual values. Start with the fixed effects
fmt_ests <- joint_ms_format(comp_obj, opt_out$par)
# the parameters for the first marker
fmt_ests$markers[[1]]
#> $fixef
#>
#> -0.4302 2.0078
#>
#> $fixef_vary
#>
#> 1.355 1.218 -1.801
fixef_marker[[1]] # true values
#> [1] -0.5 2.0
fixef_vary_marker[[1]] # true values
#> [1] 1.4 1.2 -2.1
# the parameters for the second marker
fmt_ests$markers[[2]]
#> $fixef
#>
#> 0.9901
#>
#> $fixef_vary
#>
#> 0.53124 -0.02203
fixef_marker[[2]] # true values
#> [1] 1
fixef_vary_marker[[2]] # true values
#> [1] 0.50 -0.02
# the fixed effects for the survival outcome and the association parameters
# for the terminal event
fmt_ests$survival[[1]]
#> $fixef
#>
#> -1.0857 0.2989
#>
#> $fixef_vary
#>
#> 0.74032 0.08106 -0.10742 -0.42202
#>
#> $associations
#>
#> -0.4093 0.5460 -0.4314
fixef_surv[[1]]
#> [1] -1.00 0.25
fixef_vary_surv[[1]]
#> [1] 0.50 0.10 -0.20 0.11
associations[[1]]
#> [1] -0.4 0.6 -0.4
# same for the observation process
fmt_ests$survival[[2]]
#> $fixef
#>
#> 0.2204
#>
#> $fixef_vary
#>
#> -1.2801 -0.2923
#>
#> $associations
#>
#> -0.7104 0.1854 -0.7896
fixef_surv[[2]]
#> [1] 0.2
fixef_vary_surv[[2]]
#> [1] -1.00 -0.25
associations[[2]]
#> [1] -0.7 0.2 -1.0
# the parameters for covariance matrix of the random effects
fmt_ests$vcov$vcov_vary
#> [,1] [,2] [,3] [,4]
#> [1,] 0.34379 0.2413 -0.07381 -0.03167
#> [2,] 0.24130 2.1631 -0.26000 -0.11238
#> [3,] -0.07381 -0.2600 0.27222 0.10383
#> [4,] -0.03167 -0.1124 0.10383 0.12941
vcov_vary # the true values
#> [,1] [,2] [,3] [,4]
#> [1,] 0.35 0.08 -0.05 0.01
#> [2,] 0.08 1.92 -0.24 -0.04
#> [3,] -0.05 -0.24 0.32 0.09
#> [4,] 0.01 -0.04 0.09 0.12
# the parameters for the error term covariance matrix
fmt_ests$vcov$vcov_marker
#> [,1] [,2]
#> [1,] 0.3332 0.0930
#> [2,] 0.0930 0.1621
vcov_marker
#> [,1] [,2]
#> [1,] 0.36 0.10
#> [2,] 0.10 0.16
# the parameters for the frailty covariance matrix
fmt_ests$vcov$vcov_surv
#> [,1] [,2]
#> [1,] 0.06259 -0.04691
#> [2,] -0.04691 0.03692
vcov_surv
#> [,1] [,2]
#> [1,] 0.0400 0.0225
#> [2,] 0.0225 0.0625Technical Details
We provide a few technical details in this section. The concatenated coefficient is given by
\[\vec\theta = \begin{pmatrix} \vec\gamma \\ \vec\beta\\ \vec\omega_1 \\ \vec\delta_1 \\ \vec\alpha_1 \\ \vdots \\ \vec\omega_H \\ \vec\delta_H \\ \vec\alpha_H \\ \text{vec}(\Sigma) \\ \text{vec}(\Psi) \\ \text{vec}(\Xi) \end{pmatrix}\]
where \(\text{vec}(\cdot)\) is the vectorization function that stacks the column of a matrix on top of each other. \(\vec\theta\) is with the parameters from the variational approximation. We work with log Cholesky decompositions of the covariance matrices in practice. The joint_ms_format is a utility functions which returns a list with each of parameter separately.
There are different types of terms in the lower bound in the GVA. We cover the types below. The lower bound has many parameters even with a moderate amount of parameters. However, we quickly optimize the lower bound using the psqn package.
Kullback–Leibler Divergence Term
In the GVA, we assume that the conditional distribution of \((\vec U_i^\top, \vec\xi_i^\top)^\top\) is a normal distribution with mean \(\vec\zeta_i\) and covariance matrix \(\Omega_i\). One of the terms in the lower bound of the GVA is the Kullback–Leibler (KL) divergence term between the unconditional distribution of \((\vec U_i^\top, \vec\xi_i^\top)^\top\) and the assumed conditional distribution. This term is given by
\[\begin{multline*} \frac{1}{2}\Big(\log\lvert\Omega_i\rvert - \log\lvert\Psi\rvert - \log\lvert\Xi\rvert -\vec\zeta_{i,1:R}^\top\Psi^{-1}\vec\zeta_{i,1:R}-\vec\zeta_{i,(-1:R)}^\top\Xi^{-1}\vec\zeta_{i,(-1:R)} \\- \text{tr}\Omega_{i,1:R,1:R}\Psi^{-1}- \text{tr}\Omega_{i,(-1:R),(-1:R)}\Xi^{-1}+ R + H\Big)\end{multline*}\]
where \(\text{tr}(\cdot)\) returns the trace of a matrix. Thus, the derivatives w.r.t. \(\Omega_i\), \(\Psi\), \(\Xi\), \(\vec\zeta_i\) are respectively
\[\frac{1}{2}\left(\Omega_i^{-1} - \begin{pmatrix} \Psi^{-1} & 0 \\ 0 & \Xi^{-1}\end{pmatrix}\right)\]
\[\frac{1}{2}\Psi^{-1}(\vec\zeta_{i,1:R}\vec\zeta_{i,1:R}^\top + \Omega_{i,1:R,1:R} - \Psi)\Psi^{-1}\]
\[\frac{1}{2}\Xi^{-1}(\vec\zeta_{i,(-1:R)}\vec\zeta_{i,(-1:R)}^\top+ \Omega_{i,(-1:R),(-1:R)} - \Xi)\Xi^{-1}.\]
\[-\begin{pmatrix} \Psi^{-1} & 0 \\ 0 & \Xi^{-1} \end{pmatrix}\vec\zeta_i\]
Marker Terms
The log conditional density term of observation \(j\) of individual \(i\) at time \(s_{ij}\) is
\[\int \log\left(\phi^{(L)}\left(\vec y_{ij}; X_i\vec\gamma + G(s_{ij})\vec\beta + M(s_{ij})\vec w, \Sigma\right)\right)\phi^{(R)}\left(\vec w; \vec\zeta_{i,1:R}, \Omega_{i,1:R,1:R}\right) d\vec w\] which gives
\[\begin{multline*}-\frac{L}{2}\log 2\pi -\frac {1}{2} \log\lvert\Sigma \rvert - \frac{1}{2}\left(\vec y_{ij} - X_i\vec\gamma - G(s_{ij})\vec\beta - M(s_{ij})\vec\zeta_{i,1:R}\right)^\top \\ \Sigma^{-1}\left(\vec y_{ij} - X_i\vec\gamma - G(s_{ij})\vec\beta - M(s_{ij})\vec\zeta_{i,1:R}\right) - \frac{1}{2}\text{tr}\Sigma^{-1}M(s_{ij})\Omega_{i,1:R,1:R}M(s_{ij})^\top. \end{multline*}\]
Survival Outcomes
Let \(T_{i1}, \dots, T_{iL}\) be the minimum of the censoring time or the observed event time and \(D_{i1}, \dots, D_{iL}\) be event indicators. The the lower bound terms for the \(k\)th time to event is
\[\int \left(d_{ik}\log h_{ik}(t_{ik}\mid \vec w_{1:R}, w_{R + k})-\int_0^{t_{ik}} h_{ik}(s\mid \vec w_{1:R}, w_{R + k})ds \right)\phi^{(R + H)}(\vec w; \vec\zeta_i, \Omega_i)d\vec w.\]
Left-truncation is easily handled by replacing the zero in the lower limit by the left-truncation point. The log hazard term gives the following lower bound terms
\[d_{ik}\left(\vec z_{ik}^\top\vec\delta_k + \omega_k^\top \vec b_k(t_{ik}) + \vec\alpha_k^\top M(t_{ik})\vec\zeta_{i,1:R} + \zeta_{i,R + k}\right).\]
As for the expected cumulative hazard term, assume that we can interchange the order of integration. Then the lower bound term is
\[-\exp(\vec z_{ik}^\top\vec\delta_k + \zeta_{R +k})\int_0^{t_{ij}} \exp\left(\omega_k^\top \vec b_k(s) + \vec\alpha_k^\top M(s)\vec\zeta_{i,1:R} +\frac{1}{2} (\vec\alpha_k^\top, 1) O_{ik}(s)\begin{pmatrix} \vec\alpha_k \\ 1 \end{pmatrix}\right) ds\]
where
\[O_{ik}(s) = \begin{pmatrix}M(s) & 0 \\ 0 & 1\end{pmatrix} \Omega_{(1:R, R + k), (1:R, R + k)} \begin{pmatrix}M(s)^\top & 0 \\ 0 & 1\end{pmatrix}.\]
References
Berg, Gerard J. van den, and Bettina Drepper. 2016. “Inference for Shared-Frailty Survival Models with Left-Truncated Data.” Econometric Reviews 35 (6): 1075–98. https://doi.org/10.1080/07474938.2014.975640.
Crowther, Michael J., Therese M.-L Andersson, Paul C. Lambert, Keith R. Abrams, and Keith Humphreys. 2016. “Joint Modelling of Longitudinal and Survival Data: Incorporating Delayed Entry and an Assessment of Model Misspecification.” Statistics in Medicine 35 (7): 1193–1209. https://doi.org/https://doi.org/10.1002/sim.6779.