trivPenalNL {frailtypack}R Documentation

Fit a Non-Linear Trivariate Joint Model for Recurrent Events and a Terminal Event with a Biomarker Described with an ODE Population Model

Description

Fit a non-linear trivariate joint model for a longitudinal biomarker, recurrent events and a terminal event using a semiparametric penalized likelihood estimation or a parametric estimation on the hazard functions.

The values yi(t) (i=1,...,N) for N subjects represent the individual evolution of the biomarker e.g. tumor size expressed as the sum of the longest diameters (SLD) of target lesions. The dynamics of the biomarker are described by an ordinary differential equation (ODE) that includes the effect of the natural net growth and the treatment effect:

trivNLmodel1.png

The model includes the following parameters (using the interpretation of tumor dynamics): exp(KG,0) the constant tumor growth rate, exp(KD,0) the drug-induced tumor decline rate, \lambda resistance effect to drug (exponential tumor decay change with time), exp(y0) the initial level of the biomarker and di is the treatment concentration (e.g dose). The random effects biT = (by0,i,bG,i,bD,i,b\lambda,i)T are gaussian variables with a diagonal covariance matrix Bi. In the trivariate model we use the analytical solution of the equation with the population-based approach of the non-linear mixed effects model. We can also assume a transformation for the observations of the biomarker (one parameter Box-Cox transformation) and we include a gaussian measurement error, for individual i and observation k (k=1,...,ni), \epsilonik ~ N(0,\sigma\epsilon2).

The risks of the recurrent (rij(.) the risk of the jth event of the individual i) and terminal events (\lambdai the risk of the event of the individual i) are represented by proportional hazard risk models. The joint model is constructed assuming that the processes are linked via a latent structure and includes the non-linear mixed effects model for the longitudinal data:

trivNLmodel2.png

where XG,i(t), XD,i(t), XR,ij(t) and XT,i(t) are vectors of possible time-varying fixed effects covariates and \betaG, \betaD, \betaR and \betaT are the associated coefficients. The random effects bi are independent from the measurement error. The relationship between the biomarker and recurrent events is explained via g(yi(t)) with coefficients \etaR and between the biomarker and terminal event is explained via h(yi(t)) with coefficients \etaT. Currently, only one form of the functions g(.) and h(.) is available: the random effects bi. The frailty term vi is gaussian with mean 0 and variance \sigmav. Together with bi constitutes the random effects of the model:

trivNLmodel3.png

Any combination of the random effects bi, e.g. bi=by0,i or bi = {bG,i,bD,i,b\lambda,i} can be chosen for the model.

We consider that the longitudinal outcome can be a subject to a quantification limit, i.e. some observations, below a level of detection s cannot be quantified (left-censoring).

Usage

trivPenalNL(formula, formula.terminalEvent, biomarker, formula.KG,
formula.KD, dose, time.biomarker, data, data.Longi, random, id, link =
"Random-effects", BoxCox = FALSE, left.censoring = FALSE, recurrentAG =
FALSE, n.knots, kappa, maxit = 300, hazard = "Splines", init.B, init.Random,
init.Eta, init.Alpha, init.Biomarker, method.GH = "Standard", init.GH =
FALSE, n.nodes, LIMparam = 1e-3, LIMlogl = 1e-3, LIMderiv = 1e-3,
print.times = TRUE)

Arguments

formula

a formula object, with the response on the left of a \sim operator, and the terms on the right. The response must be a survival object as returned by the 'Surv' function like in survival package. Interactions are possible using * or :.

formula.terminalEvent

A formula object, only requires terms on the right to indicate which variables are modelling the terminal event. Interactions are possible using * or :.

biomarker

Name of the variable representing the longitudinal biomarker.

formula.KG

A formula object, only requires terms on the right to indicate which covariates related to the biomarker growth are included in the longitudinal sub-model. It must follow the standard form used for linear mixed-effects models. Interactions are possible using * or :.

formula.KD

