Fit a multivariate frailty model for two types of recurrent events with a terminal event using a penalized likelihood estimation on the hazard function or a parametric estimation. Rightcensored data are allowed. Lefttruncated data and stratified analysis are not possible. Multivariate frailty models allow studying, with a joint model, three survival dependent processes for two types of recurrent events and a terminal event. Multivariate joint frailty models are applicable in mainly two settings. First, when focus is on the terminal event and we wish to account for the effect of previous endogenous recurrent event. Second, when focus is on a recurrent event and we wish to correct for informative censoring.
The multivariate frailty model for two types of recurrent events with a terminal event is (in the calendar or timetoevent timescale):
≤ft\{ \begin{array}{lll} r_{i}^{(1)}(tu_i,v_i) &= r_0^{(1)}(t)\exp({{β_1^{'}}}Z_{i}(t)+u_i) &\quad \mbox{(rec. of type 1)}\\ r_{i}^{(2)}(tu_i,v_i) &= r_0^{(2)}(t)\exp({{β_2^{'}}}Z_{i}(t)+v_i) &\quad \mbox{(rec. of type 2)}\\ λ_i(tu_i,v_i) &= λ_0(t)\exp({{β_3^{'}}}Z_{i}(t)+α_1u_i+α_2v_i) &\quad \mbox{(death)}\\ \end{array} \right.
where r_0^{(l)}(t), l\in{1,2} and λ_0(t) are respectively the recurrent and terminal event baseline hazard functions, and β_1,β_2,β_3 the regression coefficient vectors associated with Z_{i}(t) the covariate vector. The covariates could be different for the different event hazard functions and may be timedependent. We consider that death stops new occurrences of recurrent events of any type, hence given t>D, dN^{R(l)*}(t), l\in{1,2} takes the value 0. Thus, the terminal and the two recurrent event processes are not independent or even conditional upon frailties and covariates. We consider the hazard functions of recurrent events among individuals still alive. The three components in the above multivariate frailty model are linked together by two Gaussian and correlated random effects u_i,v_i:
(u_i,v_i)^{T}\sim\mathcal{N}≤ft({{0}},Σ_{uv}\right), with
Σ_{uv}=≤ft(\begin{array}{cc} θ_1 & ρ√{θ_1θ_2} \\ ρ√{θ_1θ_2}&θ_2 \end{array}\right)
Dependencies between these three types of event are taken into account by two correlated random effects and parameters θ_1,θ_2 the variance of the random effects and α_1,α_2 the coefficients for these random effects into the terminal event part. If α_1 and θ_1 are both significantly different from 0, then the recurrent events of type 1 and death are significantly associated (the sign of the association is the sign of α_1). If α_2 and θ_2 are both significantly different from 0, then the recurrent events of type 2 and death are significantly associated (the sign of the association is the sign of α_2). If ρ, the correlation between the two random effects, is significantly different from 0, then the recurrent events of type 1 and the recurrent events of type 2 are significantly associated (the sign of the association is the sign of ρ).
1 2 3 4  multivPenal(formula, formula.Event2, formula.terminalEvent, data,
initialize = TRUE, recurrentAG = FALSE, n.knots, kappa,
maxit = 350, hazard = "Splines", nb.int,
print.times = TRUE)

formula 
a formula object, with the response for the first recurrent event 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.Event2 
a formula object, with the response for the second recurrent event 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, with the response for the terminal event 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. 
data 
a 'data.frame' with the variables used in 'formula', 'formula.Event2' and 'formula.terminalEvent'. 
initialize 
Logical value to initialize regression coefficients and baseline hazard functions parameters. When the estimation is semiparametric with splines, this initialization produces also values for smoothing parameters (by cross validation). When initialization is requested, the program first fit two shared frailty models (for the two types of recurrent events) and a Cox proportional hazards model (for the terminal event). Default is TRUE. 
recurrentAG 
Logical value. Is AndersenGill 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 timedependent covariates. The default is FALSE. 
n.knots 
integer vector of length 3 (for the three outcomes) giving the number of knots to use. First is for the recurrent of type 1, second is for the recurrent of type 2 and third is for the terminal event hazard function. 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. Number of knots must be between 4 and 20. (See Note) 
kappa 
vector of length 3 (for the three outcomes) for positive smoothing parameters in the penalized likelihood estimation. First is for the recurrent of type 1, second is for the recurrent of type 2 and third is for the terminal event hazard function. The coefficient kappa of the integral of the squared second derivative of hazard function in the fit (penalized log likelihood). Initial values for the kappas can be obtained with the option "initialize=TRUE". We advise the user to identify several possible tuning parameters, note their defaults and look at the sensitivity of the results to varying them. Value required.(See Note) 
maxit 
maximum number of iterations for the Marquardt algorithm. Default is 350. 
hazard 
Type of hazard functions: "Splines" for semiparametric hazard functions with the penalized likelihood estimation, "Piecewiseper" for piecewise constant hazard function using percentile, "Piecewiseequi" for piecewise constant hazard function using equidistant intervals, "Weibull" for parametric Weibull function. Default is "Splines". 
nb.int 
An integer vector of length 3 (for the three outcomes). First is the Number of intervals (between 1 and 20) for the recurrent of type 1 parametric hazard functions ("Piecewiseper", "Piecewiseequi"). Second is the Number of intervals (between 1 and 20) for the recurrent of type 2 parametric hazard functions ("Piecewiseper", "Piecewiseequi"). Third is Number of intervals (between 1 and 20) for the death parametric hazard functions ("Piecewiseper", "Piecewiseequi") 
print.times 
a logical parameter to print iteration process. Default is TRUE. 
Parameters estimates of a multivariate joint frailty model, more generally a 'fraityPenal' object. Methods defined for 'frailtyPenal' objects are provided for print, plot and summary. The following components are included in a 'multivPenal' object for multivariate Joint frailty models.
b 
sequence of the corresponding estimation of the splines coefficients, the random effects variances, the coefficients of the frailties and the regression coefficients. 
call 
The code used for fitting the model. 
n 
the number of observations used in the fit. 
groups 
the number of subjects used in the fit. 
n.events 
the number of recurrent events of type 1 observed in the fit. 
n.events2 
the number of the recurrent events of type 2 observed in the fit. 
n.deaths 
the number of deaths observed in the fit. 
loglikPenal 
the complete marginal penalized loglikelihood in the semiparametric case. 
loglik 
the marginal loglikelihood in the parametric case. 
LCV 
the approximated likelihood crossvalidation criterion in the semi parametric case (with H minus the converged Hessian matrix, and l(.) the full loglikelihood. 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(.)) 
theta1 
variance of the frailty parameter for recurrences of type 1 (\bold{Var}(u_i)) 
theta2 
variance of the frailty parameter for recurrences of type 2 (\bold{Var}(v_i)) 
alpha1 
the coefficient associated with the frailty parameter u_i in the terminal hazard function. 
alpha2 
the coefficient associated with the frailty parameter v_i in the terminal hazard function. 
rho 
the correlation coefficient between u_i and v_i 
npar 
number of parameters. 
coef 
the regression coefficients. 
nvar 
A vector with the number of covariates of each type of hazard function as components. 
varH 
the variance matrix of all parameters before positivity constraint transformation (theta, the regression coefficients and the spline coefficients). Then, the delta method is needed to obtain the estimated variance parameters. 
varHIH 
the robust estimation of the variance matrix of all parameters (theta, the regression coefficients and the spline coefficients). 
formula 
the formula part of the code used for the model for the recurrent event. 
formula.Event2 
the formula part of the code used for the model for the second recurrent event. 
formula.terminalEvent 
the formula part of the code used for the model for the terminal event. 
x1 
vector of times for hazard functions of the recurrent events of type 1 are estimated. By default seq(0,max(time),length=99), where time is the vector of survival times. 
lam1 
matrix of hazard estimates and confidence bands for recurrent events of type 1. 
xSu1 
vector of times for the survival function of the recurrent event of type 1. 
surv1 
matrix of baseline survival estimates and confidence bands for recurrent events of type 1. 
x2 
vector of times for the recurrent event of type 2 (see x1 value). 
lam2 
the same value as lam1 for the recurrent event of type 2. 
xSu2 
vector of times for the survival function of the recurrent event of type 2 
surv2 
the same value as surv1 for the recurrent event of type 2. 
xEnd 
vector of times for the terminal event (see x1 value). 
lamEnd 
the same value as lam1 for the terminal event. 
xSuEnd 
vector of times for the survival function of the terminal event 
survEnd 
the same value as surv1 for the terminal event. 
type.of.Piecewise 
Type of Piecewise hazard functions (1:"percentile", 0:"equidistant"). 
n.iter 
number of iterations needed to converge. 
type.of.hazard 
Type of hazard functions (0:"Splines", "1:Piecewise", "2:Weibull"). 
n.knots 
a vector with number of knots for estimating the baseline functions. 
kappa 
a vector with the smoothing parameters in the penalized likelihood estimation corresponding to each baseline function as components. 
n.knots.temp 
initial value for the number of knots. 
zi 
splines knots. 
time 
knots for Piecewise hazard function for the recurrent event of type 1. 
timedc 
knots for Piecewise hazard function for the terminal event. 
time2 
knots for Piecewise hazard function for the recurrent event of type 2. 
noVar 
indicator vector for recurrent, death and recurrent 2 explanatory variables. 
nvarRec 
number of the recurrent of type 1 explanatory variables. 
nvarEnd 
number of death explanatory variables. 
nvarRec2 
number of the recurrent of type 2 explanatory variables. 
nbintervR 
Number of intervals (between 1 and 20) for the the recurrent of type 1 parametric hazard functions ("Piecewiseper", "Piecewiseequi"). 
nbintervDC 
Number of intervals (between 1 and 20) for the death parametric hazard functions ("Piecewiseper", "Piecewiseequi"). 
nbintervR2 
Number of intervals (between 1 and 20) for the the recurrent of type 2 parametric hazard functions ("Piecewiseper", "Piecewiseequi"). 
istop 
Vector of the convergence criteria. 
shape.weib 
shape parameters for the Weibull hazard function. 
scale.weib 
scale parameters for the Weibull hazard function. 
martingale.res 
martingale residuals for each cluster (recurrent of type 1). 
martingale2.res 
martingale residuals for each cluster (recurrent of type 2). 
martingaledeath.res 
martingale residuals for each cluster (death). 
frailty.pred 
empirical Bayes prediction of the first frailty term. 
frailty2.pred 
empirical Bayes prediction of the second frailty term. 
frailty.var 
variance of the empirical Bayes prediction of the first frailty term. 
frailty2.var 
variance of the empirical Bayes prediction of the second frailty term. 
frailty.corr 
Correlation between the empirical Bayes prediction of the two frailty. 
linear.pred 
linear predictor: uses Beta'X + ui in the multivariate frailty models. 
linear2.pred 
linear predictor: uses Beta'X + vi in the multivariate frailty models. 
lineardeath.pred 
linear predictor for the terminal part form the multivariate frailty models: Beta'X + alpha1 ui + alpha2 vi 
global_chisq 
Recurrent event of type 1: a vector with the values of each multivariate Wald test. 
dof_chisq 
Recurrent event of type 1: a vector with the degree of freedom for each multivariate Wald test. 
global_chisq.test 
Recurrent event of type 1: a binary variable equals to 0 when no multivariate Wald is given, 1 otherwise. 
p.global_chisq 
Recurrent event of type 1: a vector with the pvalues for each global multivariate Wald test. 
names.factor 
Recurrent event of type 1: Names of the "as.factor" variables. 
global_chisq2 
Recurrent event of type 2: a vector with the values of each multivariate Wald test. 
dof_chisq2 
Recurrent event of type 2: a vector with the degree of freedom for each multivariate Wald test. 
global_chisq.test2 
Recurrent event of type 2: a binary variable equals to 0 when no multivariate Wald is given, 1 otherwise. 
p.global_chisq2 
Recurrent event of type 2: a vector with the p_values for each global multivariate Wald test. 
names.factor2 
Recurrent event of type 2: Names of the "as.factor" variables. 
global_chisq_d 
Terminal event: a vector with the values of each multivariate Wald test. 
dof_chisq_d 
Terminal event: a vector with the degree of freedom for each multivariate Wald test. 
global_chisq.test_d 
Terminal event: a binary variable equals to 0 when no multivariate Wald is given, 1 otherwise. 
p.global_chisq_d 
Terminal event: a vector with the pvalues for each global multivariate Wald test. 
names.factordc 
Terminal event: Names of the "as.factor" variables. 
"kappa" (kappa[1], kappa[2] and kappa[3]) and "n.knots" (n.knots[1], n.knots[2] and n.knots[3]) are the arguments that the user has to change if the fitted model does not converge. "n.knots" takes integer values between 4 and 20. But with n.knots=20, the model will take a long time to converge. So, usually, begin first with n.knots=7, and increase it step by step until it converges. "kappa" only takes positive values. So, choose a value for kappa (for instance 10000), and if it does not converge, multiply or divide this value by 10 or 5 until it converges. Moreover, it may be useful to change the value of the initialize argument.
Mazroui Y., MathoulinPellissier S., MacGrogan G., Brouste V., Rondeau V. (2013). Multivariate frailty models for two types of recurrent events with an informative terminal event : Application to breast cancer data. Biometrical journal, 55(6), 866884.
terminal
,event2
,
print.multivPenal
,summary.multivPenal
,plot.multivPenal
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29  ## Not run:
### Multivariate Frailty model ###
data(dataMultiv)
# (computation takes around 60 minutes)
modMultiv.spli < multivPenal(Surv(TIMEGAP,INDICREC)~cluster(PATIENT)+v1+v2+
event2(INDICMETA)+terminal(INDICDEATH),formula.Event2=~v1+v2+v3,
formula.terminalEvent=~v1,data=dataMultiv,n.knots=c(8,8,8),
kappa=c(1,1,1),initialize=FALSE)
print(modMultiv.spli)
modMultiv.weib < multivPenal(Surv(TIMEGAP,INDICREC)~cluster(PATIENT)+v1+v2+
event2(INDICMETA)+terminal(INDICDEATH),formula.Event2=~v1+v2+v3,
formula.terminalEvent=~v1,data=dataMultiv,hazard="Weibull")
print(modMultiv.weib)
modMultiv.cpm < multivPenal(Surv(TIMEGAP,INDICREC)~cluster(PATIENT)+v1+v2+
event2(INDICMETA)+terminal(INDICDEATH),formula.Event2=~v1+v2+v3,
formula.terminalEvent=~v1,data=dataMultiv,hazard="Piecewiseper",
nb.int=c(6,6,6))
print(modMultiv.cpm)
## End(Not run)

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
Please suggest features or report bugs with the GitHub issue tracker.
All documentation is copyright its authors; we didn't write any of that.