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")