A formula object, only requires terms on the right to indicate which covariates related to the biomarker drug-induced decline are included in the longitudinal sub-model. It must follow the standard form used for linear mixed-effects models. Interactions are possible using * or :.

dose

Name of the variable representing the drug concentration indicator.

time.biomarker

Name of the variable of times of biomarker measurements.

data

A 'data.frame' with the variables used in formula.

data.Longi

A 'data.frame' with the variables used in formula.KG, formula.KD, biomarker, dose, time.biomarker and id.

random

Names of parameters for which the random effects are included in the mixed model. The names must be chosen among "y0", "KG", "KD" and "lambda". Any combination of the mentioned names is allowed.

id

Name of the variable representing the individuals.

link

Type of link functions for the dependence between the biomarker and death and between the biomarker and the recurrent events: only "Random-effects" for the association directly via the random effects of the biomarker is allowed for the moment (option for a future extension).

BoxCox

Should the Box-Cox transformation be used for the longitudinal biomarker? If there is no transformation, the argument must be equal to FALSE, otherwise the of the transformation parameter must be given, then the transformed values are y^*=(y^{\xi}-1)/\xi, where \xi is the Box-Cox parameter.

left.censoring

Is the biomarker left-censored below a threshold s? If there is no left-censoring, the argument must be equal to FALSE, otherwise the value of the threshold must be given.

recurrentAG

Logical value. Is Andersen-Gill model fitted? If so indicates that recurrent event times with the counting process approach of Andersen and Gill is used. This formulation can be used for dealing with time-dependent covariates. The default is FALSE.

n.knots

Integer giving the number of knots to use. Value required in the penalized likelihood estimation. It corresponds to the (n.knots+2) splines functions for the approximation of the hazard or the survival functions. We estimate I or M-splines of order 4. When the user set a number of knots equals to k (n.knots=k) then the number of interior knots is (k-2) and the number of splines is (k-2)+order. Number of knots must be between 4 and 20. (See Note in frailtyPenal function)

kappa

Positive smoothing parameters in the penalized likelihood estimation. The coefficient kappa of the integral of the squared second derivative of hazard function in the fit (penalized log likelihood). To obtain an initial value for kappa, a solution is to fit the corresponding Cox model using cross validation (See cross.validation in function frailtyPenal). We advise the user to identify several possible tuning parameters, note their defaults and look at the sensitivity of the results to varying them.

maxit

Maximum number of iterations for the Marquardt algorithm. Default is 300

hazard

Type of hazard functions: "Splines" for semiparametric hazard functions using equidistant intervals or "Splines-per" using percentile with the penalized likelihood estimation, "Weibull" for the parametric Weibull functions. The default is "Splines".

init.B

Vector of initial values for regression coefficients. This vector should be of the same size as the whole vector of covariates with the first elements for the covariates related to the recurrent events, then to the terminal event and then to the biomarker (interactions in the end of each component). Default is 0.5 for each.

init.Random

Initial value for variance of the elements of the matrix of the distribution of the random effects.

init.Eta

Initial values for regression coefficients for the link functions, first for the recurrent events (\bold{\eta}_R) and for the terminal event (\bold{\eta}_T).

init.Alpha

Initial value for parameter alpha

init.Biomarker

Initial values for biomarker parameters: y_0, K_{G,0}, K_{D,0} and \lambda (using this order).

method.GH

Method for the Gauss-Hermite quadrature: "Standard" for the standard non-adaptive Gaussian quadrature and "Pseudo-adaptive" for the pseudo-adaptive Gaussian quadrature (see Details). The default is "Standard". When the option "Pseudo-adaptive" is chosen, then a univariate model (non-linear mixed model for the biomarker) is fitted in order to obtain the estimations of the random effects \bold{b}_i.

init.GH

Only when the opiton "Pseudo-adaptive" of the argument method.GH is chosen. If TRUE, the estimations of the biomarker parameters (y_0, K_{G,0}, K_{D,0} and \lambda), \sigma_{\epsilon}, \bold{\beta}_G and \bold{\beta}_D from the univariate mixed model are used as the initial values of the parameters related to the biomarker.

n.nodes

Number of nodes for the Gauss-Hermite quadrature (from 5 to 32). The default is 9.

LIMparam

Convergence threshold of the Marquardt algorithm for the parameters (see Details), 10^{-3} by default.

LIMlogl

Convergence threshold of the Marquardt algorithm for the log-likelihood (see Details), 10^{-3} by default.

LIMderiv

Convergence threshold of the Marquardt algorithm for the gradient (see Details), 10^{-3} by default.

print.times

a logical parameter to print iteration process. Default is TRUE.

Details

Typical usage for the joint model

trivPenalNL(Surv(time,event)~cluster(id) + var1 + var2 +
terminal(death), formula.terminalEvent =~ var1 + var3, biomarker =
"biomarker.name", dose = "dose.name", time.biomarker = "time", formula.KG ~
var1, formula.KD ~ var2, data, data.Longi, ...)

The method of the Gauss-Hermite quadrature for approximations of the multidimensional integrals, i.e. length of random more than 2, can be chosen among the standard (non-adaptive) and pseudo-adaptive in which the quadrature points are transformed using the information from the fitted mixed-effects model for the biomarker (Rizopoulos 2012) or multivariate non-adaptive procedure proposed by Genz et al. 1996 and implemented in FORTRAN subroutine HRMSYM. The choice of the method is important for estimations. The standard non-adaptive Gauss-Hermite quadrature ("Standard") with a specific number of points gives accurate results but can be time consuming. The pseudo-adaptive quadrature uses transformed quadrature points to center and scale the integrand by utilizing estimates of the random effects from an appropriate non-linear mixed-effects model (this transformation does not include the frailty in the trivariate model, for which the standard method, with 20 quadrature points, is used). This method enables using less quadrature points while preserving the estimation accuracy and thus lead to a better computational time.

NOTE. Data frames data and data.Longi must be consistent. Names and types of corresponding covariates must be the same, as well as the number and identification of individuals.

Value

The following components are included in a 'trivPenalNL' object for each model:

b

The sequence of the corresponding estimation of the coefficients for the hazard functions (parametric or semiparametric), the random effects variances and the regression coefficients.

call

The code used for the model.

formula

The formula part of the code used for the recurrent event part of the model.

formula.terminalEvent

The formula part of the code used for the terminal event part of the model.

formula.KG

The formula part of the code used for the longitudinal part of the model, for the biomarker growth dynamics.

formula.KD

The formula part of the code used for the longitudinal part of the model, for the biomarker decline dynamics.

coef

The regression coefficients (first for the recurrent events, then for the terminal event, then for the biomarker growth and then for the biomarker decline.

groups

The number of groups used in the fit.

kappa

The values of the smoothing parameters in the penalized likelihood estimation corresponding to the baseline hazard functions for the recurrent and terminal events.

logLikPenal

The complete marginal penalized log-likelihood in the semiparametric case.

logLik

The marginal log-likelihood in the parametric case.

n.measurements

The number of biomarker observations used in the fit.

max_rep

The maximal number of repeated measurements per individual.

n

The number of observations in 'data' (recurrent and terminal events) used in the fit.

n.events

The number of recurrent events observed in the fit.

n.deaths

The number of terminal events observed in the fit.

n.iter

The number of iterations needed to converge.

n.knots

The number of knots for estimating the baseline hazard function in the penalized likelihood estimation.

n.strat

The number of stratum.

varH

The variance matrix of all parameters (before positivity constraint transformation for the variance of the measurement error, for which the delta method is used).

varHIH

The robust estimation of the variance matrix of all parameters.

xR

The vector of times where both survival and hazard function of the recurrent events are estimated. By default seq(0,max(time),length=99), where time is the vector of survival times.

lamR

The array (dim=3) of baseline hazard estimates and confidence bands (recurrent events).

survR

The array (dim=3) of baseline survival estimates and confidence bands (recurrent events).

xD

The vector of times where both survival and hazard function of the terminal event are estimated. By default seq(0,max(time),length=99), where time is the vector of survival times.

lamD

The array (dim=3) of baseline hazard estimates and confidence bands.

survD

The array (dim=3) of baseline survival estimates and confidence bands.

medianR

The value of the median survival and its confidence bands for the recurrent event.

medianD

The value of the median survival and its confidence bands for the terminal event.

typeof

The type of the baseline hazard function (0:"Splines", "2:Weibull").

npar

The number of parameters.

nvar

The vector of number of explanatory variables for the recurrent events, terminal event, biomarker growth and biomarker decline.

nvarRec

The number of explanatory variables for the recurrent events.

nvarEnd

The number of explanatory variables for the terminal event.

nvarKG

The number of explanatory variables for the biomarker growth.

nvarKD

The number of explanatory variables for the biomarker decline.

noVarRec

The indicator of absence of the explanatory variables for the recurrent events.

noVarEnd

The indicator of absence of the explanatory variables for the terminal event.

noVarKG

The indicator of absence of the explanatory variables for the biomarker growth.

noVarKD

The indicator of absence of the explanatory variables for the biomarker decline.

LCV

The approximated likelihood cross-validation criterion in the semiparametric case (with H minus the converged Hessian matrix, and l(.) the full log-likelihood).

LCV=\frac{1}{n}(trace(H^{-1}_{pl}H) - l(.))

AIC

The Akaike information Criterion for the parametric case.

AIC=\frac{1}{n}(np - l(.))

n.knots.temp

The initial value for the number of knots.

shape.weib

The shape parameter for the Weibull hazard functions (the first element for the recurrences and the second one for the terminal event).

scale.weib

The scale parameter for the Weibull hazard functions (the first element for the recurrences and the second one for the terminal event).

random.effects.pred

The empirical Bayes predictions of the random effects (ie. using conditional posterior distributions).

global_chisq.testR

The binary variable equals to 0 when no multivariate Wald is given, 1 otherwise (for the recurrent part).

global_chisq.testT

The binary variable equals to 0 when no multivariate Wald is given, 1 otherwise (for the terminal part).

global_chisq.testKG

The binary variable equals to 0 when no multivariate Wald is given, 1 otherwise (for the biomarker growth).

global_chisq.testKD

The binary variable equals to 0 when no multivariate Wald is given, 1 otherwise (for the biomarker decline).

AG

The logical value. Is Andersen-Gill model fitted?

B1

The variance matrix of the random effects for the longitudinal outcome.

sigma2

The variance of the frailty term (\sigma_v).

alpha

The coefficient \alpha associated with the frailty parameter in the terminal hazard function.

ResidualSE

The variance of the measurement error.

etaR

The regression coefficients for the link function g(\cdot).

etaT

The regression coefficients for the link function h(\cdot).

ne_re

The number of random effects b used in the fit.

names.re

The names of variables for the random effects \bold{b}_i.

link

The name of the type of the link functions.

leftCensoring

The logical value. Is the longitudinal outcome left-censored?

leftCensoring.threshold

For the left-censored biomarker, the value of the left-censoring threshold used for the fit.

prop.censored

The fraction of observations subjected to the left-censoring.

methodGH

The Gaussian quadrature method used in the fit.

n.nodes

The number of nodes used for the Gaussian quadrature in the fit.

K_G0

Value of the estimate of the biomarker growth parameter.

K_D0

Value of the estimate of the biomarker decay parameter.

lambda

Value of the estimate of the biomarker resistance to drug.

y_0

Value of the estimate of the biomarker intial level.

biomarker

Name of the variable associated with the biomarker in the data.

time.biomarker

Name of the variable associated with the time of measurements of the biomarker in the data.

dose

Name of the variable associated with the drug concentration in the data.

BoxCox

The logical value. Is the BoxCox transformation applied for the biomarker?

BoxCox_parameter

The value of the BoxCox transformation parameter.

alpha_p.value

p-value of the Wald test for the estimated coefficient \alpha.

sigma2_p.value

p-value of the Wald test for the estimated variance of the frailty term (\sigma_v).

etaR_p.value

p-values of the Wald test for the estimated regression coefficients for the link function g(\cdot).

etaT_p.value

p-values of the Wald test for the estimated regression coefficients for the link function h(\cdot).

y_0_p.value

p-value of the Wald test for the estimated biomarker intial level.

K_G0_p.value

p-value of the Wald test for the estimated biomarker growth parameter.

K_D0_p.value

p-value of the Wald test for the estimated biomarker decay parameter.

lambda_p.value

p-value of the Wald test for the estimated biomarker resistance to drug.

beta_p.value

p-values of the Wald test for the estimated regression coefficients.

Note

It is recommended to initialize the parameter values using the results from a corresponding reduced model (frailtyPenal for the recurrent and terminal part). See example.

Estimations of models with more than three random effects can be very long.

References

A. Krol, C. Tournigand, S. Michiels and V. RondeauS (2018). Multivariate joint frailty model for the analysis of nonlinear tumor kinetics and dynamic predictions of death. Statistics in Medicine.

A. Krol, L. Ferrer, JP. Pignon, C. Proust-Lima, M. Ducreux, O. Bouche, S. Michiels, V. Rondeau (2016). Joint Model for Left-Censored Longitudinal Data, Recurrent Events and Terminal Event: Predictive Abilities of Tumor Burden for Cancer Evolution with Application to the FFCD 2000-05 Trial. Biometrics 72(3) 907-16.

D. Rizopoulos (2012). Fast fitting of joint models for longitudinal and event time data using a pseudo-adaptive Gaussian quadrature rule. Computational Statistics and Data Analysis 56, 491-501.

L. Claret, P. Girard, P.M. Hoff, E. Van Cutsem, K.P. Zuideveld, K. Jorga, J. Fagerberg, R Bruno (2009). Model-based prediction of phase III overall survival in colorectal cancer on the basis of phase II tumor dynamics. Journal of Clinical Oncology 27(25), 4103-8.

See Also

plot.trivPenalNL,print.trivPenalNL,summary.trivPenalNL

Examples



## Not run: 

###--- Non-linear trivariate joint model for longitudinal data, ---###
###--- recurrent events and a terminal event ---###

data(colorectal)
data(colorectalLongi)

# No information on dose - creation of a dummy variable 
colorectalLongi$dose <- 1


# Parameters initialisation - estimation of a simplified model
# with two random effects (a frailty term and a random effect 
# related to biomarker growth (KG))
initial.model <- trivPenalNL(Surv(time0, time1, new.lesions) ~ cluster(id)
 + age + treatment + terminal(state), formula.terminalEvent =~ age + treatment, 
 biomarker = "tumor.size", formula.KG ~ 1, formula.KD ~ treatment, dose = "dose",
 time.biomarker = "year", data = colorectal, data.Longi =colorectalLongi, 
 random = "KG", id = "id", recurrentAG = TRUE, n.knots = 5, kappa = c(0.01, 2),
 method.GH = "Pseudo-adaptive")


# Trivariate joint model with initial values for parameters
# (computation takes around 40 minutes)

model <- trivPenalNL(Surv(time0, time1, new.lesions) ~ cluster(id) + age + treatment
 + terminal(state), formula.terminalEvent =~ age + treatment, biomarker = "tumor.size",
 formula.KG ~ 1, formula.KD ~ treatment, dose = "dose", time.biomarker = "year", 
 data = colorectal, data.Longi =colorectalLongi, random = c("y0", "KG"), id = "id", 
 init.B = c(-0.22, -0.16, -0.35, -0.19, 0.04, -0.41, 0.23), init.Alpha = 1.86,
 init.Eta = c(0.5, 0.57, 0.5, 2.34), init.Biomarker = c(1.24, 0.81, 1.07, -1.53),
 recurrentAG = TRUE, n.knots = 5, kappa = c(0.01, 2), method.GH = "Pseudo-adaptive")


## End(Not run)



[Package frailtypack version 3.6.2 Index]