Nothing
#' Fit a Shared or a Joint Frailty Generalized Survival Model
#'
#' @description{
#'
#' \if{html}{\bold{I. SHARED FRAILTY GENERALIZED SURVIVAL MODELS}
#'
#' Fit a gamma Shared Frailty Generalized Survival Model using
#' a parametric estimation, or a semi-parametric penalized likelihood estimation.
#' Right-censored data and strata (up to 6 levels) are allowed.
#' It allows to obtain a parametric or flexible semi-parametric smooth
#' hazard and survival functions.
#'
#' Each frailty term \eqn{u}\out{<sub>i</sub>} is assumed
#' to act multiplicatively on the hazard function, and to be drawn from a
#' Gamma distribution with unit mean and variance \eqn{\theta}.
#' Conditional on the frailty term, the hazard function for the
#' \eqn{j}\out{<sup>th</sup>} subject in the \eqn{i}\out{<sup>th</sup>} group
#' is then expressed by
#'
#' {\figure{gsm1.png}{options: width="70\%"}}
#'
#' where \bold{\eqn{x}}\out{<sub>ij</sub>}
#' is a collection of baseline covariates,
#' \bold{\eqn{\xi}} is a vector of parameters, and
#' \eqn{\lambda}\out{<sub>ij</sub>}
#' (\eqn{t} | \bold{\eqn{x}}\out{<sub>ij</sub>} ; \bold{\eqn{\xi}})
#' is the hazard function for an average value of the frailty.
#' The associated conditional survival function writes
#'
#' {\figure{gsm2.png}{options: width="70\%"}}
#'
#' where
#' \eqn{S}\out{<sub>ij</sub>}
#' (\eqn{t} | \bold{\eqn{x}}\out{<sub>ij</sub>} ; \bold{\eqn{\xi}})
#' designates the survival function for an average value of the frailty.
#'
#' Following Liu et al. (2017, 2018), the latter function is expressed in terms of
#' a link function \eqn{g}(.) and a linear predictor
#' \eqn{\eta}\out{<sub>ij</sub>}
#' (\eqn{t}, \bold{\eqn{x}}\out{<sub>ij</sub>}; \bold{\eqn{\xi}})
#' such that
#' \eqn{g}[\eqn{S}\out{<sub>ij</sub>}
#' (\eqn{t} | \bold{\eqn{x}}\out{<sub>ij</sub>} ; \bold{\eqn{\xi}})]
#' =
#' \eqn{\eta}\out{<sub>ij</sub>}
#' (\eqn{t}, \bold{\eqn{x}}\out{<sub>ij</sub>}; \bold{\eqn{\xi}}),
#' i.e.
#' \eqn{S}\out{<sub>ij</sub>}
#' (\eqn{t} | \bold{\eqn{x}}\out{<sub>ij</sub>} ; \bold{\eqn{\xi}})
#' =
#' \eqn{h}[\eqn{\eta}\out{<sub>ij</sub>}
#' (\eqn{t}, \bold{\eqn{x}}\out{<sub>ij</sub>}; \bold{\eqn{\xi}})]
#' with \eqn{h}() = \eqn{g}\out{<sup>-1</sup>}().
#'
#' The conditional survival function is finally modeled by
#'
#' {\figure{gsm3.png}{options: width="70\%"}}
#'
#' The table below summarizes the most commonly used (inverse) link functions and
#' their associated conditional survival, hazard and cumulative hazard functions.
#' PHM stands for "Proportional Hazards Model",
#' POM for "Proportional Odds Model,
#' PROM for "Probit Model" and AHM for "Additive Hazards Model".
#'
#' {\figure{gsm4.png}{options: width="100\%"}}
#'
#'
#' \bold{I.(a) Fully parametric case}
#'
#' In the fully parametric case, linear predictors considered are of the form
#'
#' {\figure{gsm5.png}{options: width="70\%"}}
#'
#' where \eqn{\rho > 0} is a shape parameter,
#' \eqn{\gamma > 0} a scale parameter,
#' \bold{\eqn{\beta}} a vector of regression coefficients,
#' and \bold{\eqn{\xi}} = (\eqn{\rho} ,\eqn{\gamma}, \bold{\eqn{\beta}}).
#' With the appropriate link function, such linear parametric predictors
#' make it possible to recover
#' a Weibull baseline survival function for PHMs and AHMs,
#' a log-logistic baseline survival function for POMs,
#' and a log-normal one for PROMs.
#'
#' \bold{I. (b) Flexible semi-parametric case}
#'
#' For PHM and AHM, a more flexible splines-based approach is proposed for
#' modeling the baseline hazard function and time-varying regression coefficients.
#' In this case, conditional on the frailty term \eqn{u}\out{<sub>i</sub>},
#' the hazard function for the \eqn{j}\out{<sup>th</sup>} subject
#' in the \eqn{i}\out{<sup>th</sup>} group is still expressed by
#' \eqn{\lambda}\out{<sub>ij</sub>}
#' (\eqn{t} | \bold{\eqn{x}}\out{<sub>ij</sub>}, \eqn{u}\out{<sub>i</sub>} ;
#' \bold{\eqn{\xi}})
#' = \eqn{u}\out{<sub>i</sub>}
#' \eqn{\lambda}\out{<sub>ij</sub>}
#' (\eqn{t} | \bold{\eqn{x}}\out{<sub>ij</sub>} ; \bold{\eqn{\xi}}),
#' but we have this time
#'
#' {\figure{gsm6.png}{options: width="70\%"}}
#'
#' The smoothness of baseline hazard function \eqn{\lambda}\out{<sub>0</sub>}()
#' is ensured by penalizing the log-likelihood by a term which has
#' large values for rough functions.
#'
#' Moreover, for parametric and flexible semi-parametric AHMs, the
#' log-likelihood is constrained to ensure the strict positivity of the hazards,
#' since the latter is not naturally guaranteed by the model.
#'
#'
#'
#' \bold{II. JOINT FRAILTY GENERALIZED SURVIVAL MODELS}
#'
#' Fit a gamma Joint Frailty Generalized Survival Model for recurrent and
#' terminal events using a parametric estimation,
#' or a semi-parametric penalized likelihood estimation.
#' Right-censored data and strata (up to 6 levels) for the recurrent event part
#' are allowed.
#' Joint frailty models allow studying, jointly, survival processes
#' of recurrent and terminal events, by considering the terminal event as an
#' informative censoring.
#'
#' This model includes an common patient-specific frailty term
#' \eqn{u}\out{<sub>i</sub>} for the two survival functions which will
#' take into account the unmeasured heterogeneity in the data,
#' associated with unobserved covariates.
#' The frailty term acts differently for the two survival functions
#' (\eqn{u}\out{<sub>i</sub>} for the recurrent survival function and
#' \eqn{u}\out{<sub>i</sub>}\out{<sup>α</sup>} for the death one).
#' The covariates could be different for the recurrent and terminal event parts.
#'
#' \bold{II.(a) Fully parametric case}
#'
#' For the \eqn{j}\out{<sup>th</sup>} recurrence (j=1,...,n\out{<sub>i</sub>})
#' and the \eqn{i}\out{<sup>th</sup>} patient (i=1,...,N),
#' the gamma Joint Frailty Generalized Survival Model
#' for recurrent event survival function
#' \eqn{S}\out{<sub>Rij</sub>}(.) and death survival function
#' \eqn{S}\out{<sub>Di</sub>}(.) is
#'
#' {\figure{gsm7.png}{options: width="70\%"}}
#'
#' - \eqn{\eta}\out{<sub>Rij</sub>} (resp. \eqn{\eta}\out{<sub>Di</sub>})
#' is the linear predictor for the recurrent (resp. terminal) event process.
#' The form of these linear predictors is the same as the one presented in I.(a).
#'
#' - \eqn{h}\out{<sub>R</sub>}(.) (resp. \eqn{h}\out{<sub>D</sub>}(.))
#' is the inverse link function associated with
#' recurrent events (resp. terminal event).
#'
#' - \bold{\eqn{x}}\out{<sub>Rij</sub>} and \bold{\eqn{x}}\out{<sub>Di</sub>}
#' are two vectors of baseline covariates associated with
#' recurrent and terminal events.
#'
#' - \bold{\eqn{\xi}}\out{<sub>R</sub>} and \bold{\eqn{\xi}}\out{<sub>D</sub>}
#' are the parameter vectors for recurrent and terminal events.
#'
#' - \eqn{\alpha} is a parameter allowing more flexibility in the association
#' between recurrent and terminal events processes.
#'
#' - The random frailties \eqn{u}\out{<sub>i</sub>} are still assumed iid and
#' drown from a \eqn{\Gamma}(1/\eqn{\theta},1/\eqn{\theta}).
#'
#'
#' \bold{II.(b) Flexible semi-parametric case}
#'
#' If one chooses to fit a PHM or an AHM for recurrent and/or terminal events,
#' a splines-based approach for modeling baseline hazard functions
#' and time-varying regression coefficients is still available.
#' In this approach, the submodel for recurrent events is expressed as
#' \eqn{\lambda}\out{<sub>Rij</sub>}
#' (\eqn{t} | \bold{\eqn{x}}\out{<sub>Rij</sub>}, \eqn{u}\out{<sub>i</sub>} ;
#' \bold{\eqn{\xi}}\out{<sub>R</sub>})
#' = \eqn{u}\out{<sub>i</sub>}
#' \eqn{\lambda}\out{<sub>Rij</sub>}
#' (\eqn{t} | \bold{\eqn{x}}\out{<sub>Rij</sub>} ;
#' \bold{\eqn{\xi}}\out{<sub>R</sub>}), where
#'
#' {\figure{gsm8.png}{options: width="70\%"}}
#'
#' The submodel for terminal event is expressed as
#' \eqn{\lambda}\out{<sub>Di</sub>}
#' (\eqn{t} | \bold{\eqn{x}}\out{<sub>Di</sub>}, \eqn{u}\out{<sub>i</sub>} ;
#' \bold{\eqn{\xi}}\out{<sub>D</sub>})
#' = \eqn{u}\out{<sub>i</sub>}\out{<sup>α</sup>}
#' \eqn{\lambda}\out{<sub>Di</sub>}
#' (\eqn{t} | \bold{\eqn{x}}\out{<sub>Di</sub>} ;
#' \bold{\eqn{\xi}}\out{<sub>D</sub>}), where
#'
#' {\figure{gsm9.png}{options: width="70\%"}}
#'
#' Baseline hazard functions
#' \eqn{\lambda}\out{<sub>R0</sub>}(.) and \eqn{\lambda}\out{<sub>D0</sub>}(.)
#' are estimated using cubic M-splines (of order 4)
#' with positive coefficients, and the time-varying coefficients
#' \eqn{\beta}\out{<sub>R</sub>}(.) and \eqn{\beta}\out{<sub>D</sub>}(.)
#' are estimated using B-splines of order q.
#'
#' The smoothness of baseline hazard functions
#' is ensured by penalizing the log-likelihood by two terms
#' which has large values for rough functions.
#'
#' Moreover,
#' if one chooses an AHM for recurrent and/or terminal event submodel,
#' the log-likelihood is constrained to ensure
#' the strict positivity of the hazards,
#' since the latter is not naturally guaranteed by the model.
#'
#'
# There is a possibility to use a weighted penalized maximum likelihood
# approach for nested case-control design, in which risk set sampling is
# performed based on a single outcome (Jazic et al., \emph{Submitted}).
#
# General Joint Frailty model Fit a general joint frailty model for recurrent
# and terminal events considering two independent frailty terms. The frailty
# term \eqn{u}\out{<sub>i</sub>} represents the unobserved association between recurrences and
# death. The frailty term \eqn{v}\out{<sub>i</sub>} is specific to the recurrent event rate.
# Thus, the general joint frailty model is:
#
# {\figure{frailtymodel9.png}{options: width="90\%"}}
#
# where the \eqn{iid} random effects
# \bold{\eqn{u}\out{<sub>i</sub>}} \out{~} \bold{\eqn{\Gamma}(1/\eqn{\theta},1/\eqn{\theta})} and the
# \eqn{iid} random effects
# \bold{\eqn{v}\out{<sub>i</sub>}} \out{~} \bold{\eqn{\Gamma}(1/\eqn{\eta},1/\eqn{\eta})} are independent
# from each other. The joint model is fitted using a penalized likelihood
# estimation on the hazard. Right-censored data and time-varying covariates
# \bold{\eqn{Z}\out{<sub>i</sub>}(t)} are allowed.
#
# \bold{Nested Frailty model}
#
# \bold{\emph{Data should be ordered according to cluster and subcluster}}
#
# Fit a nested frailty model using a Penalized Likelihood on the hazard
# function or using a parametric estimation. Nested frailty models allow
# survival studies for hierarchically clustered data by including two iid
# gamma random effects. Left-truncated and right-censored data are allowed.
# Stratification analysis is allowed (maximum of strata = 2).
# The hazard function conditional on the two frailties \eqn{v}\out{<sub>i</sub>} and
# \eqn{w}\out{<sub>ij</sub>} for the k\out{<sup>th</sup>} individual of the j\out{<sup>th</sup>} subgroup of
# the i\out{<sup>th</sup>} group is :
#
# {\figure{frailtymodel10.png}{options: width="80\%"}}
#
# where \eqn{\lambda}\out{<sub>0</sub>}(t) is the baseline hazard function, \eqn{X}\out{<sub>ijk</sub>}
# denotes the covariate vector and \eqn{\beta} the corresponding vector of
# regression parameters.
#
# \bold{Joint Nested Frailty Model}
#
# Fit a joint model for recurrent and terminal events using a penalized
# likelihood on the hazard functions or a parametric estimation.
# Right-censored data are allowed but left-truncated data and stratified
# analysis are not allowed.
#
# Joint nested frailty models allow studying, jointly, survival processes of
# recurrent and terminal events for hierarchically clustered data, by
# considering the terminal event as an informative censoring and by including
# two iid gamma random effects.
#
# The joint nested frailty model includes two shared frailty terms, one for
# the subgroup (\eqn{u}\out{<sub>fi</sub>}) and one for the group (\eqn{w}\out{<sub>f</sub>}) into the
# hazard functions. This random effects account the heterogeneity in the data,
# associated with unobserved covariates. The frailty terms act differently for
# the two rates (\eqn{u}\out{<sub>fi</sub>}, \eqn{w}\out{<sub>f</sub>}\out{<sup>\eqn{\xi}</sup>} for the recurrent rate and
# \eqn{u}\out{<sub>fi</sub>}\out{<sup>\eqn{\alpha}</sup>}, \eqn{w}\out{<sub>i</sub>} for the terminal event rate). The covariates
# could be different for the recurrent rate and death rate.
#
# For the j\out{<sup>th</sup>} recurrence (j = 1, ..., \eqn{n}\out{<sub>i</sub>}) of the i\out{<sup>th</sup>}
# individual (i = 1, ..., \eqn{m}\out{<sub>f</sub>}) of the \eqn{f}\out{<sup>th</sup>} group (f = 1,...,
# n), the joint nested gamma frailty model for recurrent event hazard function
# \eqn{r}\out{<sub>fij</sub>}(.) and for terminal event hazard function \eqn{\lambda}\out{<sub>fi</sub>}
# is :
#
# {\figure{frailtymodel11.png}{options: width="90\%"}}
#
# where \eqn{r}\out{<sub>0</sub>}(resp. \eqn{\lambda}\out{<sub>0</sub>}) is the recurrent (resp.
# terminal) event baseline hazard function, \eqn{\beta} (resp. \eqn{\gamma})
# the regression coefficient vector, \bold{\eqn{X}\out{<sub>fij</sub>}}(t) the covariates
# vector. The random effects are \eqn{\omega}\out{<sub>f</sub>} \out{<span>~</span>} \bold{\eqn{\Gamma}(1/\eqn{\eta},1/\eqn{\eta})}
# and \eqn{u}\out{<sub>fi</sub>} \out{<span>~</span>} \bold{\eqn{\Gamma}(1/\eqn{\theta},1/\eqn{\theta})}.
#' }
#' \if{latex}{\bold{Shared Frailty model}
#'
#' Fit a shared gamma or log-normal frailty model using a semiparametric
#' Penalized Likelihood estimation or parametric estimation on the hazard
#' function. Left-truncated, right-censored data, interval-censored data and
#' strata (up to 6 levels) are allowed. It allows to obtain a non-parametric
#' smooth hazard of survival function. This approach is different from the
#' partial penalized likelihood approach of Therneau et al.
#'
#' The hazard function, conditional on the frailty term \eqn{\omega_i}, of a
#' shared gamma frailty model for the \eqn{j^{th}} subject in the \eqn{i^{th}}
#' group:
#'
#' \deqn{\lambda_{ij}(t|\omega_i)=\lambda_0(t)\omega_i\exp(\bold{\beta^{'}Z_{ij}})}
#'
#' \deqn{\omega_i\sim\Gamma\left(\frac{1}{\theta},\frac{1}{\theta}\right)
#' \hspace{0.5cm} \bold{E}(\omega_i)=1
#' \hspace{0.5cm}\bold{Var}(\omega_i)=\theta}
#'
#' where \eqn{\lambda_0(t)} is the baseline hazard function, \eqn{\bold{\beta}}
#' the vector of the regression coefficient associated to the covariate vector
#' \eqn{\bold{Z_{ij}}} for the \eqn{j^{th}} individual in the \eqn{i^{th}}
#' group.
#'
#' Otherwise, in case of a shared log-normal frailty model, we have for the
#' \eqn{j^{th}} subject in the \eqn{i^{th}} group:
#'
#' \deqn{\lambda_{ij}(t|\eta_i)=\lambda_0(t)\exp(\eta_i+\bold{\beta^{'}Z_{ij}})}
#'
#' \deqn{\eta_i\sim N(0,\sigma^2)}
#'
#' From now on, you can also consider time-varying effects covariates in your
#' model, see \code{timedep} function for more details.
#'
#' \bold{Joint Frailty model}
#'
#' Fit a joint either with gamma or log-normal frailty model for recurrent and
#' terminal events using a penalized likelihood estimation on the hazard
#' function or a parametric estimation. Right-censored data and strata (up to 6
#' levels) for the recurrent event part are allowed. Left-truncated data is not
#' possible. Joint frailty models allow studying, jointly, survival processes
#' of recurrent and terminal events, by considering the terminal event as an
#' informative censoring.
#'
#' There is two kinds of joint frailty models that can be fitted with
#' \code{GenfrailtyPenal} :
#'
#' - The first one (Rondeau et al. 2007) includes a common frailty term to the
#' individuals \eqn{(\omega_i)} for the two rates which will take into account
#' the heterogeneity in the data, associated with unobserved covariates. The
#' frailty term acts differently for the two rates ( \eqn{\omega_i} for the
#' recurrent rate and \eqn{\omega_i^{\alpha}} for the death rate). The
#' covariates could be different for the recurrent rate and death rate.
#'
#' For the \eqn{j^{th}}{j^th} recurrence \eqn{(j=1,...,n_i)} and the
#' \eqn{i^{th}}{i^th} subject \eqn{(i=1,...,G)}, the joint gamma frailty model
#' for recurrent event hazard function \eqn{r_{ij}(.)} and death rate
#' \eqn{\lambda_i(.)} is :
#'
#' \deqn{\left\{ \begin{array}{ll}
#' r_{ij}(t|\omega_i)=\omega_ir_0(t)\exp(\bold{\beta_1^{'}Z_i(t)}) &
#' \mbox{(Recurrent)} \\
#' \lambda_i(t|\omega_i)=\omega_i^{\alpha}\lambda_0(t)\exp(\bold{\beta_2^{'}Z_i(t)})
#' & \mbox{(Death)} \\ \end{array} \right. }
#'
#' where \eqn{r_0(t)} (resp. \eqn{\lambda_0(t)}) is the recurrent (resp.
#' terminal) event baseline hazard function, \eqn{\bold{\beta_1}} (resp.
#' \eqn{\bold{\beta_2}}) the regression coefficient vector, \eqn{\bold{Z_i(t)}}
#' the covariate vector. The random effects of frailties
#' \eqn{\omega_i\sim\bold{\Gamma}(\frac{1}{\theta},\frac{1}{\theta})} and are
#' iid.
#'
#' The joint log-normal frailty model will be :
#'
#' \deqn{\left\{ \begin{array}{ll}
#' r_{ij}(t|\eta_i)=r_0(t)\exp(\eta_i+\bold{\beta_1^{'}Z_i(t)}) &
#' \mbox{(Recurrent)} \\ \lambda_i(t|\eta_i)=\lambda_0(t)\exp(\alpha
#' \eta_i+\bold{\beta_2^{'}Z_i(t)}) & \mbox{(Death)} \\ \end{array} \right. }
#'
#' where \deqn{\eta_i\sim N(0,\sigma^2)}
#'
#' - The second one (Rondeau et al. 2011) is quite similar but the frailty term
#' is common to the individuals from a same group. This model is useful for the
#' joint modelling two clustered survival outcomes. This joint models have been
#' developed for clustered semi-competing events. The follow-up of each of the
#' two competing outcomes stops when the event occurs. In this case, j is for
#' the subject and i for the cluster.
#'
#' \deqn{\left\{ \begin{array}{ll}
#' r_{ij}(t|u_i)=u_ir_0(t)\exp(\bold{\beta_1^{'}Z_{ij}(t)}) & \mbox{(Time to
#' event)} \\
#' \lambda_{ij}(t|u_i)=u_i^{\alpha}\lambda_0(t)\exp(\bold{\beta_2^{'}Z_{ij}(t)})
#' & \mbox{(Death)} \\ \end{array} \right. }
#'
#' It should be noted that in these models it is not recommended to include
#' \eqn{\alpha} parameter as there is not enough information to estimate it and
#' thus there might be convergence problems.
#'
#' In case of a log-normal distribution of the frailties, we will have :
#'
#' \deqn{\left\{ \begin{array}{ll}
#' r_{ij}(t|v_i)=r_0(t)\exp(v_i+\bold{\beta_1^{'}Z_{ij}(t)}) & \mbox{(Time to
#' event)} \\ \lambda_{ij}(t|v_i)=\lambda_0(t)\exp(\alpha
#' v_i+\bold{\beta_2^{'}Z_{ij}(t)}) & \mbox{(Death)} \\ \end{array} \right. }
#'
#' where \deqn{v_i\sim N(0,\sigma^2)}
#'
#' This joint frailty model can also be applied to clustered recurrent events
#' and a terminal event (example on "readmission" data below).
#'
#' From now on, you can also consider time-varying effects covariates in your
#' model, see \code{timedep} function for more details.
#'
#' There is a possibility to use a weighted penalized maximum likelihood
#' approach for nested case-control design, in which risk set sampling is
#' performed based on a single outcome (Jazic et al., \emph{Submitted}).
#'
#' General Joint Frailty model Fit a general joint frailty model for recurrent
#' and terminal events considering two independent frailty terms. The frailty
#' term \eqn{u_i} represents the unobserved association between recurrences and
#' death. The frailty term \eqn{v_i} is specific to the recurrent event rate.
#' Thus, the general joint frailty model is:
#'
#' \eqn{\left\{ \begin{array}{ll}
#' r_{ij}(t|u_i,v_i)=u_iv_ir_0(t)\exp(\bold{\beta_1^{'}Z_{ij}(t)})
#' =u_iv_ir_{ij}(t) & \mbox{(Recurrent)} \\
#' \lambda_{i}(t|u_i)=u_i\lambda_0(t)\exp(\bold{\beta_1^{'}Z_{i}(t)}) = u_i
#' \lambda_{i}(t) & \mbox{(Death)} \\ \end{array} \right. }
#'
#' where the \eqn{iid} random effects
#' \eqn{\bold{u_i}\sim\Gamma(\frac{1}{\theta},\frac{1}{\theta})} and the
#' \eqn{iid} random effects
#' \eqn{\bold{v_i}\sim\Gamma(\frac{1}{\eta},\frac{1}{\eta})} are independent
#' from each other. The joint model is fitted using a penalized likelihood
#' estimation on the hazard. Right-censored data and time-varying covariates
#' \eqn{\bold{Z}_i(t)} are allowed.
#'
#' \bold{Nested Frailty model}
#'
#' \bold{\emph{Data should be ordered according to cluster and subcluster}}
#'
#' Fit a nested frailty model using a Penalized Likelihood on the hazard
#' function or using a parametric estimation. Nested frailty models allow
#' survival studies for hierarchically clustered data by including two iid
#' gamma random effects. Left-truncated and right-censored data are allowed.
#' Stratification analysis is allowed (maximum of strata = 2).
#'
#' The hazard function conditional on the two frailties \eqn{v_i} and
#' \eqn{w_{ij}} for the \eqn{k^{th}} individual of the \eqn{j^{th}} subgroup of
#' the \eqn{i^{th}} group is :
#'
#' \deqn{\left\{ \begin{array}{ll}
#' \lambda_{ijk}(t|v_i,w_{ij})=v_iw_{ij}\lambda_0(t)exp(\bold{\beta^{'}X_{ijk}})
#' \\ v_i\sim\Gamma\left(\frac{1}{\alpha},\frac{1}{\alpha}\right)
#' \hspace{0.05cm}i.i.d. \hspace{0.2cm} \bold{E}(v_i)=1
#' \hspace{0.2cm}\bold{Var}(v_i)=\alpha \\
#' w_{ij}\sim\Gamma\left(\frac{1}{\eta},\frac{1}{\eta}\right)\hspace{0.05cm}i.i.d.
#' \hspace{0.2cm} \bold{E}(w_{ij})=1 \hspace{0.2cm} \bold{Var}(w_{ij})=\eta
#' \end{array} \right. }
#'
#' where \eqn{\lambda_0(t)} is the baseline hazard function, \eqn{X_{ijk}}
#' denotes the covariate vector and \eqn{\beta} the corresponding vector of
#' regression parameters.
#'
#' \bold{Joint Nested Frailty Model}
#'
#' Fit a joint model for recurrent and terminal events using a penalized
#' likelihood on the hazard functions or a parametric estimation.
#' Right-censored data are allowed but left-truncated data and stratified
#' analysis are not allowed.
#'
#' Joint nested frailty models allow studying, jointly, survival processes of
#' recurrent and terminal events for hierarchically clustered data, by
#' considering the terminal event as an informative censoring and by including
#' two iid gamma random effects.
#'
#' The joint nested frailty model includes two shared frailty terms, one for
#' the subgroup (\eqn{u_{fi}}) and one for the group (\eqn{w_f}) into the
#' hazard functions. This random effects account the heterogeneity in the data,
#' associated with unobserved covariates. The frailty terms act differently for
#' the two rates (\eqn{u_{fi}}, \eqn{w_f^\xi} for the recurrent rate and
#' \eqn{u_{fi}^\alpha, {w_i}} for the terminal event rate). The covariates
#' could be different for the recurrent rate and death rate.
#'
#' For the \eqn{j^{th}} recurrence (j = 1, ..., \eqn{n_i}) of the \eqn{i^{th}}
#' individual (i = 1, ..., \eqn{m_f}) of the \eqn{f^{th}} group (f = 1, ...,
#' n), the joint nested gamma frailty model for recurrent event hazard function
#' \eqn{r_{fij}}(.) and for terminal event hazard function \eqn{\lambda_{fi}}
#' is :
#'
#' \deqn{\left\{ \begin{array}{ll} r_{fij}(t|\omega_f, u_{fi}, \bold{X_{fij}})=
#' r_0(t) u_{fi} \omega_f^\xi \exp(\bold{\beta'} \bold{X_{fij}}) &
#' \mbox{(Recurrent)} \\ \lambda_{fi}(t|\omega_f, u_{fi},
#' \bold{X_{fij}})=\lambda_0(t)u_{fi}^\alpha \omega_f \exp(\bold{\gamma'}
#' \bold{X_{fi}}) & \mbox{(Death)} \\ \end{array} \right. }
#'
#' where \eqn{r_0(t)}(resp. \eqn{\lambda_0(t)}) is the recurrent (resp.
#' terminal) event baseline hazard function, \eqn{\beta} (resp. \eqn{\gamma})
#' the regression coefficient vector, \eqn{\bold{X_{fij}}(t)} the covariates
#' vector. The random effects are \deqn{\omega_f \sim \Gamma \left(
#' \frac{1}{\eta}, \frac{1}{\eta}\right)} and \deqn{ u_{fi} \sim \Gamma \left(
#' \frac{1}{\theta}, \frac{1}{\theta} \right)}
#' }
#' }
#'
#' @details{
#' \bold{TYPICAL USES}
#'
#' For a Generalized Survival Model:
#' \preformatted{GenfrailtyPenal(
#' formula=Surv(time,event)~var1+var2,
#' data, family, \dots)}
#'
#' For a Shared Frailty Generalized Survival Model:
#' \preformatted{GenfrailtyPenal(
#' formula=Surv(time,event)~cluster(group)+var1+var2,
#' data, family, \dots)}
#'
#' For a Joint Frailty Generalized Survival Model:
#' \preformatted{GenfrailtyPenal(
#' formula=Surv(time,event)~cluster(group)+var1+var2+var3+terminal(death),
#' formula.terminalEvent= ~var1+var4,
#' data, family, \dots)}
#'
#for a joint model for clustered data
#\preformatted{GenfrailtyPenal(Surv(time,event)~cluster(group)+num.id(group2)+
#var1+var2+var3+terminal(death), formula.terminalEvent=~var1+var4, data,
#\dots)}
#
#for a joint model for data from nested case-control studies
#\preformatted{GenfrailtyPenal(Surv(time,event)~cluster(group)+num.id(group2)+
#var1+var2+var3+terminal(death)+wts(wts.ncc),
#formula.terminalEvent=~var1+var4, data, \dots)}
#
#for a nested model
#\preformatted{GenfrailtyPenal(Surv(time,event)~cluster(group)+subcluster(sbgroup)+
#var1+var2, data, \dots)}
#
#for a joint nested frailty model
#\preformatted{GenfrailtyPenal(Surv(time,event)~cluster(group)+subcluster(sbgroup)+
#var1+var2++terminal(death), formula.terminalEvent=~var1+var4, data, \dots)}
#
#'
#' \bold{OPTIMIZATION ALGORITHM}
#'
#' The estimated parameters are obtained using the robust Marquardt algorithm
#' (Marquardt, 1963) which is a combination between a Newton-Raphson algorithm
#' and a steepest descent algorithm. The iterations are stopped when
#' the difference between two consecutive log-likelihoods is small
#' (\eqn{<10}\out{<sup>-3</sup>}),
#' the estimated coefficients are stable
#' (consecutive values \eqn{<10}\out{<sup>-3</sup>},
#' and the gradient small enough (\eqn{<10}\out{<sup>-3</sup>}).
#' When the frailty variance is small, numerical problems may arise.
#' To solve this problem, an alternative formula of the penalized log-likelihood
#' is used (see Rondeau, 2003 for further details).
#' For Proportional Hazards and Additive Hazards submodels,
#' cubic M-splines of order 4 can be used to estimate the hazard function.
#' In this case, I-splines (integrated M-splines) are used to compute the
#' cumulative hazard function.
#'
#' The inverse of the Hessian matrix is the variance estimator.
#' To deal with the positivity constraint of the variance component and the
#' spline coefficients, a squared transformation is used and the standard errors
#' are computed by the \eqn{\Delta}-method (Knight & Xekalaki, 2000).
#The smooth parameter can be chosen by maximizing a likelihood cross validation
#criterion (Joly and other, 1998).
#' The integrations in the full log likelihood are evaluated using
#' Gaussian quadrature. Laguerre polynomials with 20 points are used to treat
#' the integrations on \eqn{[0,\infty[}.
#'
#'
#' \bold{INITIAL VALUES}
#'
#' In case of a shared frailty model,
#' the splines and the regression coefficients are initialized to 0.1.
#' The program fits, firstly, an adjusted Cox model to give new initial values
#' for the splines and the regression coefficients.
#' The variance of the frailty term \eqn{\theta} is initialized to 0.1.
#' Then, a shared frailty model is fitted.
#'
#' In case of a joint frailty model,
#' the splines and the regression coefficients are initialized to 0.5.
#' The program fits firstly, an adjusted Cox model to have new initial values
#' for the splines and the regression coefficients.
#' The variance of the frailty term \eqn{\theta} and the association parameter
#' \eqn{\alpha} are initialized to 1.
#' Then, a joint frailty model is fitted.
#'
#In case of a general joint frailty model we need to initialize the
#\code{jointGeneral} logical value to \code{TRUE}.
#
#In case of a nested model, the program fits an adjusted Cox model to provide
#new initial values for the regression and the splines coefficients. The
#variances of the frailties are initialized to 0.1. Then, a shared frailty
#model with covariates with only subgroup frailty is fitted to give a new
#initial value for the variance of the subgroup frailty term. Then, a shared
#frailty model with covariates and only group frailty terms is fitted to give
#a new initial value for the variance of the group frailties. In a last step,
#a nested frailty model is fitted.
#
#In case of a joint nested model, the splines and the regression coefficients
#are initialized to 0.5 and the variances of the frailty terms \eqn{\eta} and
#\eqn{\xi} are initialized to 1. If the option \code{'initialize'} is
#\code{TRUE}, the program fits a joint frailty model to provide initial
#values for the splines, covariates coefficients, variance \eqn{\theta} of
#the frailty terms and \eqn{\alpha}. The variances of the second frailty term
#(\eqn{\eta}) and the second coefficient \eqn{\xi} are initialized to 1.
#Then, a joint nested frailty model is fitted.
#
#\bold{NCC DESIGN}
#
#It is possible to fit a joint frailty model for data from nested
#case-control studies using the approach of weighted penalized maximum
#likelihood. For this model, only splines can be used for baseline hazards
#and no time-varying effects of covariates can be included. To accommodate
#the nested case-control design, the formula for the recurrent events should
#simply include the special term wts(wts.ncc), where wts.ncc refers to a
#column of prespecified weights in the data set for every observation. For
#details, see Jazic et al., \emph{Submitted} (available on request from the
#package authors).
#' }
#'
#' @aliases GenfrailtyPenal
#' @usage
#'
#' GenfrailtyPenal(formula, formula.terminalEvent, data, recurrentAG = FALSE,
#' family, hazard = "Splines", n.knots, kappa, betaknots = 1, betaorder = 3,
#' RandDist = "Gamma", init.B, init.Theta, init.Alpha, Alpha, maxit = 300,
#' nb.gh, nb.gl, LIMparam = 1e-3, LIMlogl = 1e-3, LIMderiv = 1e-3, print.times = TRUE,
#' cross.validation, jointGeneral, nb.int, initialize, init.Ksi, Ksi, init.Eta)
#' @param formula A formula object, with the response on the left of a
#' \eqn{\sim} operator, and the terms on the right. The response must be a
#' survival object as returned by the '\code{Surv}' function
#' like in survival package. Interactions are possible using ' * ' or ' : '.
#' @param formula.terminalEvent Only for joint frailty models: a formula object,
#' only requires terms on the right to indicate which variables are used for
#' the terminal event. Interactions are possible using ' * ' or ' : '.
#' @param data A 'data.frame' with the variables used in '\code{formula}'.
#' @param 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.
#' @param family Type of Generalized Survival Model to fit.
#' \code{"PH"} for a proportional hazards model,
#' \code{"AH"} for an additive hazards model,
#' \code{"PO"} for a proportional odds model and
#' \code{"probit"} for a probit model.
#' A vector of length 2 is expected for joint models
#' (e.g., \code{family=c("PH","PH")}).
#' @param hazard Type of hazard functions:
#' \code{"Splines"} for semi-parametric hazard functions using equidistant
#' intervals, or \code{"parametric"} for parametric distribution functions.
#' In case of \code{family="PH"} or \code{family="AH"},
#' the \code{"parametric"} option corresponds to a Weibull distribution.
#' In case of \code{family="PO"} and \code{family="probit"},
#' the \code{"parametric"} option corresponds to a log-logistic
#' and a log-normal distribution, respectively.
#' So far, the \code{"Splines"} option is only available for PH and AH submodels.
#' Default is \code{"Splines"}.
#' @param n.knots Integer giving the number of knots to use. Value required in
#' the penalized likelihood estimation. It corresponds to the \code{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 \code{k} (i.e. \code{n.knots=k}),
#' then the number of interior knots is \code{k-2} and the number of splines is
#' \code{(k-2)+order}. Number of knots must be between 4 and 20. (See Note)
#' @param kappa Positive smoothing parameter in the penalized likelihood
#' estimation. The coefficient kappa tunes the intensity of the penalization
#' (the integral of the squared second derivative of hazard function).
#' In a stratified shared model, this argument must be a vector
#' with kappas for both strata. In a stratified joint model, this argument
#' must be a vector with kappas for both strata for recurrent events plus one
#' kappa for terminal event.
#' 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).
#' @param betaknots Number of inner knots used for the
#' B-splines time-varying coefficient estimation. Default is 1.
#' See '\code{timedep}' function for more details.
#' @param betaorder Order of the B-splines used for the
#' time-varying coefficient estimation.
#' Default is cubic B-splines (\code{order=3}).
#' See '\code{timedep}' function for more details.
#' Not implemented for Proportional Odds and Probit submodels.
#' @param RandDist Type of random effect distribution:
#' \code{"Gamma"} for a gamma distribution,
#' and \code{"LogN"} for a log-normal distribution (not implemented yet).
#' Default is \code{"Gamma"}.
#' @param init.B A 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 and then
#' to the terminal event (interactions in the end of each component). Default
#' is 0.1 for each (for Generalized Survival and Shared Frailty Models)
#' or 0.5 (for Generalized Joint Frailty Models).
#' @param init.Theta Initial value for frailty variance.
#' @param init.Alpha Only for Generalized Joint Frailty Models:
#' initial value for parameter alpha.
#' @param Alpha Only for Generalized Joint Frailty Models:
#' input "None" so as to fit a joint model without the parameter alpha.
#' @param maxit Maximum number of iterations for the Marquardt algorithm.
#' Default is 300
#' @param nb.gh Number of nodes for the Gaussian-Hermite quadrature.
#' It can be chosen among 5, 7, 9, 12, 15, 20 and 32.
#' The default is 20 if \code{hazard="Splines"}, 32 otherwise.
#' @param nb.gl Number of nodes for the Gaussian-Laguerre quadrature.
#' It can be chosen between 20 and 32.
#' The default is 20 if \code{hazard="Splines"}, 32 otherwise.
#' @param LIMparam Convergence threshold of the Marquardt algorithm for the
#' parameters (see Details), \eqn{10}\out{<sup>-3</sup>} by default.
#' @param LIMlogl Convergence threshold of the Marquardt algorithm for the
#' log-likelihood (see Details), \eqn{10}\out{<sup>-3</sup>} by default.
#' @param LIMderiv Convergence threshold of the Marquardt algorithm for the
#' gradient (see Details), \eqn{10}\out{<sup>-3</sup>} by default.
#' @param print.times A logical parameter to print iteration process. Default
#' is TRUE.
#' @param cross.validation Not implemented yet for the generalized settings.
#Logical value. Is cross validation procedure used
#for estimating smoothing parameter in the penalized likelihood estimation?
#If so a search of the smoothing parameter using cross validation is done,
#with kappa as the seed. The cross validation is not implemented for several
#strata, neither for interval-censored data. The cross validation has been
#implemented for a Cox proportional hazard model, with no covariates. The
#default is FALSE.
#' @param jointGeneral Not implemented yet for the generalized settings.
#Logical value. Does the model include two independent
#random effects? If so, this will fit a general joint frailty model with an
#association between the recurrent events and a terminal event (explained by
#the variance \eqn{\theta}) and an association amongst the recurrent events
#(explained by the variance \eqn{\eta}).
#' @param nb.int Not implemented yet for the generalized settings.
#Number of time intervals (between 1 and 20) for the parametric
#hazard functions ("Piecewise-per", "Piecewise-equi"). In a joint model, you
#need to specify a number of time interval for both recurrent hazard function
#and the death hazard function (vector of length 2).
#' @param initialize Not implemented yet for the generalized settings.
#Option \code{TRUE} indicates fitting an appropriate standard joint frailty
#model (without group effect, only the subgroup effect) to provide initial
#values for the joint nested model. Default is \code{TRUE}.
#' @param init.Ksi Not implemented yet for the generalized settings.
#Initial value for parameter \eqn{\xi}.
#' @param init.Eta Not implemented yet for the generalized settings.
#(only for general joint and joint nested frailty models) !!
#Only for general joint and joint nested frailty models :
#initial value for the variance \eqn{\eta} of the frailty \eqn{v_i} (general
#joint model) and of the frailty \eqn{\omega_i} (joint nested frailty model).
#' @param Ksi Not implemented yet for the generalized settings.
#input \code{"None"} indicates a joint nested frailty model without
#the parameter \eqn{\xi}.
#' @return
#'
#' The following components are included in a 'frailtyPenal' object for each
#' model.
#'
#' \item{b}{Sequence of the corresponding estimation of the coefficients for
#' the hazard functions (parametric or semiparametric), the random effects
#' variances and the regression coefficients.}
#' \item{call}{The code used for the model.}
#' \item{formula}{The formula part of the code used for the model.}
#' \item{n}{The number of observations used in the fit.}
#' \item{groups}{The maximum number of groups used in the fit.}
#' \item{n.events}{The number of events observed in the fit.}
#' \item{n.eventsbygrp}{A vector of length the number of groups
#' giving the number of observed events in each group.}
#' \item{loglik}{The marginal log-likelihood in the parametric case.}
#' \item{loglikPenal}{The marginal penalized log-likelihood
#' in the semiparametric case.}
#' \item{coef}{The regression coefficients.}
#' \item{varH}{The variance matrix of
#' the regression coefficients before positivity constraint transformation.
#' Then, the delta method is needed to obtain the estimated variance parameters.
#' That is why some variances don't match with
#' the printed values at the end of the model.}
#' \item{varHtotal}{The variance matrix of
#' all the parameters before positivity constraint transformation.
#' Then, the delta method is needed to obtain the estimated variance parameters.
#' That is why some variances don't match with
#' the printed values at the end of the model.}
#' \item{varHIH}{The robust estimation of
#' the variance matrix of the regression coefficients}
#' \item{varHIHtotal}{The robust estimation of
#' the variance matrix of all parameters.}
#' \item{x}{Matrix of times where the hazard functions are estimated.}
#By default, seq(0,max(time),length=99),
#where time is the vector of survival times.}
#' \item{xSu}{Matrix of times where the survival functions are estimated.}
#By default, seq(0,max(time),length=99),
#where time is the vector of survival times.}
#' \item{lam}{Array (dim=3) of baseline hazard estimates
#' and confidence bands.}
#' \item{surv}{Array (dim=3) of baseline survival estimates
#' and confidence bands.}
#' \item{type}{Character string specifying the type of censoring,
#' see the \code{Surv} function for more details.}
#' \item{n.strat}{Number of strata.}
#' \item{n.iter}{Number of iterations needed to converge.}
#' \item{median}{The value of the median survival and its confidence bands.
#' If there are two strata or more, the first value corresponds to the value
#' for the first strata, etc.}
#' \item{LCV}{The approximated likelihood cross-validation criterion in the
#' semiparametric case.
#' With H (resp. H\out{<sub>pen</sub>}) the hessian matrix
#' of log-likelihood (resp. penalized log-likelihood),
#' EDF = H\out{<sub>pen</sub>}\out{<sup>-1</sup>} H
#' the effective degrees of freedom,
#' L(\eqn{\xi},\eqn{\theta}) the log-likelihood and
#' n the number of observations,
#' \deqn{LCV = 1/n x (trace(EDF) - L(\xi,\theta)).}}
#' \item{AIC}{The Akaike information Criterion for the parametric case.
#' With p the number of parameters,
#' n the number of observations and L(\eqn{\xi},\eqn{\theta}) the log-likelihood,
#' \deqn{AIC = 1/n x (p - L(\xi,\theta)).}}
#' \item{npar}{Number of parameters.}
#' \item{nvar}{Number of explanatory variables.}
#' \item{typeof}{Indicator of the type of hazard functions computed :
#' 0 for "Splines", 2 for "parametric".}
#' \item{istop}{Convergence indicator:
#' 1 if convergence is reached,
#' 2 if convergence is not reached,
#' 3 if the hessian matrix is not positive definite,
#' 4 if a numerical problem has occurred in the likelihood calculation}
#' \item{shape.param}{Shape parameter for the parametric hazard function
#' (a Weibull distribution is used for proportional and additive hazards models,
#' a log-logistic distribution is used for proportional odds models,
#' a log-normal distribution is used for probit models).}
#' \item{scale.param}{Scale parameter for the parametric hazard function.}
#' \item{Names.data}{Name of the dataset.}
#' \item{Frailty}{Logical value. Was model with frailties fitted ?}
#\item{martingale.res}{martingale residuals for each cluster.}
#\item{martingaleCox}{martingale residuals for observation in the Cox model.}
#\item{frailty.pred}{empirical Bayes prediction of the frailty term
#(ie, using conditional posterior distributions).}
#\item{frailty.var}{variance of the empirical Bayes prediction of the frailty
#term (only for gamma frailty models).}
#\item{frailty.sd}{standard error of the frailty empirical Bayes prediction
#(only for gamma frailty models).}
#' \item{linear.pred}{Linear predictor:
#' \bold{\eqn{\beta}}'\bold{\eqn{X}} in the generalized survival models or
#' \bold{\eqn{\beta}}'\bold{\eqn{X}} + log(\eqn{u}\out{<sub>i</sub>})
#' in the shared frailty generalized survival models.}
#otherwise uses "Beta'X + w_i" for log-normal frailty distribution.}
#' \item{BetaTpsMat}{Matrix of time varying-effects and confidence bands
#' (the first column used for abscissa of times).}
#' \item{nvartimedep}{Number of covariates with time-varying effects.}
#' \item{Names.vardep}{Name of the covariates with time-varying effects.}
#' \item{EPS}{Convergence criteria concerning
#' the parameters, the likelihood and the gradient.}
#' \item{family}{Type of Generalized Survival Model fitted
#' (0 for PH, 1 for PO, 2 for probit, 3 for AH).}
#' \item{global_chisq.test}{A binary variable equals to 0 when no multivariate
#' Wald is given, 1 otherwise.}
#' \item{beta_p.value}{p-values of the Wald test for the estimated
#' regression coefficients.}
#' \item{cross.Val}{Logical value. Is cross validation procedure used for
#' estimating the smoothing parameters in the penalized likelihood estimation?}
#' \item{DoF}{Degrees of freedom associated with the smoothing parameter
#' \code{kappa}.}
#' \item{kappa}{A vector with the smoothing parameters in the penalized
#' likelihood estimation corresponding to each baseline function as components.}
#' \item{n.knots}{Number of knots for estimating the baseline functions in the
#' penalized likelihood estimation.}
#\item{nbintervR}{Number of intervals (between 1 and 20) for the
#parametric hazard functions ("Piecewise-per", "Piecewise-equi").}
#' \item{n.knots.temp}{Initial value for the number of knots.}
#' \item{global_chisq}{A vector with the values of each multivariate Wald test.}
#' \item{dof_chisq}{A vector with the degree of freedom for each multivariate
#' Wald test.}
#' \item{p.global_chisq}{A vector with the p-values for each global multivariate
#' Wald test.}
#' \item{names.factor}{Names of the "as.factor" variables.}
#' \item{Xlevels}{Vector of the values that factor might have taken.}
#\item{contrasts}{Type of contrast for factor variable.}
#'
#'
#' The following components are specific to \bold{shared} models.
#'
#' \item{equidistant}{Indicator for the intervals used in the spline estimation
#' of baseline hazard functions :
#' 1 for equidistant intervals ; 0 for intervals using percentile
#' (note: \code{equidistant = 2} in case of parametric estimation).}
#\item{intcens}{Logical value. Indicator if a joint frailty
#model with interval-censored data was fitted)}
#' \item{Names.cluster}{Cluster names.}
#' \item{theta}{Variance of the gamma frailty parameter, i.e.
#' Var(\eqn{u}\out{<sub>i</sub>}).}
#' \item{varTheta}{Variance of parameter \code{theta}.}
#' \item{theta_p.value}{p-value of the Wald test for
#' the estimated variance of the gamma frailty.}
#\item{sigma2}{Variance of the log-normal frailty parameter
#\eqn{(\bold{Var}(\eta_i))}.}
#\item{sigma2_p.value}{p-value of the Wald test for the estimated variance of
#the log-normal frailty.}
#'
#'
#' The following components are specific to \bold{joint} models.
#'
#\item{intcens}{Logical value. Indicator if a joint frailty model with
#interval-censored data was fitted)}
#' \item{formula}{The formula part of the code
#' used for the recurrent events.}
#' \item{formula.terminalEvent}{The formula part of the code
#' used for the terminal model.}
#' \item{n.deaths}{Number of observed deaths.}
#' \item{n.censored}{Number of censored individuals.}
#' \item{theta}{Variance of the gamma frailty parameter, i.e.
#' Var(\eqn{u}\out{<sub>i</sub>}).}
#\item{sigma2}{Variance of the log-normal frailty parameter
#\eqn{(\bold{Var}(\eta_i))} or \eqn{(\bold{Var}(v_i))}.}
#\item{eta}{Variance of the second gamma frailty parameter in general joint frailty models
#\eqn{(\bold{Var}(v_i))}.}
#' \item{indic_alpha}{Indicator if a joint frailty model with
#' \eqn{\alpha} parameter was fitted.}
#' \item{alpha}{The coefficient \eqn{\alpha} associated
#' with the frailty parameter in the terminal hazard function.}
#\item{nbintervR}{Number of intervals (between 1 and 20) for the
#recurrent parametric hazard functions ("Piecewise-per", "Piecewise-equi").}
#\item{nbintervDC}{Number of intervals (between 1 and 20) for the death
#parametric hazard functions ("Piecewise-per", "Piecewise-equi").}
#' \item{nvar}{A vector with the number of covariates
#' of each type of hazard function as components.}
#' \item{nvarnotdep}{A vector with the number of constant effect covariates
#' of each type of hazard function as components.}
#' \item{nvarRec}{Number of recurrent explanatory variables.}
#' \item{nvarEnd}{Number of death explanatory variables.}
#' \item{noVar1}{Indicator of recurrent explanatory variables.}
#' \item{noVar2}{Indicator of death explanatory variables.}
#' \item{Names.vardep}{Name of the covariates with time-varying effects
#' for the recurrent events.}
#' \item{Names.vardepdc}{Name of the covariates with time-varying effects
#' for the terminal event.}
#' \item{xR}{Matrix of times where both survival and hazard function
#' are estimated for the recurrent event.}
#By default seq(0,max(time),length=99),
#where time is the vector of survival times.}
#' \item{xD}{Matrix of times for the terminal event.}
#' \item{lamR}{Array (dim=3) of hazard estimates and confidence bands
#' for recurrent event.}
#' \item{lamD}{The same value as \code{lamR} for the terminal event.}
#' \item{survR}{Array (dim=3) of baseline survival estimates and
#' confidence bands for recurrent event.}
#' \item{survD}{The same value as \code{survR} for the terminal event.}
#' \item{nb.gh}{Number of nodes for the Gaussian-Hermite quadrature.}
#' \item{nb.gl}{Number of nodes for the Gaussian-Laguerre quadrature.}
#' \item{medianR}{The value of the median survival for the recurrent events
#' and its confidence bands.}
#' \item{medianD}{The value of the median survival for the terminal event
#' and its confidence bands.}
#' \item{names.factor}{Names of the "as.factor" variables
#' for the recurrent events.}
#' \item{names.factordc}{Names of the "as.factor" variables
#' for the terminal event.}
#' \item{Xlevels}{Vector of the values that factor might have taken
#' for the recurrent events.}
#' \item{Xlevels2}{Vector of the values that factor might have taken
#' for the terminal event.}
#\item{martingale.res}{martingale residuals for each cluster (recurrent).}
#\item{martingaledeath.res}{martingale residuals for
#each cluster (death).}
#' \item{linear.pred}{Linear predictor for the recurrent part:
#' \bold{\eqn{\beta}}'\bold{\eqn{X}} + log(\eqn{u}\out{<sub>i</sub>}).}
#' \item{lineardeath.pred}{Linear predictor for the terminal part:
#' \bold{\eqn{\beta}}'\bold{\eqn{X}} + \eqn{\alpha} x log(\eqn{u}\out{<sub>i</sub>}).}
#' \item{Xlevels}{Vector of the values that factor might have taken
#' for the recurrent part.}
#\item{contrasts}{type of contrast for factor variable
#for the recurrent part.}
#' \item{Xlevels2}{vector of the values that factor might have taken
#' for the death part.}
#\item{contrasts2}{type of contrast for factor variable
#for the death part.}
#' \item{BetaTpsMat}{Matrix of time varying-effects and confidence bands for
#' recurrent event (the first column used for abscissa of times of recurrence).}
#' \item{BetaTpsMatDc}{Matrix of time varying-effects and confidence bands for
#' terminal event (the first column used for abscissa of times of death).}
#' \item{alpha_p.value}{p-value of the Wald test for the estimated \eqn{\alpha}.}
#\item{ncc}{Logical value whether nested case-control design
#with weights was used for the joint model.}
#'
#'
#The following components are specific to \bold{nested} models.
#\item{alpha}{variance of the cluster effect \eqn{(\bold{Var}(v_{i}))}}
#\item{eta}{variance of the subcluster effect \eqn{(\bold{Var}(w_{ij}))}}
#\item{subgroups}{the maximum number of subgroups used in the fit.}
#\item{frailty.pred.group}{empirical Bayes prediction of the frailty term by
#group.}
#\item{frailty.pred.subgroup}{empirical Bayes prediction of the
#frailty term by subgroup.}
#\item{linear.pred}{linear predictor: uses "Beta'X + log v_i.w_ij".}
#\item{subgbyg}{subgroup by group.}
#\item{n.strat}{A vector with the number of covariates of each type of hazard
#function as components.}
#\item{alpha_p.value}{p-value of the Wald test for the estimated
#variance of the cluster effect.}
#\item{eta_p.value}{p-value of the Wald test
#for the estimated variance of the subcluster effect.}
#'
#'
#The following components are specific to \bold{joint nested frailty} models.
#\item{theta}{variance of the subcluster effect \eqn{(\bold{Var}(u_{fi}))}}
#\item{eta}{variance of the cluster effect \eqn{(\bold{Var}(\omega_f))}}
#\item{alpha}{the power coefficient \eqn{\alpha} associated with the frailty
#parameter (\eqn{u_{fi}}) in the terminal event hazard function.}
#\item{ksi}{the power coefficient \eqn{\xi} associated with the frailty
#parameter (\eqn{\omega_f}) in the recurrent event hazard function.}
#\item{indic_alpha}{indicator if a joint frailty model with \eqn{\alpha}
#parameter was fitted or not.}
#\item{indic_ksi}{indicator if a joint frailty
#model with \eqn{\xi} parameter was fitted or not.}
#\item{frailty.fam.pred}{empirical Bayes prediction of the frailty term by
#family.} \item{eta_p.value}{p-value of the Wald test for the estimated
#variance of the cluster effect.}
#\item{alpha_p.value}{p-value of the Wald
#test for the estimated power coefficient \eqn{\alpha}.}
#\item{ksi_p.value}{p-value of the Wald test for the estimated power
#coefficient \eqn{\xi}.}
#'
#'
#' @note
#From a prediction aim, we recommend you to input a data sorted by the
#group variable with numerical numbers from 1 to n (number of groups).
#In case of a nested model, we recommend you to input a data sorted by the
#group variable then sorted by the subgroup variable both with numerical
#numbers from 1 to n (number of groups) and from 1 to m (number or subgroups).
#' In the flexible semiparametric case, smoothing parameters \code{kappa} and
#' number of knots \code{n.knots} are the arguments that the user have to change
#' if the fitted model does not converge.
#' \code{n.knots} takes integer values between 4 and 20.
#' But with \code{n.knots=20}, the model would take a long time to converge.
#' So, usually, begin first with \code{n.knots=7}, and increase it step by step
#' until it converges.
#' \code{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.
#'
#'
#' @seealso
#' \code{\link[survival]{Surv}},
#' \code{\link{terminal}},
#' \code{\link{timedep}}
#'
#'
#' @references
#' J. Chauvet and V. Rondeau (2021). A flexible class of generalized
#' joint frailty models for the analysis of survival endpoints. In revision.
#'
#' Liu XR, Pawitan Y, Clements M. (2018)
#' Parametric and penalized generalized survival models.
#' \emph{Statistical Methods in Medical Research} \bold{27}(5), 1531-1546.
#'
#' Liu XR, Pawitan Y, Clements MS. (2017)
#' Generalized survival models for correlated time-to-event data.
#' \emph{Statistics in Medicine} \bold{36}(29), 4743-4762.
#'
#' A. Krol, A. Mauguen, Y. Mazroui, A. Laurent, S. Michiels and V. Rondeau
#' (2017). Tutorial in Joint Modeling and Prediction: A Statistical Software
#' for Correlated Longitudinal Outcomes, Recurrent Events and a Terminal Event.
#' \emph{Journal of Statistical Software} \bold{81}(3), 1-52.
#'
#' V. Rondeau, Y. Mazroui and J. R. Gonzalez (2012). Frailtypack: An R package
#' for the analysis of correlated survival data with frailty models using
#' penalized likelihood estimation or parametric estimation. \emph{Journal of
#' Statistical Software} \bold{47}, 1-28.
#'
#' V. Rondeau, J.P. Pignon, S. Michiels (2011). A joint model for the
#' dependance between clustered times to tumour progression and deaths: A
#' meta-analysis of chemotherapy in head and neck cancer. \emph{Statistical
#' methods in medical research} \bold{897}, 1-19.
#'
#' V. Rondeau, S. Mathoulin-Pellissier, H. Jacqmin-Gadda, V. Brouste, P.
#' Soubeyran (2007). Joint frailty models for recurring events and death using
#' maximum penalized likelihood estimation:application on cancer events.
#' \emph{Biostatistics} \bold{8},4, 708-721.
#'
#' V. Rondeau, D. Commenges, and P. Joly (2003). Maximum penalized likelihood
#' estimation in a gamma-frailty model. \emph{Lifetime Data Analysis} \bold{9},
#' 139-153.
#'
#' C.A. McGilchrist, and C.W. Aisbett (1991). Regression with frailty in
#' survival analysis. \emph{Biometrics} \bold{47}, 461-466.
#'
#' D. Marquardt (1963). An algorithm for least-squares estimation of nonlinear
#' parameters. \emph{SIAM Journal of Applied Mathematics}, 431-441.
#'
#'
#' @keywords models
#' @export
#'
#'
#' @examples
#' \donttest{
#'
#' #############################################################################
#' # ----- GENERALIZED SURVIVAL MODELS (without frailties) ----- #
#' #############################################################################
#'
#' adult.retino = retinopathy[retinopathy$type == "adult", ]
#' adult.retino[adult.retino$futime >= 50, "status"] = 0
#' adult.retino[adult.retino$futime >= 50, "futime"] = 50
#'
#' ### --- Parametric PH, AH, PO and probit models --- ###
#'
#' GenfrailtyPenal(formula=Surv(futime,status)~trt, data=adult.retino,
#' hazard="parametric", family="PH")
#' GenfrailtyPenal(formula=Surv(futime,status)~trt, data=adult.retino,
#' hazard="parametric", family="AH")
#' GenfrailtyPenal(formula=Surv(futime,status)~trt, data=adult.retino,
#' hazard="parametric", family="PO")
#' GenfrailtyPenal(formula=Surv(futime,status)~trt, data=adult.retino,
#' hazard="parametric", family="probit")
#'
#' ### --- Semi-parametric PH and AH models --- ###
#'
#' GenfrailtyPenal(formula=Surv(futime,status)~timedep(trt), data=adult.retino,
#' family="PH", hazard="Splines", n.knots=8, kappa=10^6, betaknots=1, betaorder=2)
#' GenfrailtyPenal(formula=Surv(futime,status)~timedep(trt), data=adult.retino,
#' family="AH", hazard="Splines", n.knots=8, kappa=10^10, betaknots=1, betaorder=2)
#'
#'
#'
#' #############################################################################
#' # ----- SHARED FRAILTY GENERALIZED SURVIVAL MODELS ----- #
#' #############################################################################
#'
#' adult.retino = retinopathy[retinopathy$type == "adult", ]
#' adult.retino[adult.retino$futime >= 50, "status"] = 0
#' adult.retino[adult.retino$futime >= 50, "futime"] = 50
#'
#' ### --- Parametric PH, AH, PO and probit models --- ###
#'
#' GenfrailtyPenal(formula=Surv(futime,status)~trt+cluster(id), data=adult.retino,
#' hazard="parametric", family="PH")
#' GenfrailtyPenal(formula=Surv(futime,status)~trt+cluster(id), data=adult.retino,
#' hazard="parametric", family="AH")
#' GenfrailtyPenal(formula=Surv(futime,status)~trt+cluster(id), data=adult.retino,
#' hazard="parametric", family="PO")
#' GenfrailtyPenal(formula=Surv(futime,status)~trt+cluster(id), data=adult.retino,
#' hazard="parametric", family="probit")
#'
#' ### --- Semi-parametric PH and AH models --- ###
#'
#' GenfrailtyPenal(formula=Surv(futime,status)~cluster(id)+timedep(trt),
#' data=adult.retino, family="PH", hazard="Splines",
#' n.knots=8, kappa=10^6, betaknots=1, betaorder=2)
#' GenfrailtyPenal(formula=Surv(futime,status)~cluster(id)+timedep(trt),
#' data=adult.retino, family="AH", hazard="Splines",
#' n.knots=8, kappa=10^10, betaknots=1, betaorder=2)
#'
#'
#'
#' #############################################################################
#' # ----- JOINT FRAILTY GENERALIZED SURVIVAL MODELS ----- #
#' #############################################################################
#'
#' data("readmission")
#' readmission[, 3:5] = readmission[, 3:5]/365.25
#'
#' ### --- Parametric dual-PH, AH, PO and probit models --- ###
#'
#' GenfrailtyPenal(
#' formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+chemo,
#' formula.terminalEvent=~sex+dukes+chemo, data=readmission, recurrentAG=TRUE,
#' hazard="parametric", family=c("PH","PH"))
#' GenfrailtyPenal(
#' formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+chemo,
#' formula.terminalEvent=~sex+dukes+chemo, data=readmission, recurrentAG=TRUE,
#' hazard="parametric", family=c("AH","AH"))
#' GenfrailtyPenal(
#' formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+chemo,
#' formula.terminalEvent=~sex+dukes+chemo, data=readmission, recurrentAG=TRUE,
#' hazard="parametric", family=c("PO","PO"))
#' GenfrailtyPenal(
#' formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+chemo,
#' formula.terminalEvent=~sex+dukes+chemo, data=readmission, recurrentAG=TRUE,
#' hazard="parametric", family=c("probit","probit"))
#'
#' ### --- Semi-parametric dual-PH and AH models --- ###
#'
#' GenfrailtyPenal(
#' formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+timedep(chemo),
#' formula.terminalEvent=~sex+dukes+timedep(chemo), data=readmission, recurrentAG=TRUE,
#' hazard="Splines", family=c("PH","PH"),
#' n.knots=5, kappa=c(100,100), betaknots=1, betaorder=3)
#' GenfrailtyPenal(
#' formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+timedep(chemo),
#' formula.terminalEvent=~sex+dukes+timedep(chemo), data=readmission, recurrentAG=TRUE,
#' hazard="Splines", family=c("AH","AH"),
#' n.knots=5, kappa=c(600,600), betaknots=1, betaorder=3)
#'
#' }
#'
#'
"GenfrailtyPenal" <-
function (formula, formula.terminalEvent, data, recurrentAG=FALSE,
family, hazard="Splines", n.knots, kappa, betaknots=1,betaorder=3,
RandDist="Gamma", init.B, init.Theta, init.Alpha, Alpha, maxit=300,
nb.gh, nb.gl, LIMparam=1e-3, LIMlogl=1e-3, LIMderiv=1e-3,
print.times=TRUE,
cross.validation=FALSE, jointGeneral, nb.int, initialize=TRUE,
init.Ksi, Ksi, init.Eta){
# Ajout de la fonction minmin issue de print.survfit, permettant de calculer la mediane
minmin <- function(y, x) {
tolerance <- .Machine$double.eps^.5 #same as used in all.equal()
keep <- (!is.na(y) & y <(.5 + tolerance))
if (!any(keep)) NA
else {
x <- x[keep]
y <- y[keep]
if (abs(y[1]-.5) <tolerance && any(y< y[1]))
(x[1] + x[min(which(y<y[1]))])/2
else x[1]
}
}
# JOCELYN : options possibles pour 'family'
if (!all(family %in% c("PH","AH","PO","probit"))){
stop("Available options for argument 'family' are 'PH', 'AH', PO' and 'probit'")
}
# JOCELYN : pas de parametres "initialize", "init.Ksi", "Ksi", "init.Eta"
if (!missing(initialize) | !missing(init.Ksi) | !missing(Ksi) | !missing(init.Eta)){
stop ("Options 'initialize', 'init.Ksi', 'Ksi' and 'init.Eta' are only available for joint nested frailty models, which are not implemented in the generalized framework.")
}
# JOCELYN : Pas de jointGeneral (pour l'instant) dans GenfrailtyPenal
if (!missing(jointGeneral) & missing(formula.terminalEvent)) stop("You have to be in a JOINT frailty model to consider two independent frailty terms!")
if (!missing(jointGeneral)) stop("Considering two independent frailty terms is not yet an option in the framework of generalized joint frailty models.")
if (missing(jointGeneral)) jointGeneral<-FALSE
if (!missing(init.Eta) & jointGeneral) init.Alpha <- init.Eta
# al suppression de l'argument joint
if (!missing(formula.terminalEvent)) joint <- TRUE
else joint <- FALSE
if ((!missing(Alpha) | !missing(init.Alpha)) & !joint) stop("init.Alpha and Alpha parameters belong to joint frailty model")
#ad 15/02/12 :add Audrey
m2 <- match.call()
m2$formula <- m2$formula.terminalEvent <- m2$recurrentAG <- m2$cross.validation <- m2$n.knots <- m2$kappa <- m2$maxit <- m2$hazard <- m2$nb.int <- m2$RandDist <- m2$betaorder <- m2$betaknots <- m2$init.B <- m2$LIMparam <- m2$LIMlogl <- m2$LIMderiv <- m2$print.times <- m2$init.Theta <- m2$init.Alpha <- m2$Alpha <- m2$init.Ksi <- m2$Ksi <- m2$init.Eta <- m2$Eta <- m2$initialize <- m2$nb.gh <- m2$nb.gl <- m2$family <- m2$... <- NULL
Names.data <- m2$data
#### Betaknots et betaorder ####
if (betaknots > 10) stop("Number of knots for beta(t) greater than 10 is useless, please choose a number between 0 and 10 (3 is optimal)")
if ((betaorder == 0)|(betaorder > 4)) stop("B-splines order for beta(t) must be a number between 1 and 4 (3 is optimal)")
#### Frailty distribution specification ####
# JOCELYN : Uniquement "Gamma" pour l'instant
if (!RandDist=="Gamma") stop("For now, only 'Gamma' distribution for frailties is allowed")
if (!(RandDist %in% c("Gamma","LogN"))) { stop("Only 'Gamma' and 'LogN' distributions for frailties are allowed") }
logNormal <- switch(RandDist,"Gamma"=0,"LogN"=1)
# JOCELYN : hazard="parametric" or "Splines", et c'est tout !
if( !(hazard %in% c("parametric","Splines")) ){
stop("Only 'parametric' or 'Splines' hazard can be specified in hazard argument.")
}
if( (hazard=="Splines") & any(family %in% c("PO","probit")) ){
stop("Splines hazard is only available for PH and AH models")
}
if (RandDist=="LogN" & jointGeneral==TRUE) stop("Log normal distribution is not available for the Joint General Model !")
if ((hazard!="Splines") & jointGeneral== TRUE) stop("No general joint frailty model allowed here! Only 'Splines' is allowed")
##### hazard specification ######
haztemp <- hazard
hazard <- strsplit(hazard,split="-")
hazard <- unlist(hazard)
if(!(length(hazard) %in% c(1,2))){stop("Please check and revise the hazard argument according to the format specified in the help.")}
### longueur hazard = 1
if((all.equal(length(hazard),1)==T)==T){
if(!(hazard %in% c("parametric","Piecewise","Splines"))){
stop("Only 'parametric', 'Splines' or 'Piecewise' hazard can be specified in hazard argument.")
}else{
typeof <- switch(hazard,"Splines"=0,"Piecewise"=1,"parametric"=2)
### Splines (equidistant par defaut)
if (typeof == 0){
if (!missing(nb.int)){
stop("When the hazard function equals 'Splines' or 'parametric', 'nb.int' argument must be deleted.")
}
size1 <- 100
size2 <- 100
equidistant <- 1
nbintervR <- 0
nbintervDC <- 0
}
### Weibull
if (typeof == 2){
if (!missing(nb.int)){
stop("When the hazard function equals 'Splines' or 'parametric', 'nb.int' argument must be deleted.")
}
size1 <- 100
size2 <- 100
equidistant <- 2
nbintervR <- 0
nbintervDC <- 0
}
if (typeof == 1){
stop ("The hazard argument is incorrectly specified. Type of hazard are required ('per' or 'equi'). Please refer to the help file of frailtypack.")
}
}
}else{
#### longueur hazard > 1
if(all(!(c("Splines","Piecewise") %in% hazard))){
stop("Only 'Splines' and 'Piecewise' hazard can be specified in hazard argument in this case")
}
### Splines percentile
if ("Splines" %in% hazard){
typeof <- 0
if(!(all(hazard %in% c("Splines","per")))){
stop ("The hazard argument is incorrectly specified. Only 'per' is allowed with 'Splines'. Please refer to the help file of frailtypack.")
}else{
if (haztemp == "Splines-per") equidistant <- 0
else equidistant <- 1
size1 <- 100
size2 <- 100
nbintervR <- 0
nbintervDC <- 0
}
}
### Piecewise (per or equi)
if ("Piecewise" %in% hazard){
typeof <- 1
if(!(all(hazard %in% c("Piecewise","per","equi")))){
stop ("The hazard argument is incorrectly specified. Type of hazard are required ('per' or 'equi'). Please refer to the help file of frailtypack.")
}else{
if (!(haztemp %in% c("Piecewise-per","Piecewise-equi"))){
stop ("The hazard argument is incorrectly specified. Please refer to the help file of frailtypack.")
}
equidistant <- switch(haztemp,"Piecewise-per"=0,"Piecewise-equi"=1)
}
}
}
#ad Julien pour nb.gh et nb.gl
if (missing(nb.gh)) {
if (typeof == 0) {nb.gh <- 20}
else {nb.gh <- 32}
}
if (!(nb.gh %in% c(5,7,9,12,15,20,32))) stop("nb.gh must be chosen among 5,7,9,12,15,20 and 32")
if (missing(nb.gl)) {
if (typeof == 0) {nb.gl <- 20}
else {nb.gl <- 32}
}
if (!(nb.gl %in% c(20,32))) stop("nb.gl must be chosen among 20 and 32")
#AD:
if (missing(formula))stop("The argument formula must be specified in any model")
if (!inherits(formula, "formula"))stop("The argument formula must be a formula")
if(typeof == 0){
#AD:
if (missing(n.knots))stop("number of knots are required")
#AD:
n.knots.temp <- n.knots
#AD
if (n.knots<4) n.knots<-4
if (n.knots>20) n.knots<-20
if (missing(kappa))stop("smoothing parameter (kappa) is required")
#AD:
# JOCELYN : pas de validation croisee pour les modeles de survie generalises
if (cross.validation){
stop("The cross validation is not implemented for Generalized (Joint) Frailty Models")
}
if (length(kappa)>1 & cross.validation){
stop("The cross validation is not implemented for two strata or more")
}
if (joint & cross.validation){
stop("The cross validation is not implemented for the joint model")
}
}else{
if (!(missing(n.knots)) || !(missing(kappa)) || !(missing(cross.validation))){
stop("When parametric hazard function is specified, 'kappa', 'n.knots' and 'cross.validation' arguments must be deleted.")
}
n.knots <- 0
kappa <- 0
crossVal <- 0
}
call <- match.call()
# browser()
m <- match.call(expand.dots = FALSE) # recupere l'instruction de l'utilisateur
m$formula.terminalEvent <- m$n.knots <- m$recurrentAG <- m$cross.validation <- m$jointGeneral <- m$kappa <- m$maxit <- m$hazard <- m$nb.int <- m$RandDist <- m$betaorder <- m$betaknots <- m$init.B <- m$LIMparam <- m$LIMlogl <- m$LIMderiv <- m$print.times <- m$init.Theta <- m$init.Alpha <- m$Alpha <- m$init.Ksi <- m$Ksi <- m$init.Eta <- m$Eta <- m$initialize <- m$nb.gh <- m$nb.gl <- m$family <- m$... <- NULL
special <- c("strata", "cluster", "subcluster", "terminal","num.id","timedep", "wts") #wts for weights (ncc design) ncc - nested case-control
Terms <- if (missing(data)){
terms(formula, special)
}else{
terms(formula, special, data = data)
}
ord <- attr(Terms, "order") # longueur de ord=nbre de var.expli
#if (length(ord) & any(ord != 1))stop("Interaction terms are not valid for this function")
#si pas vide tous si il ya au moins un qui vaut 1 on arrete
m$formula <- Terms
m[[1]] <- as.name("model.frame") # m[[1]]=GenfrailtyPenal, il le remplace par model.frame en fait
m <- eval(m, sys.parent()) #ici la classe de m est un data.frame donc il recupere ce qu'on lui donne en argument
cluster <- attr(Terms, "specials")$cluster # (indice) nbre de var qui sont en fonction de cluster()
# al 13/02/14 : suppression de l'argument Frailty
if (length(cluster)) Frailty <- TRUE
else Frailty <- FALSE
subcluster <- attr(Terms, "specials")$subcluster #nbre de var qui sont en fonction de subcluster()
wts <- attr(Terms, "specials")$wts #nbre de var qui sont en fonction de wts()
# booleen pour voir si l'objet Y est reponse avant tri des donnees Surv ou SurvIC
classofY <- attr(model.extract(m, "response"),"class")
# attention le package pec rajoute un element dans l'attribut "class" des objets de survie
if (length(classofY)>1) classofY <- classofY[2]
typeofY <- attr(model.extract(m, "response"),"type") # type de reponse : interval etc..
# # # We retrieve here the t0 name give in the 'data' dataframe (needed for checking if left truncated data)
# # # vart0 <- dimnames(m)[[2]][1]
# # # vart0 <- gsub('(.*\\()(.*)', '\\2', strsplit(vart0,",")[[1]][1])
#Al : tri du jeu de donnees par cluster croissant
if (length(cluster)){
tempc <- untangle.specials(Terms, "cluster", 1:10)
ord <- attr(Terms, "order")[tempc$terms] # le type ordinal de la variable (0,1 ou 2)
if (any(ord > 1))stop("Cluster can not be used in an interaction")
m <- m[order(m[,tempc$vars]),] # soit que des nombres, soit des caracteres
ordre <- as.integer(row.names(m)) # recupere l'ordre du data set
cluster <- strata(m[, tempc$vars], shortlabel = TRUE)
uni.cluster <- unique(cluster)
}
#Verification de l'utilisation du parametre Ksi pour le nested joint modele
if ((!missing(Ksi) | !missing (init.Ksi)) & (!joint & !length(subcluster))) stop("init.Ksi and Ksi parameters belong only to nested joint frailty model")
# verification de la structure nested si besoin
if (length(subcluster) && Frailty == TRUE){
tempsub <- untangle.specials(Terms, "subcluster", 1:10)
ordsub <- attr(Terms, "order")[tempsub$terms] # type ordinal de la variable
if (any(ordsub > 1))stop("subcluster can not be used in an interaction")
if (any(ifelse(apply(ifelse(table(m[,tempsub$vars],m[,tempc$vars])>0,1,0),1,sum)==1,FALSE,TRUE))){
stop("nested structure is necessary to fit a nested model")
}
# tri par ordre croissant de subcluster a l'interieur des clusters
m <- m[order(m[,tempc$vars],m[,tempsub$vars]),]
ordretmp <- order(m[,tempsub$vars])
subcluster <- strata(m[, tempsub$vars], shortlabel = TRUE)
ordre <- as.integer(row.names(m))
subcluster <- as.integer(subcluster) # a determiner si il y en a besoin
curr <- subcluster[1]
subcluster[1] <- 1
for (i in 2:length(subcluster)) {
if (subcluster[i] == curr) { subcluster[i] <- subcluster[i-1] }
else {
curr <- subcluster[i]
subcluster[i] <- subcluster[i-1] + 1
}
}
}
#Al
if (NROW(m) == 0)stop("No (non-missing) observations") #nombre ligne different de 0
Y <- model.extract(m, "response") # objet de type Surv =Time
if (classofY == "SurvIC") intcens <- TRUE # booleen censure par intervalle
else intcens <- FALSE
if (intcens == TRUE) {
if (classofY != "SurvIC") stop("When interval censoring, must use the SurvIC fonction")
} else {
#if (!inherits(Y, "Surv")) stop("Response must be a survival object") #test si c bien un objet de type "Surv"
if (classofY != "Surv") stop("Response must be a survival object")
}
ll <- attr(Terms, "term.labels")#liste des variables explicatives
##########################################################################################################################
# add Myriam 11/05/2016 Checking left-truncating data
# # # troncat <- function(cluster, time0){
# # # mini = time0[1]
# # # indiv <- cluster[1]
# # # for (i in 1:length(time0)){
# # # if (cluster[i] != indiv){
# # # if (mini != 0) stop('Sorry but left-troncature is not allowed for joint or joint nested frailty models')
# # # else {
# # # indiv <- cluster[i]
# # # mini <- time0[i]
# # # }
# # # }
# # # else{
# # # if (mini > time0[i]) mini <- time0[i]
# # # }
# # # }
# # # }
# # # Timet0 <- data[vart0]
# # # if(length(cluster)){
# # # varclust <- gsub('(.*\\()(.*)\\)', '\\2', tempc$vars) # pour recuperer le nom de la variable definie par cluster()
# # # }
#verification pour le modele joint nested
# # # if (length(subcluster) & joint){
# # # varsubclust <- gsub('(.*\\()(.*)\\)', '\\2', tempsub$vars)
# # # Timet0 <- Timet0[order(data[,varclust], data[,varsubclust]),]
# # # if (typeofY == "counting") troncat(as.numeric(as.vector(data[order(data[,varclust], data[,varsubclust]),varsubclust])), Timet0)
# # # }
# verification pour le modele joint
# # # if (!length(subcluster) & joint){
# # # Timet0 <- list(Timet0)[[1]][order(data[,varclust]),]
# # # #data[,varclust] = data[order(data[,varclust]), varclust]
# # # if (typeofY == "counting") troncat(as.numeric(as.vector(data[order(data[,varclust]), varclust])), Timet0)
# # # }
# # # NOTE : 07-03-17 : this operation makes errors in GenfrailtyPenal execution when the user put in surv() formula
# # # covariates with a such name : 'data$var' or 'data[,2]' or 'data[,"var"]'
# # # So, for the moment we disable this process, any left-truncated data will not be detected by frailtypack
# # # and calculation will done with bias because left truncating is not took into account joint and joint nested modelling
##############################################################################################################################
#=========================================================>
mt <- attr(m, "terms") #m devient de class "formula" et "terms"
X <- if (!is.empty.model(mt)) model.matrix(mt, m) #, contrasts) #idem que mt sauf que ici les factor sont divise en plusieurs variables
ind.place <- unique(attr(X,"assign")[duplicated(attr(X,"assign"))]) ### unique : changement au 25/09/2014
vec.factor <- NULL
vec.factor <- c(vec.factor,ll[ind.place])
#=========================================================>
# On determine le nombre de categorie pour chaque var categorielle
strats <- attr(Terms, "specials")$strata #nbre de var qui sont en fonction de strata()
cluster <- attr(Terms, "specials")$cluster #nbre de var qui sont en fonction de cluster()
num.id <- attr(Terms, "specials")$num.id #nbre de var qui sont en fonction de patkey()
vartimedep <- attr(Terms, "specials")$timedep #nbre de var en fonction de timedep()
wts <- attr(Terms, "specials")$wts #nbre de var en fonction de wts()....redundant but keeping with cluster()
#booleen pour savoir si au moins une var depend du tps
if (is.null(vartimedep)) timedep <- 0
else timedep <- 1
# JOCELYN : Dans quels cas peut-on ajuster un modele avec certains effets dependant du temps ?
if( (missing(formula.terminalEvent)) & (hazard=="parametric") & (timedep==1)){
if (!(family=="PH")){
stop("For AH, PO and probit models, considering time-varying effects is not allowed when hazard='parametric'")
}
}
if( (missing(formula.terminalEvent)) & (hazard=="Splines") & (timedep==0)){
if (!(family=="PH")){
stop("For AH models, considering no time-varying effects is not allowed when hazard='Splines'")
}
}
if( (!missing(formula.terminalEvent)) & (hazard=="parametric") & (timedep==1)){
if (!all(family==c("PH","PH"))){
stop("For AH, PO and probit models, considering time-varying effects is not allowed when hazard='parametric'")
}
}
if( !(missing(formula.terminalEvent)) & (hazard=="Splines") & (timedep==0)){
if (!all(family==c("PH","PH"))){
stop("For AH models, considering no time-varying effects is allowed when hazard='Splines', but only using 'timedep(covariate)' with options 'betaknots=0' and 'betaorder=1'")
}
}
if (intcens & (equidistant == 0)) stop("You can not fit a model with a baseline hazard function estimated using percentiles and interval censoring")
if (intcens & timedep) stop("You can not use time-varying effect covariates with interval censoring")
if (intcens & cross.validation) stop("You can not do cross validation with interval censoring")
if (intcens & logNormal) stop("It is currently impossible to fit a model with interval censoring and a log normal distribution of the frailties")
if (timedep & logNormal) stop("You can not use time-varying effect covariates with a log normal distribution of the frailties")
if (intcens & joint & length(subcluster)) stop("Cox model doesn't support intervall-censored survival data")
if (intcens & !joint & length(subcluster)) stop ("Nested model doesn't support intervall-censored survival data")
if(is.null(num.id)){
joint.clust <- 1
}else{
joint.clust <- 0
if (!joint) stop("num.id function can only be used with joint models")
}
if (jointGeneral==TRUE) joint.clust <- 2
if (joint & length(subcluster)) joint.clust <- 3
subcluster <- attr(Terms, "specials")$subcluster #nbre de var qui sont en fonction de subcluster()
if (length(subcluster)){
ll <- ll[-grep("subcluster",ll)]
}
if (length(cluster)){
ll_tmp <- ll[grep("cluster",ll)]
ll <- ll[-grep("cluster",ll)]
pos1 <- grep("r",unlist(strsplit(ll_tmp,split="")))[1]+2
pos2 <- length(unlist(strsplit(ll_tmp,split="")))-1
Names.cluster <- substr(ll_tmp,start=pos1,stop=pos2) # nom du cluster
}
if (length(strats)){
ll <- ll[-grep("strata",ll)]
}
if (length(wts)){
ll <- ll[-grep("wts",ll)]
}
# plus besoin de as.factor() pour afficher le test de Wald global
if (length(grep("strata",vec.factor))) vec.factor <- vec.factor[-grep("strata",vec.factor)]
if (length(grep("cluster",vec.factor))) vec.factor <- vec.factor[-grep("cluster",vec.factor)]
if (length(grep("subcluster",vec.factor))) vec.factor <- vec.factor[-grep("subcluster",vec.factor)]
if (length(grep("num.id",vec.factor))) vec.factor <- vec.factor[-grep("num.id",vec.factor)]
if (length(grep("wts",vec.factor))) vec.factor <- vec.factor[-grep("wts",vec.factor)]
mat.factor <- matrix(vec.factor,ncol=1,nrow=length(vec.factor))
# Fonction servant a prendre les termes entre "as.factor" et (AK 04/11/2015) interactions
vec.factor <-apply(mat.factor,MARGIN=1,FUN=function(x){
if (length(grep("factor",x))>0){
if(length(grep(":",x))>0){
if(grep('\\(',unlist(strsplit(x,split="")))[1]<grep(":",unlist(strsplit(x,split="")))[1] && length(grep('\\(',unlist(strsplit(x,split=""))))==1){
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- grep(":",unlist(strsplit(x,split="")))[1]-2
pos3 <- grep(":",unlist(strsplit(x,split="")))[1]
pos4 <- length(unlist(strsplit(x,split="")))
return(paste(substr(x,start=pos1,stop=pos2),substr(x,start=pos3,stop=pos4),sep=""))
}else if(grep("\\(",unlist(strsplit(x,split="")))[1]>grep(":",unlist(strsplit(x,split="")))[1] && length(grep('\\(',unlist(strsplit(x,split=""))))==1){
pos2 <- grep(":",unlist(strsplit(x,split="")))[1]
pos3 <- grep("\\(",unlist(strsplit(x,split="")))[1]+1
pos4 <- length(unlist(strsplit(x,split="")))-1
return(paste(substr(x,start=1,stop=pos2),substr(x,start=pos3,stop=pos4),sep=""))
}else{#both factors
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- grep(":",unlist(strsplit(x,split="")))[1]-2
pos3 <- grep("\\(",unlist(strsplit(x,split="")))[2]+1
pos4 <- length(unlist(strsplit(x,split="")))-1
return(paste(substr(x,start=pos1,stop=pos2),":",substr(x,start=pos3,stop=pos4),sep=""))
}
}else{
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- length(unlist(strsplit(x,split="")))-1
return(substr(x,start=pos1,stop=pos2))}
}else{
return(x)
}
} )
if(length(grep("terminal",ll))>0){
ind.place <- grep(paste(vec.factor,collapse="|"),ll[-grep("terminal",ll)])
}
else{
ind.place <- grep(paste(vec.factor,collapse="|"),ll)
}
if(length(vec.factor) > 0){
vect.fact <- attr(X,"dimnames")[[2]]
#vect.fact <- vect.fact[grep("factor",vect.fact)]
vect.fact <- vect.fact[grep(paste(vec.factor,collapse="|"),vect.fact)]
# vect.fact <-apply(matrix(vect.fact,ncol=1,nrow=length(vect.fact)),MARGIN=1,FUN=function(x){
# pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
# pos2 <- grep(")",unlist(strsplit(x,split="")))[1]-1
# return(substr(x,start=pos1,stop=pos2))})
occur <- rep(0,length(vec.factor))
interaction<-as.vector(apply(matrix(vect.fact,nrow=length(vect.fact)),MARGIN=1,FUN=function(x){length(grep(":",unlist(strsplit(x,split=""))))}))
which.interaction <- which(interaction==1)
for(i in 1:length(vec.factor)){
if(length(grep(":",unlist(strsplit(vec.factor[i],split=""))))>0){
pos <- grep(":",unlist(strsplit(vec.factor[i],split="")))
length.grep <- 0
for(j in 1:length(vect.fact)){
if(j%in%which.interaction){
if(length(grep(substr(vec.factor[i],start=1,stop=pos-1),vect.fact[j]))>0 && length(grep(substr(vec.factor[i],start=pos+1,stop=length(unlist(strsplit(vec.factor[i],split="")))),vect.fact[j]))>0){
length.grep <- length.grep + 1
which <- i
}
}
}
occur[i] <- length.grep
}else{
if(length(vect.fact[-which.interaction])>0){occur[i] <- length(grep(vec.factor[i],vect.fact[-which.interaction]))
}else{occur[i] <- length(grep(vec.factor[i],vect.fact))}
}
}
}
if (joint & length(subcluster)){
if(initialize) initialize <- 1
else initialize <- 0
}
else{
if (!missing (initialize)) stop ('Sorry but \'initialize\' option is only available for joint nested frailty model')
}
#=========================================================>
terminalEvent <- attr(Terms, "specials")$terminal #nbre de var qui sont en fonction de terminal()
dropx <- NULL
if (length(cluster) & Frailty == TRUE){
tempc <- untangle.specials(Terms, "cluster", 1:10)
ord <- attr(Terms, "order")[tempc$terms]
if (any(ord > 1))stop("Cluster can not be used in an interaction")
cluster <- strata(m[, tempc$vars], shortlabel = TRUE)
dropx <- tempc$terms
uni.cluster<-unique(cluster)
}
else if (!length(cluster) & Frailty == TRUE){
stop("grouping variable is needed")
}
else if (length(cluster) & Frailty == FALSE){
stop("cluster not necessary for proportional hazard model")
}
else if (!length(cluster) & Frailty == FALSE){
cluster <- 1:nrow(m) #nrow(data) # valeurs inutiles pour un modele de Cox
uni.cluster <- 1:nrow(m) #nrow(data)
}
if (!missing(RandDist) & (Frailty == FALSE)){
stop("RandDist not necessary for proportional hazard model")
}
if (length(num.id)){
temppat <- untangle.specials(Terms, "num.id", 1:10)
num.id <- m[,temppat$vars]
dropx <- c(dropx,temppat$terms)
}
if(length(uni.cluster)==1){
stop("grouping variable must have more than 1 level")
}
if (length(subcluster)){
tempsub <- untangle.specials(Terms, "subcluster", 1:10)
ordsub <- attr(Terms, "order")[tempsub$terms]
if (any(ordsub > 1))stop("subcluster can not be used in an interaction")
subcluster <- strata(m[, tempsub$vars], shortlabel = TRUE)
dropx <- c(dropx,tempsub$terms)
uni.subcluster<-unique(subcluster)
#if (joint)stop("joint model is not implemented for nested model")
if(length(uni.subcluster)==1){
stop("subcluster variable must have more than 1 level")
}
}
#AD:
if (length(cluster) == length(subcluster)){
if (all(all.equal(cluster,subcluster)==T)){
stop("'Subgroup' variable and 'group' variable need to be different")
}
}
#AD:
if (length(strats)){
temp <- untangle.specials(Terms, "strata", 1)
dropx <- c(dropx, temp$terms)
if (length(temp$vars) == 1)strata.keep <- m[[temp$vars]]
else strata.keep <- strata(m[, temp$vars], shortlabel = TRUE)
strats <- as.numeric(strata.keep)
uni.strat<-length(unique(strats))
if (missing(kappa)) stop("smoothing parameter (kappa) is required")
if (length(subcluster)){
if (uni.strat > 2) {
if (joint) stop("maximum number of strata for nested frailty joint model is 2")
else stop("maximum number of strata for nested frailty joint model is 2")
}
}
else{
if (uni.strat > 6) stop("maximum number of strata is 6")
}
if ((uni.strat > 2) & (intcens)) stop("maximum number of strata for interval censored data is 2")
if ((uni.strat > 2) & (timedep)) stop("maximum number of strata for time-varying effect of covariates is 2")
}else{
uni.strat<-1
strats <- rep(1,nrow(data))
}
if (!joint & (typeof==0) & (length(kappa)!=uni.strat)) stop("wrong length of argument 'kappa' for the current stratification")
#AD: indicator of terminal()
ind.terminal <- length(terminalEvent)
#AD:
if (length(terminalEvent)){
tempterm <- untangle.specials(Terms, "terminal", 1:10)
#ici on comme terme tempterm$vars qui est le nom dans l'appel(ex;"terminal(death)"
#et tempterm$terms qui est la position de la variable dans l'appel, ici elle vient a la position 6
ord <- attr(Terms, "order")[tempterm$terms] # ord[6]=1 ici dans notre exemple
if (any(ord > 1))stop("Terminal can not be used in an interaction")
dropx <- c(dropx,tempterm$terms) # vecteur de position
terminal <- strata(m[, tempterm$vars], shortlabel = TRUE)
terminal <- as.numeric(as.character(terminal))
}
#IJ 2018:
if (length(wts)){
tempwts <- untangle.specials(Terms, "wts", 1:10)
ord <- attr(Terms, "order")[tempwts$terms] # ord[6]=1 ici dans notre exemple
if (any(ord > 1))stop("Weights can not be used in an interaction")
dropx <- c(dropx,tempwts$terms) # vecteur de position
weights.vec <- strata(m[, tempwts$vars], shortlabel = TRUE)
weights.vec <- as.numeric(as.character(weights.vec))
} else{
weights.vec <- rep(1, nrow(data))
}
weights.agg <- aggregate(weights.vec, by=list(cluster), FUN=function(x) x[length(x)])[,2]
#type <- attr(Y, "type")
type <- typeofY
if (type != "right" && type != "counting" && type != "interval" && type != "intervaltronc") { # Cox supporte desormais la censure par intervalle
stop(paste("Cox model doesn't support \"", type, "\" survival data", sep = ""))
}
# if ((type == "interval" || type == "interval2" || type == "intervaltronc") && intcens == FALSE) { # rajout
# stop("You are trying to do interval censoring without intcens = TRUE")
# }
if (type != "counting" && recurrentAG) {
stop("recurrentAG needs counting process formulation")
}
if (intcens == TRUE & recurrentAG == TRUE) {
stop("recurrentAG invalid for interval censored data")
}
#drop contient les position liees au fonction() ic ex:cluster(id) et terminal(death)
if (length(dropx)){
newTerms <- Terms[-dropx]
}else{
newTerms <- Terms
}
#newTerm vaut Terms - les variables dont les position sont dans drop
X <- model.matrix(newTerms, m)
assign <- lapply(attrassign(X, newTerms)[-1], function(x) x - 1)
Xlevels <- .getXlevels(newTerms, m)
contr.save <- attr(X, 'contrasts')
# assigne donne la position pour chaque variables
#ncol(X) : nombre de variable sans sans les fonction speciaux comme terminal()...+id
if(length(vec.factor) > 0){
#========================================>
position <- unlist(assign,use.names=F)
}
#========================================>
if (ncol(X) == 1){
X<-X-1
noVar1 <- 1
}else{
X <- X[, -1, drop = FALSE]
noVar1 <- 0
}
# on enleve ensuite la premiere colonne correspondant a id
nvar<-ncol(X) #nvar==1 correspond a 2 situations:
# au cas ou on a aucune var explicative dans la partie rec, mais X=0
# cas ou on a 1seul var explicative, ici X est en general different de 0
varnotdep <- colnames(X)[-grep("timedep",colnames(X))]
vardep <- colnames(X)[grep("timedep",colnames(X))]
vardep <- apply(matrix(vardep,ncol=1,nrow=length(vardep)),1,timedep.names)
if (length(intersect(varnotdep,vardep)) != 0) {
stop("A variable is both used as a constant and time-varying effect covariate")
}
nvartimedep <- length(vardep)
filtretps <- rep(0,nvar)
filtretps[grep("timedep",colnames(X))] <- 1
var<-matrix(c(X),nrow=nrow(X),ncol=nvar)
n<-nrow(X)
#add Alexandre 04/06/2012
#lire les donnees differemment si censure par intervalle
if (intcens==TRUE) {
if (type=="intervaltronc") {
tt0 <- Y[,1]
tt1 <- Y[,2]
ttU <- Y[,3]
cens <- Y[,4]
}
else {
tt0 <- rep(0,n)
tt1 <- Y[,1]
ttU <- Y[,2]
cens <- Y[,3]
tt1[tt1==0] <- 0.1
}
}
else {
if (type=="right"){
tt0 <- rep(0,n)
tt1 <- Y[,1]
cens <- Y[,2]
ttU <- Y[,1] # rajouter quand meme dans frailPenal mais ne sera pas utilise
}
else {
tt0 <- Y[,1]
tt1 <- Y[,2]
cens <- Y[,3]
ttU <- Y[,2] # ne sera pas pris en compte dans le tri des temps de survie dans frailtypack.f90
} # attention ne pas mettre de 0 sinon en cas de left trunc probleme dans la logV
}
if (min(cens)==0) cens.data<-1
if (min(cens)==1 && max(cens)==1) cens.data<-0
AG<-ifelse(recurrentAG,1,0)
if (typeof == 0){
crossVal<-ifelse(cross.validation,0,1)
}
#=======================================>
#======= Construction du vecteur des indicatrice
if(length(vec.factor) > 0){
# ind.place <- ind.place -1
k <- 0
for(i in 1:length(vec.factor)){
ind.place[i] <- ind.place[i]+k
k <- k + occur[i]-1
}
}
#==================================
# Begin SHARED MODEL
#
if (!joint & !length(subcluster)){
if((equidistant %in% c(0,1)) & (typeof == 1)){
if (missing(nb.int)) stop("Number of time interval 'nb.int' is required")
if (length(nb.int) != 1) stop("Wrong length of number of time interval argument 'nb.int'")
if (!inherits(nb.int, "numeric"))stop("The argument 'nb.int' must be a numeric")
if ((nb.int < 1)) stop("Number of time interval 'nb.int' must be between 1 and 20")
if (nb.int > 20){
nb.int <- 20
indic.nb.int <- 1 # equals 1 for nb.int > 20
}
else{
indic.nb.int <- 0 # equals 0 for nb.int < 20
}
nbintervR <- nb.int
size1 <- 3*nbintervR
}
if ((typeof == 0) | (typeof == 2)) indic.nb.int <- 0
if (sum(as.double(var))==0) nvar <- 0
if (timedep==0){
npbetatps <- 0
}else{
npbetatps <- (betaknots+betaorder-1)*nvartimedep
}
np <- switch(as.character(typeof),
"0"=(as.integer(uni.strat) * (as.integer(n.knots) + 2) + as.integer(nvar) + as.integer(Frailty)) + npbetatps,
"1"=(as.integer(uni.strat) * nbintervR + nvar + as.integer(Frailty)) + npbetatps,
"2"=(as.integer(uni.strat) * 2 + nvar + as.integer(Frailty)) + npbetatps)
# traitement de l'initialisation du Beta rentre par l'utilisateur
Beta <- rep(0,np)
if (!missing(init.B)) {
if (length(init.B) != nvar) stop("Wrong number of regression coefficients in init.B")
if (timedep) stop("You can hardly know initial regression coefficient while time-varying effect")
Beta <- c(rep(0,np-nvar),init.B)
}
if (!missing(init.Theta)) {
if (!is.numeric(init.Theta)) stop("init.Theta must be numeric")
if (Frailty==FALSE) stop("init.Theta does not exist in a Cox proportional hazard model")
Beta[np-nvar] <- init.Theta
}
xSuT <- matrix(0,nrow=100,ncol=uni.strat)
if (typeof==0){
mt1 <- size1
}else{
mt1 <- 100
}
size2 <- mt1
flush.console()
if (print.times){
ptm<-proc.time()
cat("\n")
cat("Be patient. The program is computing ... \n")
}
# typeof <- switch(hazard,"Splines"=0,"Piecewise"=1,"Weibull"=2)
familyrisk <- switch(family[1],"PH"=0,"PO"=1,"probit"=2, "AH"=3,"AH2"=4)
#cat("fichier GENFRAILTYPENAL.R : appel de la subroutine FRAILPENALGEN (dans frailtypackgen.f90) \n")
ans <- .Fortran(C_frailpenalgen,
as.integer(n),
as.integer(length(uni.cluster)),
as.integer(cens.data),
as.integer(uni.strat),
as.integer(Frailty),
as.integer(n.knots),
as.double(kappa),
as.double(tt0),
as.double(tt1),
as.integer(cens),
as.integer(cluster),
as.integer(nvar),
as.double(strats),
as.double(var),
as.integer(AG),
as.integer(noVar1),
as.integer(maxit),
as.integer(crossVal),
np=as.integer(np),
b=as.double(Beta),
as.double(matrix(0,nrow=np,ncol=np)),
as.double(matrix(0,nrow=np,ncol=np)),
as.double(0),
LCV=as.double(rep(0,2)),
as.double(matrix(0,nrow=size1,ncol=uni.strat)),
as.double(array(0,dim=c(size1,3,uni.strat))), #lam
xSuT=as.double(xSuT),
as.double(array(0,dim=c(size2,3,uni.strat))),
as.integer(typeof),
as.integer(equidistant),
as.integer(nbintervR),
as.integer(size1),
as.integer(0),
as.integer(0),
as.integer(0),
as.double(c(0,0)),
as.double(0),
istop=as.integer(0),
shape.weib=as.double(rep(0,2)),
scale.weib=as.double(rep(0,2)),
as.integer(mt1),
zi=as.double(rep(0,(n.knots+6))),
martingale.res=as.double(rep(0,as.integer(length(uni.cluster)))),
martingaleCox=as.double(rep(0,n)),
frailty.pred=as.double(rep(0,as.integer(length(uni.cluster)))),
frailty.var=as.double(rep(0,as.integer(length(uni.cluster)))),
frailty.sd=as.double(rep(0,as.integer(length(uni.cluster)))),
linear.pred=as.double(rep(0,n)),
time=as.double(rep(0,(nbintervR+1))),
as.integer(intcens), # rajout
as.double(ttU), # rajout
logNormal=as.integer(logNormal),
timedep=as.integer(timedep),
as.integer(betaknots),
as.integer(betaorder),
as.integer(filtretps),
BetaTpsMat=as.double(matrix(0,nrow=101,ncol=1+4*nvartimedep)),
EPS=as.double(c(LIMparam,LIMlogl,LIMderiv)),
nbgh = as.integer(nb.gh),
as.integer(familyrisk)
)#,
#PACKAGE = "frailtypack") # 60 arguments
#AD:
if (ans$istop == 4){
warning("Problem in the loglikelihood computation. The program stopped abnormally. Please verify your dataset. \n")
}
if (ans$istop == 2){
warning("Model did not converge.")
}
if (ans$istop == 3){
warning("Matrix non-positive definite.")
}
#AD:
if (noVar1 == 1) nvar<-0
np <- ans[[19]]
fit <- NULL
fit$b <- ans$b
fit$na.action <- attr(m, "na.action")
fit$call <- call
fit$n <- n
fit$groups <- length(uni.cluster)
fit$n.events <- ans[[34]]
#Al:
if(dim(table(cens,cluster))[1]==2) fit$n.eventsbygrp <- table(cens,cluster)[2,]
if(as.character(typeof)=="0"){
fit$logLikPenal <- ans[[23]]
}else{
fit$logLik <- ans[[23]]
}
if (Frailty) {
fit$coef <- ans[[20]][(np - nvar - npbetatps + 1):np]
if (logNormal == 0){fit$theta <- (ans[[20]][np - nvar - npbetatps])^2
}else{fit$sigma2 <- (ans[[20]][np - nvar - npbetatps])^2
}
}
if (!Frailty) {
if (logNormal == 0) fit$theta <- NULL
else fit$sigma2 <- NULL
}
if (noVar1 == 1) {
fit$coef <- NULL
}
else{
fit$coef <- ans[[20]][(np - nvar - npbetatps + 1):np]
noms <- factor.names(colnames(X))
if(length(grep(":",noms))>0)noms <- factor.names(noms)
if (timedep == 1){ # on enleve les parametres des B-splines qui ne serviront pas a l'utilisateur
while (length(grep("timedep",noms))!=0){
pos <- grep("timedep",noms)[1]
noms <- noms[-pos]
fit$coef <- fit$coef[-(pos:(pos+betaorder+betaknots-1))]
}
}
names(fit$coef) <- noms
}
temp1 <- matrix(ans[[21]], nrow = np, ncol = np)
temp2 <- matrix(ans[[22]], nrow = np, ncol = np)
if (Frailty) {
fit$varTheta <- c(temp1[(np - nvar - npbetatps),(np - nvar - npbetatps)],temp2[(np - nvar - npbetatps),(np - nvar - npbetatps)])
fit$varTheta <- ((2*ans[[20]][np - nvar - npbetatps])^2)*fit$varTheta # delta-method
if (logNormal == 0) fit$theta_p.value <- 1 - pnorm(fit$theta/sqrt(fit$varTheta[1]))
else fit$sigma2_p.value <- 1 - pnorm(fit$sigma2/sqrt(fit$varTheta[1]))
}
#AD:modification des dimensions des tableaux
if(nvar > 0){
fit$varH <- temp1[(np - nvar - npbetatps + 1):np, (np - nvar - npbetatps + 1):np]
fit$varHIH <- temp2[(np - nvar - npbetatps + 1):np, (np - nvar - npbetatps + 1):np]
noms <- factor.names(colnames(X))
if(length(grep(":",noms))>0)noms <- factor.names(noms)
if (timedep == 1){ # on enleve les variances des parametres des B-splines
while (length(grep("timedep",noms))!=0){
pos <- grep("timedep",noms)[1]
noms <- noms[-pos]
fit$varH <- fit$varH[-(pos:(pos+betaorder+betaknots-1)),-(pos:(pos+betaorder+betaknots-1))]
fit$varHIH <- fit$varHIH[-(pos:(pos+betaorder+betaknots-1)),-(pos:(pos+betaorder+betaknots-1))]
}
}
}
fit$varHtotal <- temp1 # new Al: 20/06/13
fit$varHIHtotal <- temp2
fit$formula <- formula(Terms)
fit$x <- matrix(ans[[25]], nrow = size1, ncol = uni.strat)
fit$lam <- if(typeof == 1){array(ans[[26]][seq(1,length(ans[[26]]),3)], dim = c(nb.int,3,uni.strat))} else{array(ans[[26]], dim = c(size1,3,uni.strat))} # Le lam s'ecrit differemment selon la fonction de hasard (Piecewise selon le nombre d'intervalles specifies, Weibull et Splines selon size1 = 100)
fit$xSu <- matrix(ans$xSuT, nrow = 100, ncol = uni.strat)
fit$surv <- array(ans[[28]], dim = c(size2,3,uni.strat))
fit$type <- type
fit$n.strat <- uni.strat
fit$n.iter <- ans[[33]]
median <- NULL
for (i in (1:fit$n.strat)) median[i] <- ifelse(typeof==0, minmin(fit$surv[,1,i],fit$x), minmin(fit$surv[,1,i],fit$xSu))
lower <- NULL
for (i in (1:fit$n.strat)) lower[i] <- ifelse(typeof==0, minmin(fit$surv[,2,i],fit$x), minmin(fit$surv[,2,i],fit$xSu))
upper <- NULL
for (i in (1:fit$n.strat)) upper[i] <- ifelse(typeof==0, minmin(fit$surv[,3,i],fit$x), minmin(fit$surv[,3,i],fit$xSu))
fit$median <- cbind(lower,median,upper)
if (typeof == 0){
fit$n.knots<-n.knots
if (uni.strat > 1) fit$kappa <- ans[[36]]
else fit$kappa <- ans[[36]][1]
fit$DoF <- abs(ans[[37]])
fit$cross.Val<-cross.validation
fit$n.knots.temp <- n.knots.temp
fit$zi <- ans$zi
}
if(typeof == 1) fit$time <- ans$time
#AD:
fit$LCV <- ans$LCV[1]
fit$AIC <- ans$LCV[2]
fit$npar <- np
fit$nvar <- nvar
fit$noVar1 <- noVar1
fit$indic.nb.int <- indic.nb.int
#AD:
if(ans[[35]]==2000) stop("The cross validation procedure cannot be finished. Try to change
either the number of knots or the seed for kappa parameter")
fit$typeof <- typeof
fit$equidistant <- equidistant
#fit$nbintervR <- nbintervR
fit$istop <- ans$istop
fit$AG <- recurrentAG
fit$intcens <- intcens # rajout
fit$logNormal <- ans$logNormal
fit$shape.param <- ans$shape.weib
fit$scale.param <- ans$scale.weib
fit$Names.data <- Names.data
if (Frailty) fit$Names.cluster <- Names.cluster
fit$Frailty <- Frailty
# if (Frailty){
# fit$martingale.res <- ans$martingale.res
# fit$frailty.pred <- ans$frailty.pred
# if (logNormal==0){
# fit$frailty.var <- ans$frailty.var
# fit$frailty.sd <- ans$frailty.sd
# }
# }
# else{
# fit$martingaleCox <- ans$martingaleCox
# }
if (Frailty) fit$linear.pred <- ans$linear.pred[order(ordre)] # pour remettre dans le bon ordre
else fit$linear.pred <- ans$linear.pred
fit$BetaTpsMat <- matrix(ans$BetaTpsMat,nrow=101,ncol=1+4*nvartimedep)
fit$nvartimedep <- nvartimedep
fit$Names.vardep <- vardep
fit$EPS <- ans$EPS
fit$family <- familyrisk
#
#========================= Test de Wald pour shared
if(ans$istop==1){
if ((length(vec.factor) > 0) & (timedep == 0)){
Beta <- ans[[20]][(np - nvar + 1):np]
VarBeta <- fit$varH
nfactor <- length(vec.factor)
p.wald <- rep(0,nfactor)
if(fit$istop == 1) fit$global_chisq <- waldtest(N=nvar,nfact=nfactor,place=ind.place,modality=occur,b=Beta,Varb=VarBeta)
else fit$global_chisq <- 0
fit$dof_chisq <- occur
fit$global_chisq.test <- 1
# Calcul de pvalue globale
for(i in 1:length(vec.factor)){
p.wald[i] <- signif(1 - pchisq(fit$global_chisq[i], occur[i]), 3)
}
fit$p.global_chisq <- p.wald
fit$names.factor <- vec.factor
}
else{
fit$global_chisq.test <- 0
}
}
if(!is.null(fit$coef)){
if(nvar != 1){
seH <- sqrt(diag(fit$varH))
}else{
seH <- sqrt(fit$varH)
}
fit$beta_p.value <- 1 - pchisq((fit$coef/seH)^2, 1)
}
#===============================================
if (length(Xlevels) >0)fit$Xlevels <- Xlevels
fit$contrasts <- contr.save
attr(fit,"joint")<-joint
attr(fit,"subcluster")<-FALSE
class(fit) <- "frailtyPenal"
} # End SHARED MODEL
#
# Begin JOINT MODEL
#
if (joint & !length(subcluster)){
# Preparing data ...
#AD:
if(Frailty =="FALSE"){
stop("For joint frailty models, 'Frailty' must be equal to 'TRUE' ")
}
#AD
if (classofY == "Surv"){
if (!recurrentAG){
if(joint.clust==0){
tempdc <- aggregate(tt1,by=list(num.id),FUN=sum)[,2]
lignedc0 <- length(tempdc)
tempdc <- cbind(rep(0,lignedc0),tempdc)
clusterdc <- aggregate(cluster,by=list(num.id),FUN=function(x) x[length(x)])[,2]
tt1.death <- 0
tt0.death <- 0
tt1 <- aggregate(tt1,by=list(num.id), FUN=function(x) x[1])[,2]
tt0 <- aggregate(tt0,by=list(num.id), FUN=function(x) x[1])[,2]
cluster <- aggregate(cluster,by=list(num.id), FUN=function(x) x[1])[,2]
table <- as.data.frame(cbind(tt0,tt1,cluster))
table <- table[order(table$cluster),]
cluster <- table$cluster
tt1 <- table$tt1
tt0 <- table$tt0
n <- length(tt0)
uni.cluster<-unique(num.id)#unique(cluster)
}
else{
tt1.death<-aggregate(tt1,by=list(cluster),FUN=sum)[,2]
tt0.death<-rep(0,length(tt1.death))
clusterdc <- 0
lignedc0 <- 0
tempdc <- 0
}
}
else{
if(joint.clust==0){
#tempdc <- aggregate(tt1,by=list(num.id,cluster),FUN=sum)[,2]
tempdc<-aggregate(tt1,by=list(num.id),FUN=function(x) x[length(x)])[,2]
lignedc0 <- length(tempdc)
tempdc <- cbind(rep(0,lignedc0),tempdc)
clusterdc <- aggregate(cluster,by=list(num.id),FUN=function(x) x[length(x)])[,2]
tt1.death <- 0
tt0.death <- 0
tt1 <- aggregate(tt1,by=list(num.id), FUN=function(x) x[1])[,2]
tt0 <- aggregate(tt0,by=list(num.id), FUN=function(x) x[1])[,2]
cluster <- aggregate(cluster,by=list(num.id), FUN=function(x) x[1])[,2]
table <- as.data.frame(cbind(tt0,tt1,cluster))
table <- table[order(table$cluster),]
cluster <- table$cluster
tt1 <- table$tt1
tt0 <- table$tt0
n <- length(tt0)
uni.cluster<-unique(num.id)#unique(cluster)
}else{
tt1.death<-aggregate(tt1,by=list(cluster),FUN=function(x) x[length(x)])[,2]
tt0.death<-rep(0,length(tt1.death))
clusterdc <- 0
lignedc0 <- 0
tempdc <- 0
}
}
}
else{ # Interval censoring
if (recurrentAG == TRUE) stop("You can't fit joint models on interval-censored data with recurrentAG = TRUE")
if (joint.clust==0){
tempdc0 <- aggregate(tt0,by=list(num.id),FUN=function(x) x[length(x)])[,2]
tempdc <- aggregate(tt1,by=list(num.id),FUN=function(x) x[length(x)])[,2]
lignedc0 <- length(tempdc)
tempdc <- cbind(tempdc0,tempdc)
clusterdc <- aggregate(cluster,by=list(num.id),FUN=function(x) x[length(x)])[,2]
tt1.death <- 0
tt0.death <- 0
# prendre en compte seulement un evenement pour le joint cluster
tt0 <- aggregate(tt0,by=list(num.id),FUN=function(x) x[1])[,2]
tt1 <- aggregate(tt1,by=list(num.id),FUN=function(x) x[1])[,2]
ttU <- aggregate(ttU,by=list(num.id),FUN=function(x) x[1])[,2]
cens <- aggregate(cens,by=list(num.id),FUN=function(x) x[1])[,2]
cluster <- aggregate(cluster,by=list(num.id),FUN=function(x) x[1])[,2]
if (!is.null(ncol(var))){ # if more than one covariate
varAG<-aggregate(var[,1],by=list(num.id), FUN=function(x) x[1])[,2]
if (ncol(var)>1){
for (i in 2:ncol(var)){
varAG.i<-aggregate(var[,i],by=list(num.id), FUN=function(x) x[1])[,2]
varAG<-cbind(varAG,varAG.i)
}
}
var<-varAG
}
else{
var<-aggregate(var,by=list(num.id), FUN=function(x) x[1])[,2]
}
nobs <- n
n <- length(tt0)
}
else{
tt1.death<-aggregate(tt1,by=list(cluster),FUN=function(x) x[length(x)])[,2]
tt0.death<-rep(0,length(tt1.death))
clusterdc <- 0
lignedc0 <- 0
tempdc <- 0
# prendre en compte seulement un evenement pour le joint
tt0 <- aggregate(tt0,by=list(cluster),FUN=function(x) x[1])[,2]
tt1 <- aggregate(tt1,by=list(cluster),FUN=function(x) x[1])[,2]
ttU <- aggregate(ttU,by=list(cluster),FUN=function(x) x[1])[,2]
cens <- aggregate(cens,by=list(cluster),FUN=function(x) x[1])[,2]
if (!is.null(ncol(var))){ # if more than one covariate
varAG<-aggregate(var[,1],by=list(cluster), FUN=function(x) x[1])[,2]
if (ncol(var)>1){
for (i in 2:ncol(var)){
varAG.i<-aggregate(var[,i],by=list(cluster), FUN=function(x) x[1])[,2]
varAG<-cbind(varAG,varAG.i)
}
}
var<-varAG
}
else{
var<-aggregate(var,by=list(cluster), FUN=function(x) x[1])[,2]
}
nobs <- n
n <- length(tt0)
}
}
Terms2 <- if (missing(data)){
if (!missing(formula.terminalEvent))terms(formula.terminalEvent, special)
}
else{
if (!missing(formula.terminalEvent))terms(formula.terminalEvent, special, data = data)
}
#AD:
if (!missing(formula.terminalEvent)){
ord2 <- attr(Terms2, "order")
if (length(ord2) & any(ord2 != 1)){
# stop("Interaction terms are not valid for terminal event formula")
}
}
#AD:
#AD:Joint model needs "terminal()"
if (ind.terminal){
if(joint.clust==0 ){
icdc00 <- aggregate(terminal,by=list(num.id),FUN=function(x) x[length(x)])[,2] #+aggregate(cens,by=list(num.id),FUN=function(x) x[length(x)])[,2]
terminalEvent <- 0
}
else{
terminalEvent<-aggregate(terminal,by=list(cluster),FUN=function(x) x[length(x)])[,2]
icdc00 <- 0
}
}
else{
stop(" Joint frailty model miss specified ")
}
#AD:
# terminalEvent might be 0-1
if(joint.clust==0){
if (!all(icdc00%in%c(2,1,0))){
stop("terminal must contain a variable coded 0-1 and a non-factor variable")
}
}
else{
if (!all(terminalEvent%in%c(2,1,0))){
stop("terminal must contain a variable coded 0-1 and a non-factor variable")
}
}
m2 <- match.call(expand.dots = FALSE)
## AD: modified 20 06 2011, for no covariates on terminal event part
if (missing(formula.terminalEvent)){
m2$n.knots <- m2$recurrentAG <- m2$cross.validation <- m2$kappa <- m2$maxit <- m2$hazard <- m2$nb.int <- m2$RandDist <- m2$betaorder <- m2$betaknots <- m2$init.B <- m2$LIMparam <- m2$LIMlogl <- m2$LIMderiv <- m2$print.times <- m2$init.Theta <- m2$init.Alpha <- m2$Alpha <- m2$init.Ksi <- m2$Ksi <- m2$init.Eta <- m2$Eta <- m2$initialize <- m2$nb.gh <- m2$nb.gl <- m2$family <- m2$... <- NULL
}else{
m2$formula.terminalEvent <- m2$n.knots <- m2$recurrentAG <- m2$cross.validation <- m2$jointGeneral<- m2$kappa <- m2$maxit <- m2$hazard <- m2$nb.int <- m2$RandDist <- m2$betaorder <- m2$betaknots <- m2$init.B <- m2$LIMparam <- m2$LIMlogl <- m2$LIMderiv <- m2$print.times <- m2$init.Theta <- m2$init.Alpha <- m2$Alpha <- m2$init.Ksi <- m2$Ksi <- m2$init.Eta <- m2$Eta <- m2$initialize <- m2$nb.gh <- m2$nb.gl <- m2$family <- m2$... <- NULL
}
m2$formula <- Terms2
m2[[1]] <- as.name("model.frame")
m2 <- eval(m2, sys.parent()) #ici il prend les colonne associe au var.exp, si rien il prend tout
match.noNA<-dimnames(m2)[[1]]%in%dimnames(m)[[1]]#masque logique pour ne pas avoir de NA
m2<-m2[match.noNA, ,drop=FALSE]#m2 inchanger si pas de NA
if (!missing(formula.terminalEvent))newTerms2<-Terms2
#=========================================================>
if (!missing(formula.terminalEvent)){
X2 <- model.matrix(newTerms2, m2)
lldc <- attr(newTerms2,"term.labels")
#ind.placedc <- grep("factor",lldc)
ind.placedc <- unique(attr(X2,"assign")[duplicated(attr(X2,"assign"))])#changement unique le 26/09/2014
vec.factordc <- NULL
vec.factordc <- c(vec.factordc,lldc[ind.placedc])
mat.factordc <- matrix(vec.factordc,ncol=1,nrow=length(vec.factordc))
# Fonction servant a prendre les termes entre "as.factor"
vec.factordc <-apply(mat.factordc,MARGIN=1,FUN=function(x){
if (length(grep("factor",x))>0){
if(length(grep(":",x))>0){
if(grep('\\(',unlist(strsplit(x,split="")))[1]<grep(":",unlist(strsplit(x,split="")))[1] && length(grep('\\(',unlist(strsplit(x,split=""))))==1){
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- grep(":",unlist(strsplit(x,split="")))[1]-2
pos3 <- grep(":",unlist(strsplit(x,split="")))[1]
pos4 <- length(unlist(strsplit(x,split="")))
return(paste(substr(x,start=pos1,stop=pos2),substr(x,start=pos3,stop=pos4),sep=""))
}
else if(grep("\\(",unlist(strsplit(x,split="")))[1]>grep(":",unlist(strsplit(x,split="")))[1] && length(grep('\\(',unlist(strsplit(x,split=""))))==1){
pos2 <- grep(":",unlist(strsplit(x,split="")))[1]
pos3 <- grep("\\(",unlist(strsplit(x,split="")))[1]+1
pos4 <- length(unlist(strsplit(x,split="")))-1
return(paste(substr(x,start=1,stop=pos2),substr(x,start=pos3,stop=pos4),sep=""))
}
else{#both factors
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- grep(":",unlist(strsplit(x,split="")))[1]-2
pos3 <- grep("\\(",unlist(strsplit(x,split="")))[2]+1
pos4 <- length(unlist(strsplit(x,split="")))-1
return(paste(substr(x,start=pos1,stop=pos2),":",substr(x,start=pos3,stop=pos4),sep=""))
}
}
else{
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- length(unlist(strsplit(x,split="")))-1
return(substr(x,start=pos1,stop=pos2))
}
}
else return(x)
} )
# On determine le nombre de categorie pour chaque var categorielle
if(length(vec.factordc) > 0){
vect.factdc <- attr(X2,"dimnames")[[2]]
vect.factdc <- vect.factdc[grep(paste(vec.factordc,collapse="|"),vect.factdc)]
occurdc <- rep(0,length(vec.factordc))
# for(i in 1:length(vec.factordc)){
# occurdc[i] <- length(grep(vec.factordc[i],vect.factdc))
# }
# }
interactiondc<-as.vector(apply(matrix(vect.factdc,nrow=length(vect.factdc)),MARGIN=1,FUN=function(x){length(grep(":",unlist(strsplit(x,split=""))))}))
which.interactiondc <- which(interactiondc==1)
for(i in 1:length(vec.factordc)){
if(length(grep(":",unlist(strsplit(vec.factordc[i],split=""))))>0){
pos <- grep(":",unlist(strsplit(vec.factordc[i],split="")))
length.grep <- 0
for(j in 1:length(vect.factdc)){
if(j%in%which.interactiondc){
if(length(grep(substr(vec.factordc[i],start=1,stop=pos-1),vect.factdc[j]))>0 && length(grep(substr(vec.factordc[i],start=pos+1,stop=length(unlist(strsplit(vec.factordc[i],split="")))),vect.factdc[j]))>0){
length.grep <- length.grep + 1
which <- i
}
}
}
occurdc[i] <- length.grep
}
else{
if(length(vect.factdc[-which.interactiondc])>0) occurdc[i] <- length(grep(vec.factordc[i],vect.factdc[-which.interactiondc]))
else occurdc[i] <- length(grep(vec.factordc[i],vect.factdc))
}
}
}
#=========================================================>
assign <- lapply(attrassign(X2, newTerms2)[-1], function(x) x - 1)
Xlevels2 <- .getXlevels(newTerms2, m2)
contr.save2 <- attr(X2, 'contrasts')
#========================================>
if(length(vec.factordc) > 0){
positiondc <- unlist(assign,use.names=F)
}
#========================================>
if (ncol(X2) == 1){
X2<-X2-1
noVar2 <- 1
}
else{
X2 <- X2[, -1, drop = FALSE]
noVar2 <- 0
}
nvar2 <- ncol(X2)
# if(sum(ord)>length(ord)){
# for(i in 1:length(ord)){
# if(ord[i]>1){
# name_v1 <- strsplit(as.character(lldc[i]),":")[[1]][1]
# name_v2 <- strsplit(as.character(lldc[i]),":")[[1]][2]
# if(length(grep("as.factor",name_v1))>0){name_v1<-substring(name_v1,11,nchar(name_v1)-1)
# v1 <- as.factor(data[,names(data)==name_v1])}
# else{v1 <- data[,names(data)==name_v1]}
# if(length(grep("as.factor",name_v2))>0){name_v2<-substring(name_v2,11,nchar(name_v2)-1)
# v2 <- as.factor(data[,names(data)==name_v2])}
# else{v2 <- data[,names(data)==name_v2]}
#
# if(is.factor(v1) && length(levels(v1))>2)stop("Interactions not allowed for factors with 3 or more levels (yet)")
# if(is.factor(v2) && length(levels(v2))>2)stop("Interactions not allowed for factors with 3 or more levels (yet)")
#
#
# }
# }
# }
vartimedep2 <- attr(Terms2, "specials")$timedep #nbre de var en fonction de timedep()
# verifier qu'il y ait du timedep dans la deuxieme formule
if (!is.null(vartimedep2)) timedep <- 1
varnotdep2 <- colnames(X2)[-grep("timedep",colnames(X2))]
vardep2 <- colnames(X2)[grep("timedep",colnames(X2))]
vardep2 <- apply(matrix(vardep2,ncol=1,nrow=length(vardep2)),1,timedep.names)
if (length(intersect(varnotdep2,vardep2)) != 0) {
stop("A variable is both used as a constant and time-varying effect covariate in the formula of terminal event")
}
nvartimedep2 <- length(vardep2)
filtretps2 <- rep(0,nvar2)
filtretps2[grep("timedep",colnames(X2))] <- 1
vardc.temp<-matrix(c(X2),nrow=nrow(X2),ncol=nvar2)
if(is.null(nrow(m2))){
if (length(m2) != nrow(m)){
stop(" There are missing values in the covariates modelling the terminal event. \n Prepare data only with complete cases")
}
}
else{
if (nrow(m2) != nrow(m)){
stop(" There are missing values in the covariates modelling the terminal event. \n Prepare data only with complete cases")
}
}
if(joint.clust==0 ){
if (!is.null(ncol(vardc.temp))){
vaxdc00<-aggregate(vardc.temp[,1],by=list(num.id), FUN=function(x) x[length(x)])[,2] # num.id au lieu de cluster
if (ncol(vardc.temp)>1){
for (i in 2:ncol(vardc.temp)){
vaxdc00.i<-aggregate(vardc.temp[,i],by=list(num.id), FUN=function(x) x[length(x)])[,2]
vaxdc00<-cbind(vaxdc00,vaxdc00.i)
}
}
}
else{
vaxdc00<-aggregate(vardc.temp,by=list(num.id), FUN=function(x) x[length(x)])[,2]
}
vardc <- 0
}
else{
if (!is.null(ncol(vardc.temp))){
vardc<-aggregate(vardc.temp[,1],by=list(cluster), FUN=function(x) x[length(x)])[,2]
if (ncol(vardc.temp)>1){
for (i in 2:ncol(vardc.temp)){
vardc.i<-aggregate(vardc.temp[,i],by=list(cluster), FUN=function(x) x[length(x)])[,2]
vardc<-cbind(vardc,vardc.i)
}
}
}
else{
vardc<-aggregate(vardc.temp,by=list(cluster), FUN=function(x) x[length(x)])[,2]
}
vaxdc00 <- 0
}
}
else{
noVar2 <- 1
vardc<-0
}
if ((classofY == "SurvIC") & (joint.clust==1 || joint.clust== 2) & (recurrentAG=FALSE)) cluster <- aggregate(cluster,by=list(cluster),FUN=function(x) x[1])[,2]
# cluster <- aggregate(cluster,by=list(cluster),FUN=function(x) x[1])[,2]
nvarRec<-nvar
if (!missing(formula.terminalEvent)){
#=======================================>
#======= Construction du vecteur des indicatrice
if(length(vec.factordc) > 0){
k <- 0
for(i in 1:length(vec.factordc)){
ind.placedc[i] <- ind.placedc[i]+k
k <- k + occurdc[i]-1
}
}
#==================================
if(joint.clust==1 || joint.clust==2){
if(is.null(nrow(vardc))){
nvarEnd<-1
}
else{
nvarEnd<-ncol(vardc)
}
}
else{
if(is.null(nrow(vaxdc00))){
nvarEnd<-1
}
else{
nvarEnd<-ncol(vaxdc00)
}
}
}
else nvarEnd<-0
if (sum(as.double(var))==0) nvarRec <- 0
if ((joint.clust==0) & sum(as.double(vaxdc00))==0) nvarEnd <- 0
if ((joint.clust==1 || joint.clust==2 ) & sum(as.double(vardc))==0) nvarEnd <- 0
nvar<-nvarRec+nvarEnd
# ... end preparing data
#AD:
effet <- 1
indic_alpha <- 1
indic_xi <- 3
if (!missing(Alpha)){ # new : joint more flexible alpha = 0
if (Alpha=="None") indic_alpha <- 0
else stop("Alpha can only take 'None' as a value in this version of frailtypack package")
}
nst <- uni.strat #2
indices <- c(indic_alpha, indic_xi)
if((equidistant %in% c(0,1)) & (typeof == 1)){
if (missing(nb.int)) stop("Number of time interval for recurrences and terminal event 'nb.int' is required")
if (length(nb.int) != 2) stop("The length of argument 'nb.int' should be 2. Must indicate for both recurrent events and terminal event.")
if (!inherits(nb.int, "numeric"))stop("The argument 'nb.int' must be a numeric")
if (nb.int[1] < 1) stop("Number of time interval 'nb.int' must be between 1 and 20")
if (nb.int[2] < 1) stop("Number of time interval 'nb.int' must be between 1 and 20")
if (nb.int[1] > 20){
nb.int[1] <- 20
indic.nb.int1 <- 1 # equals 1 for nb.int1 > 20
}
else{
indic.nb.int1 <- 0 # equals 1 for nb.int1 < 20
}
if (nb.int[2] > 20){
nb.int[2] <-20
indic.nb.int2 <- 1 # equals 1 for nb.int1 > 20
}
else{
indic.nb.int2 <- 0 # equals 1 for nb.int1 < 20
}
nbintervR <- nb.int[1]
size1 <- 3*nbintervR
nbintervDC <- nb.int[2]
size2 <- 3*nbintervDC
}
if ((typeof == 0) | (typeof == 2)){
indic.nb.int1 <- 0
indic.nb.int2 <- 0
}
if (timedep==0){
npbetatps1 <- 0
npbetatps2 <- 0
}
else{
npbetatps1 <- (betaknots+betaorder-1)*nvartimedep
npbetatps2 <- (betaknots+betaorder-1)*nvartimedep2
}
npbetatps <- npbetatps1 + npbetatps2
np <- switch(as.character(typeof),
"0"=((nst+1) * (n.knots + 2) + nvarRec + nvarEnd + effet + indic_alpha + npbetatps),
"1"=(nst*nbintervR + nbintervDC + nvarRec + nvarEnd + effet + indic_alpha + npbetatps),
"2"=(2*(nst+1) + nvarRec + nvarEnd + effet + indic_alpha + npbetatps))
if (all(all.equal(as.numeric(cens),terminal)==T)){
stop("'Recurrent event' variable and 'Terminal event' variable need to be different")
}
# traitement de l'initialisation du Beta rentre par l'utilisateur
Beta <- rep(0,np)
if (!missing(init.B)){
if (length(init.B) != nvar) stop("Wrong number of regression coefficients in init.B")
if (timedep) stop("You can hardly know initial regression coefficient while time-varying effect")
Beta <- c(rep(0,np-nvar),init.B)
}
if (!missing(init.Theta)){
if (!is.numeric(init.Theta)) stop("init.Theta must be numeric")
Beta[np-nvar-indic_alpha] <- init.Theta
}
if (!missing(init.Alpha)){
if (!missing(Alpha)) stop("You can not both initialize alpha parameter and fit a joint model without it")
if (!is.numeric(init.Alpha)) stop("init.Alpha must be numeric")
Beta[np-nvar] <- init.Alpha
}
xSu1 <- rep(0,100)
xSu2 <- rep(0,100)
if (typeof==0){
mt11 <- size1
mt12 <- size2
}
else{
mt11 <- 100
mt12 <- 100
}
initialize <- 1
#npinit <- switch(as.character(typeof),
# "0"=((n.knots + 2) + nvarRec + effet),
# "1"=(nbintervR + nvarRec + nvarEnd + effet),
# "2"=(2 + nvarRec + effet))
if ((uni.strat > 1 || joint.clust==2) & (joint.clust==0)) stop("stratification for clustered joint model is not yet allowed")
if ((uni.strat > 1 || joint.clust==2) & (intcens)) stop("stratification for joint model with interval censored data is not yet allowed")
if ((uni.strat > 1 || joint.clust==2) & (timedep)) stop("stratification for joint model and time-varying effect of covariates are not yet allowed")
if ((typeof==0) & (length(kappa)!=(uni.strat+1))) stop("wrong length of argument kappa")
ordretmp <- rep(0, n)
flush.console()
if (print.times){
ptm<-proc.time()
cat("\n")
cat("Be patient. The program is computing ... \n")
}
familyriskdc <- switch(family[1],"PH"=0,"PO"=1,"probit"=2, "AH"=3, "AH2"=4)
familyriskrec <- switch(family[2],"PH"=0,"PO"=1,"probit"=2, "AH"=3, "AH2"=4)
familyrisk <- as.integer(c(familyriskdc,familyriskrec))
ans <- .Fortran(C_jointgen,
as.integer(n),
as.integer(c(length(uni.cluster),0, uni.strat)), #IJ
#as.integer(uni.strat),
as.integer(strats),
as.integer(lignedc0),###
as.integer(n.knots),
#k0=as.double(kappa), # joint intcens,tps,cluster
axT=as.double(kappa), # joint avec generalisation de strate
as.double(tt0),
as.double(tt1), #8
as.integer(cens),
as.integer(cluster),
as.integer(clusterdc),
as.integer(0),###
as.double(tt0.death),
as.double(tt1.death),
as.integer(terminalEvent),
as.double(tempdc),###
as.integer(icdc00),###
as.integer(nvarRec),
as.double(var), #19
as.integer(nvarEnd),
as.double(vardc),
as.double(vaxdc00),###
as.integer(c(noVar1,noVar2)),
#as.integer(noVar2),
as.double(weights.agg),
as.integer(maxit),
np=as.integer(np),
b=as.double(Beta),
H=as.double(matrix(0,nrow=np,ncol=np)),
HIH=as.double(matrix(0,nrow=np,ncol=np)), #29
loglik=as.double(0),
LCV=as.double(rep(0,2)),
xR=as.double(matrix(0,nrow=size1,ncol=uni.strat)),
lamR=as.double(array(0,dim=c(size1,3,uni.strat))),
xSuR=as.double(xSu1),
survR=as.double(array(0,dim=c(mt11,3,uni.strat))),
xD=as.double(rep(0,size2)),
lamD=as.double(matrix(0,nrow=size2,ncol=3)),
xSuD=as.double(xSu2),
survD=as.double(matrix(0,nrow=mt12,ncol=3)), #39
as.integer(c(typeof, equidistant)),
#as.integer(equidistant),
as.integer(c(nbintervR, nbintervDC)),
#as.integer(nbintervDC),
as.integer(c(size1,size2,mt11,mt12)),###
counts=as.integer(c(0,0,0,0)),
IerIstop = as.integer(c(0,0)),
#ier=as.integer(0),
#istop=as.integer(0),
paraweib=as.double(rep(0,4)),
# shape.weib=as.double(rep(0,2)),
# scale.weib=as.double(rep(0,2)),
MartinGale=as.double(matrix(0,nrow=length(uni.cluster),ncol=5)),###46
linear.pred=as.double(rep(0,n)),
lineardc.pred=as.double(rep(0,as.integer(length(uni.cluster)))),
zi=as.double(rep(0,(n.knots+6))),
time=as.double(rep(0,(nbintervR+1))),
timedc=as.double(rep(0,(nbintervDC+1))),
# kendall=as.double(matrix(0,nrow=4,ncol=2)),
# as.integer(initialize),
# as.integer(npinit),
# Bshared=as.double(rep(0,npinit)),
linearpredG=as.double(rep(0,lignedc0)),
joint.clust=as.integer(joint.clust),
as.integer(intcens),
as.integer(indices),
#as.integer(indic_alpha),
#as.integer(indic_xi),
# censure par intervalle, indic_alpha
as.double(ttU),
as.integer(ordretmp),
as.integer(c(initialize,logNormal)),#58
#logNormal=as.integer(logNormal),
paratps=as.integer(c(timedep,betaknots,betaorder)),
as.integer(c(filtretps,filtretps2)),
BetaTpsMat=as.double(matrix(0,nrow=101,ncol=1+4*nvartimedep)),
BetaTpsMatDc=as.double(matrix(0,nrow=101,ncol=1+4*nvartimedep2)),
EPS=as.double(c(LIMparam,LIMlogl,LIMderiv)),
nbgauss = as.integer(c(nb.gh,nb.gl)),
as.integer(familyrisk)
)#,
#PACKAGE = "frailtypack") # 65 arguments
MartinGale <- matrix(ans$MartinGale,nrow=as.integer(length(uni.cluster)),ncol=5)
istop <- ans$IerIstop[2]
if (istop == 4){
warning("Problem in the loglikelihood computation. The program stopped abnormally. Please verify your dataset. \n")
}
if (istop == 2){
warning("Model did not converge.")
}
if (istop == 3){
warning("Matrix non-positive definite.")
}
#AD:
if (noVar1==1 & noVar2==1) nvar<-0
#AD:
np <- ans$np
fit <- NULL
fit$b <- ans$b
fit$na.action <- attr(m, "na.action")
fit$call <- call
if (classofY == "SurvIC"){
fit$n <- nobs
if (typeofY == "intervaltronc") fit$indic.trunc <- 1
else fit$indic.trunc <- 0
}else{
fit$n <- n
}
if (joint.clust == 0) fit$ind <- lignedc0
fit$groups <- length(uni.cluster)
fit$n.events <- ans$counts[2]
fit$n.deaths <- ans$counts[3]
fit$n.censored<-ans$counts[4]
if(as.character(typeof)=="0"){
fit$logLikPenal <- ans$loglik
}
else{
fit$logLik <- ans$loglik
}
#AD:
fit$LCV <- ans$LCV[1]
fit$AIC <- ans$LCV[2]
#AD:
if (logNormal == 0) fit$theta <- ans$b[np - nvar - npbetatps - indic_alpha]^2
else fit$sigma2 <- ans$b[np - nvar - npbetatps - indic_alpha]^2
if (indic_alpha == 1) fit$alpha <- ans$b[np - nvar - npbetatps]
if (joint.clust==2) fit$eta <- ans$b[np - nvar]^2
fit$npar <- np
#AD:
if ((noVar1==1 & noVar2==1)){
fit$coef <- NULL
}else{
fit$coef <- ans$b[(np - nvar - npbetatps + 1):np]
noms <- c(factor.names(colnames(X)),factor.names(colnames(X2)))
if(length(grep(":",noms))>0)noms <- factor.names(noms)
if (timedep == 1){
while (length(grep("timedep",noms))!=0){
pos <- grep("timedep",noms)[1]
noms <- noms[-pos]
fit$coef <- fit$coef[-(pos:(pos+betaorder+betaknots-1))]
}
}
names(fit$coef) <- noms
# if (missing(formula.terminalEvent)){
# names(fit$coef) <- c(factor.names(colnames(X)))
# }else{
# names(fit$coef) <- c(factor.names(colnames(X)), factor.names(colnames(X2)))
# }
}
#AD:
temp1 <- matrix(ans$H, nrow = np, ncol = np)
temp2 <- matrix(ans$HIH, nrow = np, ncol = np)
#Al:
fit$varHtotal <- temp1
fit$varHIHtotal <- temp2
fit$varH <- temp1[(np - nvar - npbetatps - indic_alpha):np, (np - nvar - npbetatps - indic_alpha):np]
fit$varHIH <- temp2[(np - nvar - npbetatps - indic_alpha):np, (np - nvar - npbetatps - indic_alpha):np]
if (indic_alpha == 1)fit$alpha_p.value <- 1 - pchisq((fit$alpha/sqrt(diag(fit$varH))[2])^2,1)
if (logNormal == 0){seH.frail <- sqrt(((2 * (fit$theta^0.5))^2) * diag(fit$varH)[1])
fit$theta_p.value <- 1 - pnorm(fit$theta/seH.frail)
}else{ seH.frail <- sqrt(((2 * (fit$sigma2^0.5))^2) * diag(fit$varH)[1])
fit$sigma2_p.value <- 1 - pnorm(fit$sigma2/seH.frail)
}
if (indic_alpha == 1) noms <- c("theta","alpha",factor.names(colnames(X)),factor.names(colnames(X2)))
else noms <- c("theta",factor.names(colnames(X)),factor.names(colnames(X2)))
if(length(grep(":",noms))>0)noms <- factor.names(noms)
if (timedep == 1){ # on enleve les variances des parametres des B-splines
while (length(grep("timedep",noms))!=0){
pos <- grep("timedep",noms)[1]
noms <- noms[-pos]
fit$varH <- fit$varH[-(pos:(pos+betaorder+betaknots-1)),-(pos:(pos+betaorder+betaknots-1))]
fit$varHIH <- fit$varHIH[-(pos:(pos+betaorder+betaknots-1)),-(pos:(pos+betaorder+betaknots-1))]
}
}
fit$nvar<-c(nvarRec,nvarEnd)
fit$nvarnotdep<-c(nvarRec-nvartimedep,nvarEnd-nvartimedep2)
# fit$formula <- formula(Terms)
fit$xR <- matrix(ans$xR, nrow = size1, ncol = uni.strat)
fit$lamR <- if(typeof == 1){array(ans$lamR[seq(1,length(ans$lamR),3)], dim = c(nb.int[1],3,uni.strat))} else{array(ans$lamR, dim = c(size1,3,uni.strat))}
fit$xSuR <- matrix(ans$xSuR, nrow = 100, ncol = uni.strat)
fit$survR <- array(ans$survR, dim = c(mt11,3,uni.strat))
fit$xD <- ans$xD
fit$lamD <- if(typeof == 1) {matrix(ans$lamD[seq(1,length(ans$lamD),3)], nrow = nb.int[2], ncol = 3)} else{matrix(ans$lamD, nrow = size2, ncol = 3)}
fit$xSuD <- ans$xSuD
fit$survD <- matrix(ans$survD, nrow = mt12, ncol = 3)
fit$type <- type
fit$n.strat <- uni.strat
fit$n.iter <- ans$counts[1]
fit$typeof <- typeof
if (typeof == 0){
fit$n.knots<-n.knots
fit$kappa <- ans$axT
fit$cross.Val<-cross.validation
fit$n.knots.temp <- n.knots.temp
fit$zi <- ans$zi
}
if(typeof == 1){
fit$time <- ans$time
fit$timedc <- ans$timedc
}
fit$nb.gh <- nb.gh
fit$nb.gl <- nb.gl
medianR <- NULL
for (i in (1:fit$n.strat)) medianR[i] <- ifelse(typeof==0, minmin(fit$survR[,1,i],fit$xR), minmin(fit$survR[,1,i],fit$xSuR))
lowerR <- NULL
for (i in (1:fit$n.strat)) lowerR[i] <- ifelse(typeof==0, minmin(fit$survR[,2,i],fit$xR), minmin(fit$survR[,2,i],fit$xSuR))
upperR <- NULL
for (i in (1:fit$n.strat)) upperR[i] <- ifelse(typeof==0, minmin(fit$survR[,3,i],fit$xR), minmin(fit$survR[,3,i],fit$xSuR))
fit$medianR <- cbind(lowerR,medianR,upperR)
medianD <- ifelse(typeof==0, minmin(fit$survD[,1],fit$xD), minmin(fit$survD[,1],fit$xSuD))
lowerD <- ifelse(typeof==0, minmin(fit$survD[,3],fit$xD), minmin(fit$survD[,3],fit$xSuD))
upperD <- ifelse(typeof==0, minmin(fit$survD[,2],fit$xD), minmin(fit$survD[,2],fit$xSuD))
fit$medianD <- cbind(lowerD,medianD,upperD)
#AD:
fit$noVar1 <- noVar1
fit$noVar2 <- noVar2
#fit$nbintervR <- nbintervR
#fit$nbintervDC <- nbintervDC
fit$nvarRec <- nvarRec
fit$nvarEnd <- nvarEnd
fit$istop <- istop
fit$indic.nb.intR <- indic.nb.int1
fit$indic.nb.intD <- indic.nb.int2
fit$shape.param <- ans$paraweib[1:2]#ans$shape.weib
fit$scale.param <- ans$paraweib[3:4]#ans$scale.weib
#AD:
# verif que les martingales ont ete bien calculees
msg <- "Problem in the estimation of the random effects (perhaps high number of events in some clusters)"
if (Frailty){
if (any(MartinGale[,1]==0)){
#fit$martingale.res <- msg
#fit$martingaledeath.res <- msg
#fit$frailty.pred <- msg
# fit$frailty.var <- msg
fit$linear.pred <- msg
fit$lineardeath.pred <- msg
}
else{
#fit$martingale.res <- MartinGale[,1]#ans$martingale.res
#fit$martingaledeath.res <- MartinGale[,2]#ans$martingaledc.res
#fit$frailty.pred <- MartinGale[,3]#ans$frailty.pred
# fit$frailty.var <- MartinGale[,4]#ans$frailty.var
fit$linear.pred <- ans$linear.pred[order(ordre)]
if (joint.clust==0) fit$lineardeath.pred <- ans$linearpredG
else fit$lineardeath.pred <- ans$lineardc.pred
}
}
# if (joint.clust==0){
# fit$kendall <- matrix(ans$kendall,nrow=4,ncol=2)
# }
fit$joint.clust <- ans$joint.clust
fit$AG <- recurrentAG
fit$intcens <- intcens # rajout
fit$indic_alpha <- indic_alpha
if(joint.clust==2)fit$indic_alpha <- 0
fit$logNormal <- logNormal#ans$logNormal
fit$BetaTpsMat <- matrix(ans$BetaTpsMat,nrow=101,ncol=1+4*nvartimedep)
fit$BetaTpsMatDc <- matrix(ans$BetaTpsMatDc,nrow=101,ncol=1+4*nvartimedep2)
fit$nvartimedep <- c(nvartimedep,nvartimedep2)
fit$Names.vardep <- vardep
fit$Names.vardepdc <- vardep2
fit$EPS <- ans$EPS
fit$family <- familyrisk
#================================> For the reccurrent
#========================= Test de Wald
if ((length(vec.factor) > 0) & (timedep == 0)){
Beta <- ans$b[(np - nvar + 1):np]
if (indic_alpha == 1) VarBeta <- diag(diag(fit$varH)[-c(1,2)])
else VarBeta <- diag(diag(fit$varH)[-c(1)])
nfactor <- length(vec.factor)
p.wald <- rep(0,nfactor)
ntot <- nvarEnd + nvarRec
if(fit$istop == 1) fit$global_chisq <- waldtest(N=nvarRec,nfact=nfactor,place=ind.place,modality=occur,b=Beta,Varb=VarBeta,Llast=nvarEnd,Ntot=ntot)
else fit$global_chisq <- 0
fit$dof_chisq <- occur
fit$global_chisq.test <- 1
# Calcul de pvalue globale
for(i in 1:length(vec.factor)){
p.wald[i] <- signif(1 - pchisq(fit$global_chisq[i], occur[i]), 3)
}
fit$p.global_chisq <- p.wald
fit$names.factor <- vec.factor
}
else{
fit$global_chisq.test <- 0
}
#================================> For the death
#========================= Test de Wald
if (!missing(formula.terminalEvent)){
if ((length(vec.factordc) > 0) & (timedep == 0)){
Beta <- ans$b[(np - nvar + 1):np]
#if (indic_alpha == 1) VarBeta <- diag(diag(fit$varH)[-c(1,2)])
#else VarBeta <- diag(diag(fit$varH)[-c(1)])
if(indic_alpha == 1) VarBeta <- fit$varH[3:dim(fit$varH)[1],3:dim(fit$varH)[2]]
else VarBeta <- fit$varH[2:dim(fit$varH)[1],2:dim(fit$varH)[2]]
nfactor <- length(vec.factordc)
p.walddc <- rep(0,nfactor)
ntot <- nvarEnd + nvarRec
if(fit$istop == 1) fit$global_chisq_d <- waldtest(N=nvarEnd,nfact=nfactor,place=ind.placedc,modality=occurdc,b=Beta,Varb=VarBeta,Lfirts=nvarRec,Ntot=ntot)
else fit$global_chisq_d <- 0
fit$dof_chisq_d <- occurdc
fit$global_chisq.test_d <- 1
# Calcul de pvalue globale
for(i in 1:length(vec.factordc)){
p.walddc[i] <- signif(1 - pchisq(fit$global_chisq_d[i], occurdc[i]), 3)
}
fit$p.global_chisq_d <- p.walddc
fit$names.factordc <- vec.factordc
}
else{
fit$global_chisq.test_d <- 0
}
}
else{
fit$global_chisq.test_d <- 0
}
if(!is.null(fit$coef)){
if (indic_alpha == 1 || fit$joint.clust==2) {
seH <- sqrt(diag(fit$varH))[-c(1,2)]
}else{
seH <- sqrt(diag(fit$varH))[-1]
}
fit$beta_p.value <- 1 - pchisq((fit$coef/seH)^2, 1)
}
if (length(Xlevels) >0) fit$Xlevels <- Xlevels
fit$contrasts <- contr.save
if (length(Xlevels2) >0) fit$Xlevels2 <- Xlevels2
fit$contrasts2 <- contr.save2
#NCC design
fit$ncc <- FALSE
if(length(wts))fit$ncc <- TRUE
fit$formula <- formula
fit$formula.terminalEvent <- formula.terminalEvent
attr(fit,"joint")<-joint
attr(fit,"subcluster")<-FALSE
class(fit) <- "jointPenal"
} # End JOINT MODEL
#
# Begin NESTED MODEL
#
effet <- 1
# Modified ML 24/03/2015 for Nested Joint model
if (length(subcluster) & !joint){
if (logNormal == 1) stop("Nested model not implemented yet for log normal distribution of frailties")
if((equidistant %in% c(0,1)) & (typeof == 1)){
if (missing(nb.int)) stop("Number of time interval 'nb.int' is required")
if (length(nb.int) != 1) stop("Wrong length of number of time interval argument 'nb.int'")
if (!inherits(nb.int, "numeric"))stop("The argument 'nb.int' must be a numeric")
if (nb.int < 1) stop("Number of time interval 'nb.int' must be between 1 and 20")
if (nb.int > 20){
nb.int <-20
indic.nb.int <- 1 # equals 1 for nb.int > 20
}
else{
indic.nb.int <- 0 # equals 1 for nb.int < 20
}
nbintervR <- nb.int
size1 <- 3*nbintervR
}
if ((typeof == 0) | (typeof == 2)) indic.nb.int <- 0
if (sum(as.double(var))==0) nvar <- 0
np <- switch(as.character(typeof),
"0"=(as.integer(uni.strat) * (as.integer(n.knots) + 2) + as.integer(nvar) + 2 * as.integer(Frailty)),
"1"=(as.integer(uni.strat) * nbintervR + nvar + 2 * as.integer(Frailty)),
"2"=(as.integer(uni.strat) * 2 + nvar + 2 * as.integer(Frailty)))
xSu1 <- rep(0,100)
xSu2 <- rep(0,100)
if (typeof==0){
mt1 <- size1
}
else{
mt1 <- 100
}
size2 <- mt1
if (length(kappa)==1) kappa <- c(kappa,0)
########### group and subgroup
grpe <- function(g){
grp <- unique(g)
res <- rep(0,length(grp))
for(i in 1:length(res)){
res[i] = sum(grp[i]==g)
}
return(res)
}
grp <- grpe(as.integer(cluster))
subgrpe <- function(g,sg){
j <- 0
k <- 0
res <- rep(0,length(g))
for(i in 1:length(g)){
k <- k + g[i]
j <- j + 1
temp <- sg[j:k]
res[i] <- length(grpe(temp))
j <- k
}
return(res)
}
subgbyg <- subgrpe(grp,as.integer(subcluster))
maxng <- max(subgbyg)
ngg <- length(uni.cluster)
flush.console()
if (print.times){
ptm<-proc.time()
cat("\n")
cat("Be patient. The program is computing ... \n")
}
ans <- .Fortran(C_nested,
as.integer(n),
as.integer(length(uni.cluster)),
as.integer(length(uni.subcluster)),
as.integer(uni.strat),
as.integer(n.knots),
as.double(kappa),
as.double(tt0),
as.double(tt1),
as.integer(cens),
as.integer(cluster),
as.integer(subcluster),
as.integer(nvar),
as.double(strats),
as.double(var),
as.integer(AG),
as.integer(noVar1),
as.integer(maxit),
as.integer(crossVal),
as.integer(np),
as.integer(maxng),
b=as.double(rep(0,np)),
H=as.double(matrix(0,nrow=np,ncol=np)),
HIH=as.double(matrix(0,nrow=np,ncol=np)),
loglik=as.double(0),
LCV=as.double(rep(0,2)),
x1=as.double(rep(0,size1)),
lam=as.double(matrix(0,nrow=size1,ncol=3)),
xSu1=as.double(xSu1),
surv=as.double(matrix(0,nrow=size2,ncol=3)),
x2=as.double(rep(0,size1)),
lam2=as.double(matrix(0,nrow=size1,ncol=3)),
xSu2=as.double(xSu2),
surv2=as.double(matrix(0,nrow=size2,ncol=3)),
as.integer(typeof),
as.integer(equidistant),
as.integer(nbintervR),
as.integer(size1),
ni=as.integer(0),
cpt=as.integer(0),
ier=as.integer(0),
k0=as.double(c(0,0)),
ddl=as.double(0),
istop=as.integer(0),
shape.weib=as.double(rep(0,2)),
scale.weib=as.double(rep(0,2)),
as.integer(mt1),
zi=as.double(rep(0,(n.knots+6))),
time=as.double(rep(0,(nbintervR+1))),
martingale.res=as.double(rep(0,as.integer(length(uni.subcluster)))),
frailty.pred.group=as.double(rep(0,as.integer(length(uni.cluster)))),
frailty.pred.subgroup=as.double(matrix(0,nrow=ngg,ncol=maxng)),
frailty.var.group=as.double(rep(0,as.integer(length(uni.cluster)))),
frailty.var.subgroup=as.double(matrix(0,nrow=ngg,ncol=maxng)),
frailty.sd.group=as.double(rep(0,as.integer(length(uni.cluster)))),
frailty.sd.subgroup=as.double(matrix(0,nrow=ngg,ncol=maxng)),
linear.pred=as.double(rep(0,n)),
EPS=as.double(c(LIMparam,LIMlogl,LIMderiv)),
nbgl = as.integer(nb.gl)
)#,
#PACKAGE = "frailtypack") # 57 arguments
if (ans$istop == 4){
warning("Problem in the loglikelihood computation. The program stopped abnormally. Please verify your dataset. \n")
}
if (ans$istop == 2){
warning("Model did not converge.")
}
if (ans$istop == 3){
warning("Matrix non-positive definite.")
}
nst <- as.integer(uni.strat)
if (noVar1 == 1) nvar<-0
np <- np
fit <- NULL
fit$b <- ans$b
fit$na.action <- attr(m, "na.action")
fit$call <- call
fit$n <- n
fit$groups <- length(uni.cluster)
fit$subgroups <- length(uni.subcluster)
fit$n.events <- ans$cpt
if(as.character(typeof)=="0"){
fit$logLikPenal <- ans$loglik
}
else{
fit$logLik <- ans$loglik
}
fit$alpha<-ans$b[np-nvar-1]^2
fit$eta<-ans$b[np-nvar]^2
if (noVar1 == 1) {
fit$coef <- NULL
}
else{
fit$coef <- ans$b[(np - nvar + 1):np]
names(fit$coef) <- factor.names(colnames(X))
}
temp1 <- matrix(ans$H, nrow = np, ncol = np)
temp2 <- matrix(ans$HIH, nrow = np, ncol = np)
fit$varH <- temp1[(np - nvar - 1):np, (np - nvar - 1):np]
fit$varHIH <- temp2[(np - nvar - 1):np, (np - nvar - 1):np]
seH.alpha <- sqrt(((2 * (fit$alpha^0.5))^2) * diag(fit$varH)[1])
fit$alpha_p.value <- 1 - pnorm(fit$alpha/seH.alpha)
seH.eta <- sqrt(((2 * (fit$eta^0.5))^2) * diag(fit$varH)[2])
fit$eta_p.value <- 1 - pnorm(fit$eta/seH.eta)
# fit$formula <- formula(Terms)
fit$x <- cbind(ans$x1,ans$x2)
fit$lam <- if(typeof == 1) {array(c(ans$lam[seq(1,length(ans$lam),3)],ans$lam2[seq(1,length(ans$lam2),3)]), dim=c(nb.int,3,2))} else{array(c(ans$lam,ans$lam2), dim=c(size1,3,2))}
fit$xSu <- cbind(ans$xSu1,ans$xSu2)
fit$surv <- array(c(ans$surv,ans$surv2), dim=c(size2,3,2))
fit$type <- type
fit$n.strat <- uni.strat
fit$n.iter <- ans$ni
fit$typeof <- typeof
fit$nb.gl <- nb.gl
fit$noVar1 <- noVar1
median <- NULL
for (i in (1:fit$n.strat)) median[i] <- ifelse(typeof==0, minmin(fit$surv[,1,i],fit$x), minmin(fit$surv[,1,i],fit$xSu))
lower <- NULL
for (i in (1:fit$n.strat)) lower[i] <- ifelse(typeof==0, minmin(fit$surv[,3,i],fit$x), minmin(fit$surv[,3,i],fit$xSu))
upper <- NULL
for (i in (1:fit$n.strat)) upper[i] <- ifelse(typeof==0, minmin(fit$surv[,2,i],fit$x), minmin(fit$surv[,2,i],fit$xSu))
fit$median <- cbind(lower,median,upper)
if (typeof == 0){
fit$n.knots<-n.knots
if (uni.strat > 1) fit$kappa <- ans$k0
else fit$kappa <- ans$k0[1]
fit$DoF <- ans$ddl
fit$cross.Val<-cross.validation
fit$zi <- ans$zi
}
if(typeof == 1)fit$time <- ans$time
#AD:
fit$nbintervR <- nbintervR
fit$nvar <- nvar
fit$LCV <- ans$LCV[1]
fit$AIC <- ans$LCV[2]
fit$npar <- np
fit$nst <- nst
if (typeof == 0){
# fit$indic.Kappa2 <- indic.Kappa2
fit$n.knots.temp <- n.knots.temp
}
fit$indic.nb.int <- indic.nb.int
fit$istop <- ans$istop
fit$shape.weib <- ans$shape.weib
fit$scale.weib <- ans$scale.weib
fit$AG <- recurrentAG
fit$EPS <- ans$EPS
# if (Frailty){
msg <- "Problem in the estimation of the random effects (perhaps high number of events in some clusters)"
if (any(ans$martingale.res==0)){
fit$martingale.res <- msg
fit$frailty.pred.group <- msg
fit$frailty.pred.subgroup <- msg
fit$linear.pred <- msg
}
else{
fit$martingale.res <- ans$martingale.res
fit$frailty.pred.group <- ans$frailty.pred.group
nom1 <- paste("g_",c(1:ngg),sep="")
nom2 <- paste("sub_g",c(1:maxng))
frailty.pred.subgroup <- as.data.frame(matrix(round(ans$frailty.pred.subgroup,6),ncol=maxng))
rownames(frailty.pred.subgroup) <- nom1
colnames(frailty.pred.subgroup) <- nom2
for (i in 1:ngg) {
if (subgbyg[i] < max(subgbyg)) {
frailty.pred.subgroup[i,(subgbyg[i]+1):max(subgbyg)] <- "."
}
}
fit$frailty.pred.subgroup <- frailty.pred.subgroup
# fit$frailty.var.group <- ans$frailty.var.group
# frailty.var.subgroup <- as.data.frame(matrix(round(ans$frailty.var.subgroup,6),nc=maxng))
# rownames(frailty.var.subgroup) <- nom1
# colnames(frailty.var.subgroup) <- nom2
#
# if(sum(which(subgbyg < max(subgbyg)))>0)frailty.var.subgroup[which(subgbyg < max(subgbyg)),(subgbyg[which(subgbyg < max(subgbyg))]+1):max(subgbyg)] <- "."
# fit$frailty.var.subgroup <- frailty.var.subgroup
# fit$frailty.sd.group <- ans$frailty.sd.group
# frailty.sd.subgroup <- as.data.frame(matrix(round(ans$frailty.sd.subgroup,6),nc=maxng))
# rownames(frailty.sd.subgroup) <- nom1
# colnames(frailty.sd.subgroup) <- nom2
# if(sum(which(subgbyg < max(subgbyg)))>0)frailty.sd.subgroup[which(subgbyg < max(subgbyg)),(subgbyg[which(subgbyg < max(subgbyg))]+1):max(subgbyg)] <- "."
# fit$frailty.sd.subgroup <- frailty.sd.subgroup
fit$linear.pred <- ans$linear.pred[order(ordre)]
}
fit$subgbyg <- subgbyg
# }
if(ans$ier==2000)
stop("The cross validation procedure cannot be finished. Try to change
either the number of knots or the seed for kappa parameter")
#========================= Test de Wald pour nested
if(length(vec.factor) > 0){
Beta <- ans$b[(np - nvar + 1):np]
#VarBeta <- fit$varH[2:(nvar+1),2:(nvar+1)]
VarBeta <- fit$varH[3:dim(fit$varH)[1], 3:dim(fit$varH)[2]]
nfactor <- length(vec.factor)
p.wald <- rep(0,nfactor)
if (fit$istop == 1) fit$global_chisq <- waldtest(N=nvar,nfact=nfactor,place=ind.place,modality=occur,b=Beta,Varb=VarBeta)
else fit$global_chisq <- 0
fit$dof_chisq <- occur
fit$global_chisq.test <- 1
# Calcul de pvalue globale
for(i in 1:length(vec.factor)){
p.wald[i] <- signif(1 - pchisq(fit$global_chisq[i], occur[i]), 3)
}
fit$p.global_chisq <- p.wald
fit$names.factor <- vec.factor
}
else{
fit$global_chisq.test <- 0
}
if(!is.null(fit$coef)){
if(length(fit$coef)>1) seH <- sqrt(diag(fit$varH))[-c(1:2)]
else seH <- sqrt(fit$varH)[-c(1:2)]
fit$beta_p.value <- 1 - pchisq((fit$coef/seH)^2, 1)
}
if (length(Xlevels) >0) fit$Xlevels <- Xlevels
fit$contrasts <- contr.save
fit$formula <- formula
attr(fit,"joint")<-joint
attr(fit,"subcluster")<-TRUE
class(fit) <- "nestedPenal"
} # End NESTED MODEL
# ===================================================
# BEGIN JOINT NESTED MODEL
# ===================================================
if (length(subcluster) & joint){
#if (!Frailty) stop("For joint nested frailty models, 'Frailty' must be equal to 'TRUE' !")
if (missing(formula.terminalEvent)) stop ("For joint nested frailty model, 'formula.terminalEvent' is required !")
if (joint.clust != 3) stop("Argument is mispecified for joint nested frailty model. Please look at the GenfrailtyPenal documentation.")
if (logNormal) stop("Sorry but log normal distribution is not available for joint nested frailty models")
if (uni.strat > 1) stop("Sorry but stratification for joint nested frailty model is not allowed")
if (classofY == "Surv")
{
if (!recurrentAG)
{
tt1.death <- aggregate(tt1, by=list(subcluster), FUN=sum)[,2]
tt0.death <- rep(0,length(tt1.death))
clusterdc <- 0
lignedc0 <- 0
tempdc <- 0
}
else{
tt1.death <- aggregate(tt1, by=list(subcluster), FUN=function(x) x[length(x)])[,2] # Myriam ! modifie sort()
tt0.death <- rep(0,length(tt1.death))
clusterdc <- 0
lignedc0 <- 0
tempdc <- 0
}
#stop ("The counting process approach of Andersen and Gill with a calendar timescale for recurrent event times is not yet allowed")}
}
else stop("Interval-censored data are not allowed for joint nested model.")
Terms2 <- if(missing(data))
{
terms(formula.terminalEvent, special)
}else{
terms(formula.terminalEvent, special, data = data)
}
ord2 <- attr(Terms2, "order")
if (length(ord2) & any(ord2 != 1)) stop("Interaction terms are not valid for terminalEvent formula")
if (ind.terminal)
{
terminalEvent<-aggregate(terminal,by=list(subcluster),FUN=function(x) x[length(x)])[,2]
icdc00 <- 0
}else{
stop("Nested joint frailty model mispecified. Please look at the GenfrailtyPenal documentation.")
}
# TerminalEvent doit etre de la forme 0-1
if(!all(terminalEvent %in% c(2,1,0))) stop("'terminal' must contain a variable coded 0-1 and a non-factor variable")
m2 <- match.call(expand.dots = FALSE)
m2$formula.terminalEvent <- m2$n.knots <- m2$recurrentAG <- m2$cross.validation <- m2$jointGeneral<- m2$kappa <- m2$maxit <- m2$hazard <- m2$nb.int <- m2$RandDist <- m2$betaorder <- m2$betaknots <- m2$init.B <- m2$LIMparam <- m2$LIMlogl <- m2$LIMderiv <- m2$print.times <- m2$init.Theta <- m2$init.Alpha <- m2$Alpha <- m2$init.Ksi <- m2$Ksi <- m2$init.Eta <- m2$Eta <- m2$initialize <- m2$nb.gl <- m2$... <- NULL
m2$formula <- Terms2
m2[[1]] <- as.name("model.frame")
m2 <- eval(m2, sys.parent())
match.noNA <- dimnames(m2)[[1]]%in%dimnames(m)[[1]]
m2 <- m2[match.noNA, ,drop=FALSE]
newTerms2 <- Terms2
#if (!missing(formula.terminalEvent)) {
X2 <- model.matrix(newTerms2, m2)
lldc <- attr(newTerms2,"term.labels")
ind.placedc <- unique(attr(X2,"assign")[duplicated(attr(X2,"assign"))])#changement unique le 26/09/2014
vec.factordc <- NULL
vec.factordc <- c(vec.factordc,lldc[ind.placedc])
mat.factordc <- matrix(vec.factordc,ncol=1,nrow=length(vec.factordc))
# Fonction servant a recuperer les termes dependants de "as.factor" en ce qui concerne l'evenement terminal
vec.factordc <-apply(mat.factordc,MARGIN=1,FUN=function(x){
if (length(grep("factor",x))>0){
if(length(grep(":",x))>0){
if(grep('\\(',unlist(strsplit(x,split="")))[1]<grep(":",unlist(strsplit(x,split="")))[1] && length(grep('\\(',unlist(strsplit(x,split=""))))==1){
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- grep(":",unlist(strsplit(x,split="")))[1]-2
pos3 <- grep(":",unlist(strsplit(x,split="")))[1]
pos4 <- length(unlist(strsplit(x,split="")))
return(paste(substr(x,start=pos1,stop=pos2),substr(x,start=pos3,stop=pos4),sep=""))
}else if(grep("\\(",unlist(strsplit(x,split="")))[1]>grep(":",unlist(strsplit(x,split="")))[1] && length(grep('\\(',unlist(strsplit(x,split=""))))==1){
pos2 <- grep(":",unlist(strsplit(x,split="")))[1]
pos3 <- grep("\\(",unlist(strsplit(x,split="")))[1]+1
pos4 <- length(unlist(strsplit(x,split="")))-1
return(paste(substr(x,start=1,stop=pos2),substr(x,start=pos3,stop=pos4),sep=""))
}else{#both factors
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- grep(":",unlist(strsplit(x,split="")))[1]-2
pos3 <- grep("\\(",unlist(strsplit(x,split="")))[2]+1
pos4 <- length(unlist(strsplit(x,split="")))-1
return(paste(substr(x,start=pos1,stop=pos2),":",substr(x,start=pos3,stop=pos4),sep=""))
}
}else{
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- length(unlist(strsplit(x,split="")))-1
return(substr(x,start=pos1,stop=pos2))
}
}else{
return(x)
}
}
)
# On determine le nombre de categories pour chaque var categorielle
if(length(vec.factordc) > 0)
{
vect.factdc <- attr(X2,"dimnames")[[2]]
vect.factdc <- vect.factdc[grep(paste(vec.factordc,collapse="|"),vect.factdc)]
occurdc <- rep(0,length(vec.factordc))
interactiondc<-as.vector(apply(matrix(vect.factdc,nrow=length(vect.factdc)),MARGIN=1,FUN=function(x){length(grep(":",unlist(strsplit(x,split=""))))}))
which.interactiondc <- which(interactiondc==1)
for(i in 1:length(vec.factordc))
{
if(length(grep(":",unlist(strsplit(vec.factordc[i],split=""))))>0)
{
pos <- grep(":",unlist(strsplit(vec.factordc[i],split="")))
length.grep <- 0
for(j in 1:length(vect.factdc))
{
if(j%in%which.interactiondc)
{
if(length(grep(substr(vec.factordc[i],start=1,stop=pos-1),vect.factdc[j]))>0 && length(grep(substr(vec.factordc[i],start=pos+1,stop=length(unlist(strsplit(vec.factordc[i],split="")))),vect.factdc[j]))>0)
{
length.grep <- length.grep + 1
which <- i
}
}
}
occurdc[i] <- length.grep
}else{
if(length(vect.factdc[-which.interactiondc])>0) occurdc[i] <- length(grep(vec.factordc[i],vect.factdc[-which.interactiondc]))
else occurdc[i] <- length(grep(vec.factordc[i],vect.factdc))
}
}
}
#=========================================================>
assign <- lapply(attrassign(X2, newTerms2)[-1], function(x) x - 1)
Xlevels2 <- .getXlevels(newTerms2, m2)
contr.save2 <- attr(X2, 'contrasts')
#========================================>
if(length(vec.factordc) > 0){
positiondc <- unlist(assign,use.names=F)
}
#========================================>
if (ncol(X2) == 1)
{
X2<-X2-1
noVar2 <- 1
}else{
X2 <- X2[, -1, drop = FALSE]
noVar2 <- 0
}
nvar2 <- ncol(X2)
if ((!is.null(attr(Terms2, "specials")$timedep))|timedep) stop("Time-dependant covariates are not allowed for joint nested frailty model.")
#vartimedep2 <- if (!is.null(attr(Terms2, "specials")$timedep)) stop("Nested frailty joint model not yet able to account time-dependant covariate")
# if (!is.null(vartimedep2)) timedep <- 1
#varnotdep2 <- colnames(X2)[-grep("timedep", colnames(X2))]
vardep2 <- colnames(X2)[grep("timedep", colnames(X2))]
vardep2 <- apply(matrix(vardep2, ncol=1, nrow=length(vardep2)), 1, timedep.names)
nvartimedep2 <- length(vardep2)
filtretps2 <- rep(0, nvar2)
#filtretps2[grep("timedep", colnames(X2))] <- 1
vardc.temp<-matrix(c(X2),nrow=nrow(X2),ncol=nvar2)
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Traitement des valeurs manquantes
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(is.null(nrow(m2)))
{
if (length(m2) != nrow(m)){
stop("There are missing values in the covariates modelling the terminal event. \n Prepare data only with complete cases")
}
}else{
if (nrow(m2) != nrow(m)){
stop(" There are missing values in the covariates modelling the terminal event. \n Prepare data only with complete cases")
}
}
if (!is.null(ncol(vardc.temp)))
{
vardc<-aggregate(vardc.temp[,1],by=list(sort(subcluster)), FUN=function(x) x[length(x)])[,2]
if (ncol(vardc.temp)>1)
{
for (i in 2:ncol(vardc.temp))
{
vardc.i<-aggregate(vardc.temp[,i],by=list(sort(subcluster)), FUN=function(x) x[length(x)])[,2]
vardc<-cbind(vardc,vardc.i)
}
}
}else vardc<-aggregate(vardc.temp,by=list(sort(subcluster)), FUN=function(x) x[length(x)])[,2]
vaxdc00 <- 0
#} # if ! missing formula.terminalEvent
nvarRec <- nvar
#=======================================>
#======= Construction du vecteur des indicatrice
if(length(vec.factordc) > 0){
k <- 0
for(i in 1:length(vec.factordc)){
ind.placedc[i] <- ind.placedc[i]+k
k <- k + occurdc[i]-1
}
}
if(is.null(nrow(vardc))) nvarEnd <- 1
else nvarEnd <- ncol(vardc)
if (sum(as.double(var)) == 0) nvarRec <- 0
if (sum(as.double(vardc)) == 0) nvarEnd <- 0
nvar <- nvarRec + nvarEnd
# ... End preparing data
effet <- 2
indic_alpha <- 1
indic_xi <- 1
if (!missing(Alpha))
{
if (Alpha == "None") indic_alpha <- 0
else stop("Alpha can only take 'None' as a value in this version of frailtypack package ")
}
if (!missing(Ksi))
{
if (Ksi == "None") indic_xi <- 0
else stop("Ksi can only take 'None' as a value in this version of frailtypack package")
}
indices <- c(indic_alpha, indic_xi)
nst <- uni.strat
if (typeof == 1) stop ("!Warning! Piecewise baseline hazard function is not yet allowed for joint nested frailty model.")
else {
indic.nb.int1 <- 0
indic.nb.int2 <- 0
}
np <- switch(as.character(typeof),
'0' = ((nst+1) * (n.knots + 2) + nvar + effet + indic_alpha + indic_xi), #splines
'2' = (2*(nst+1) + nvar + effet + indic_alpha + indic_xi)) # weibull
if (all(all.equal(as.numeric(cens), terminal) == T)) stop("'Recurrent event' variable and 'Terminal event' variable need to be different")
# Initialisation des parametres beta
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Beta <- rep(0,np)
if(!missing(init.B)){
if (length(init.B) != nvar) stop ("Wrong number of regression coefficient in init.B")
Beta <- c(rep(0, np-nvar), init.B)
}
if (!missing(init.Ksi)){
if(!missing(Ksi)) stop("You can't both initialize Ksi parameter and fit a joint model without Ksi in joint nested frailty model.")
if(!is.numeric(init.Ksi)) stop ("'init.Ksi' must be numeric")
Beta[np-nvar] <- init.Ksi
}
if (!missing(init.Alpha)){
if(!missing(Alpha)) stop("You can't both initialize alpha parameter and fit a joint model without it")
if(!is.numeric(init.Alpha)) stop ("'init.Alpha' must be numeric")
Beta[np-nvar-indic_xi] <- init.Alpha
}
if (!missing(init.Eta)){
if (!is.numeric(init.Eta)) stop("'init.Eta' must be numeric")
Beta[np-nvar-indic_xi-indic_alpha] <- init.Eta
}
if (!missing(init.Theta)){
if (!is.numeric(init.Theta)) stop("'init.Theta' must be numeric")
Beta[np-nvar-indic_xi-indic_alpha-1] <- init.Theta
}
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Clusterfamdc pour les deces : doit avoir autant de valeur que d'individus :
clusterfamdc <- aggregate(cluster, by=list(subcluster), FUN = function(x) x[length(x)])[,2]
xSu1 <- rep(0,100)
xSu2 <- rep(0,100)
if(typeof==0){
mt11 <- size1
mt12 <- size2
}else{
mt11 <- 100
mt12 <- 100
}
if ((typeof == 0) & (length(kappa) != (uni.strat+1))) stop ("Wrong length of argument kappa")
weights.vec <- rep(1, nrow(data))
weights.agg <- aggregate(weights.vec, by=list(cluster), FUN=function(x) x[length(x)])[,2]
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Fin de preparation des arguments
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
flush.console()
if (print.times){
ptm<-proc.time()
cat("\n")
cat("Be patient. The program is computing ... \n")
}
ans <- .Fortran(C_joint,
as.integer(n),
as.integer(c(length(uni.subcluster),length(uni.cluster),uni.strat)),
as.integer(strats),
as.integer(lignedc0),
as.integer(n.knots),
#k0=as.double(kappa)
axT = as.double(kappa),
as.double(tt0),
as.double(tt1),
as.integer(cens),
as.integer(subcluster),
as.integer(clusterdc),
as.integer(clusterfamdc),
as.double(tt0.death),
as.double(tt1.death),
as.integer(terminalEvent),
as.integer(tempdc),
as.integer(icdc00),
as.integer(nvarRec),
as.double(var),
as.integer(nvarEnd),
as.double(vardc),
as.double(vaxdc00),
as.integer(c(noVar1, noVar2)),
as.double(weights.agg),
#as.integer(noVar2),
as.integer(maxit),
np = as.integer(np),
b = as.double(Beta),
H_hessOut = as.double(matrix(0, nrow = np, ncol = np)),
HIHOut = as.double(matrix(0, nrow=np, ncol=np)),
loglik = as.double(0),
LCV = as.double(rep(0,2)),
xR = as.double(matrix(0, nrow=size1, ncol=uni.strat)),
lamR = as.double(array(0, dim=c(size1, 3, uni.strat))),
xSuR = as.double(xSu1),
survR = as.double(array(0, dim=c(mt11, 3, uni.strat))),
xD = as.double(rep(0, size2)),
lamD = as.double(matrix(0, nrow=size2, ncol=3)),
xSuD = as.double(xSu2),
survD = as.double(matrix(0, nrow=mt12, ncol=3)),
as.integer(typeof,equidistant),
as.integer(c(nbintervR,nbintervDC)),
#as.integer(nbintervDC),
as.integer(c(size1,size2, mt11, mt12)),
counts = as.integer(c(0,0,0,0)),
IerIstop = as.integer(c(0,0)),
#ier = as.integer(0),
#istop = as.integer(0),
paraweib = as.double(rep(0,4)),
MartinGales = as.double(matrix(0,nrow=length(uni.subcluster), ncol=5)),
linear.pred = as.double(rep(0,n)),
lineardc.pred = as.double(rep(0,as.integer(length(uni.subcluster)))),
zi = as.double(rep(0,(n.knots+6))),
time = as.double(rep(0,(nbintervR+1))),
timedc = as.double(rep(0,(nbintervDC+1))),
linearpredG = as.double(rep(0, lignedc0+1)),
typeJoint0 = as.integer(joint.clust),
as.integer(intcens),
as.integer(indices),
#as.integer(indic_alpha),
#as.integer(indic_xi),
as.double(ttU),
as.integer(ordretmp),
as.integer(initialize),
logNormal0 = as.integer(logNormal),
paratps = as.integer(c(timedep, betaknots, betaorder)),
as.integer(c(filtretps, filtretps2)),
BetaTpsMat = as.double(matrix(0,nrow=101, ncol=1+4*nvartimedep)),
BetaTpsMatDc = as.double(matrix(0,nrow=101, ncol=1+4*nvartimedep2)),
EPS = as.double(c(LIMparam, LIMlogl, LIMderiv)),
nbgauss = as.integer(c(nb.gh,nb.gl))
)#,
#PACKAGE = "frailtypack") #65 arguments
###########################################
### Verification de l'execution du code ###
###########################################
ier <- ans$IerIstop[1]
istop <- ans$IerIstop[2]
if (istop == 4) warning("Problem in the likelihood computation. The program stopped abnormally. Please verify your dataset.")
if (istop == 2) warning("Model did not converge.")
if (istop == 3) warning("Matrix non-positive definite.")
MartinGale <- matrix(ans$MartinGales, nrow=as.integer(length(uni.subcluster)), ncol=5)
fit <- NULL
fit$b <- ans$b
fit$timedep <- timedep
fit$typejoint0 <- joint.clust
fit$na.action <- attr(m,"na.action")
fit$call <- call
fit$n <- n
fit$groups <- length(uni.cluster)
fit$subgroups <- length(uni.subcluster)
fit$n.events <- ans$counts[2]
fit$n.death <- ans$counts[3]
fit$n.censored <- ans$counts[4]
fit$npar <- np
fit$nst <- nst
fit$LCV <- ans$LCV[1]
fit$AIC <- ans$LCV[2]
fit$nvar<-c(nvarRec,nvarEnd)
#fit$nvarnotdep<-c(nvarRec-nvartimedep,nvarEnd-nvartimedep2)
# fit$formula <- formula(Terms)
fit$xR <- matrix(ans$xR, nrow = size1, ncol = uni.strat)
fit$lamR <- array(ans$lamR, dim = c(size1,3,uni.strat))
fit$xSuR <- matrix(ans$xSuR, nrow = 100, ncol = uni.strat)
fit$survR <- array(ans$survR, dim = c(mt11,3,uni.strat))
fit$xD <- ans$xD
fit$lamD <- matrix(ans$lamD, nrow = size2, ncol = 3)
fit$xSuD <- ans$xSuD
fit$survD <- matrix(ans$survD, nrow = mt12, ncol = 3)
fit$type <- type
fit$n.strat <- uni.strat
fit$n.iter <- ans$counts[1]
fit$typeof <- typeof
fit$nb.gh <- nb.gh
fit$nb.gl <- nb.gl
medianR <- NULL
for (i in (1:fit$n.strat)) medianR[i] <- ifelse(typeof==0, minmin(fit$survR[,1,i],fit$xR), minmin(fit$survR[,1,i],fit$xSuR))
lowerR <- NULL
for (i in (1:fit$n.strat)) lowerR[i] <- ifelse(typeof==0, minmin(fit$survR[,2,i],fit$xR), minmin(fit$survR[,2,i],fit$xSuR))
upperR <- NULL
for (i in (1:fit$n.strat)) upperR[i] <- ifelse(typeof==0, minmin(fit$survR[,3,i],fit$xR), minmin(fit$survR[,3,i],fit$xSuR))
fit$medianR <- cbind(lowerR,medianR,upperR)
medianD <- ifelse(typeof==0, minmin(fit$survD[,1],fit$xD), minmin(fit$survD[,1],fit$xSuD))
lowerD <- ifelse(typeof==0, minmin(fit$survD[,3],fit$xD), minmin(fit$survD[,3],fit$xSuD))
upperD <- ifelse(typeof==0, minmin(fit$survD[,2],fit$xD), minmin(fit$survD[,2],fit$xSuD))
fit$medianD <- cbind(lowerD,medianD,upperD)
fit$noVar1 <- noVar1
fit$noVar2 <- noVar2
fit$nbintervR <- nbintervR
fit$nbintervDC <- nbintervDC
fit$nvarRec <- nvarRec
fit$nvarEnd <- nvarEnd
fit$istop <- istop
fit$indic.nb.intR <- indic.nb.int1
fit$indic.nb.intD <- indic.nb.int2
fit$shape.weib <- ans$paraweib[1:2]#ans$shape.weib
fit$scale.weib <- ans$paraweib[3:4]#ans$scale.weib
fit$joint.clust <- ans$typeJoint0
fit$AG <- recurrentAG
fit$intcens <- intcens
fit$indic_alpha <- indic_alpha
fit$indic_ksi <- indic_xi
fit$logNormal <- ans$logNormal0
#fit$BetaTpsMat <- matrix(ans$BetaTpsMat,nrow=101,ncol=1+4*nvartimedep) !! Pas encore MEP pour time-varying covariables
#fit$BetaTpsMatDc <- matrix(ans$BetaTpsMatDc,nrow=101,ncol=1+4*nvartimedep2)
#fit$nvartimedep <- c(nvartimedep,nvartimedep2)
fit$Names.vardep <- vardep
fit$Names.vardepdc <- vardep2
fit$EPS <- ans$EPS
if(typeof == 0) fit$logLikPenal <- ans$loglik
else fit$logLik <- ans$loglik
fit$theta <- ans$b[np-nvar-indic_xi-indic_alpha-1]^2
fit$eta <- ans$b[np-nvar-indic_alpha-indic_xi]^2
if (indic_alpha == 1) fit$alpha <- ans$b[np-nvar-indic_xi]
if (indic_xi == 1) fit$ksi <- ans$b[np-nvar]
if (noVar1==1 & noVar2==1) fit$coef <- NULL
else {
fit$coef <- ans$b[(np-nvar+1):np]
noms <- c(factor.names(colnames(X)),factor.names(colnames(X2)))
if(length(grep(":", noms))>0) noms <- factor.names(noms)
#if timedep == 1
names(fit$coef) <- noms
}
temp1 <- matrix(ans$H_hessOut, nrow = np, ncol = np)
temp2 <- matrix(ans$HIHOut, nrow = np, ncol = np)
fit$varHtotal <- temp1
fit$varHIHtotal <- temp2
fit$varH <- temp1[(np-nvar-indic_xi-indic_alpha-1):np, (np-nvar-indic_xi-indic_alpha-1):np]
fit$varHIH <- temp2[(np-nvar-indic_xi-indic_alpha-1):np, (np-nvar-indic_xi-indic_alpha-1):np]
seH.theta <- sqrt(((2 * (fit$theta^0.5))^2) * diag(fit$varH)[1])
fit$theta_p.value <- 1 - pnorm(fit$theta/seH.theta)
seH.eta <- sqrt(((2 * (fit$eta^0.5))^2) * diag(fit$varH)[2])
fit$eta_p.value <- 1 - pnorm(fit$eta/seH.eta)
if (indic_alpha == 1) fit$alpha_p.value <- 1 - pchisq((fit$alpha/sqrt(diag(fit$varH))[3])^2,1)
if (indic_xi == 1)fit$ksi_p.value <- 1 - pchisq((fit$ksi/sqrt(diag(fit$varH))[3+indic_alpha])^2,1)
if (indic_alpha == 1) noms <- c("theta","alpha",factor.names(colnames(X)),factor.names(colnames(X2)))
else noms <- c("theta",factor.names(colnames(X)),factor.names(colnames(X2)))
if(length(grep(":",noms))>0)noms <- factor.names(noms)
if (typeof == 0){
fit$n.knots<-n.knots
fit$kappa <- ans$axT
fit$cross.Val<-cross.validation
fit$n.knots.temp <- n.knots.temp
fit$zi <- ans$zi
}
# if constante par morceaux
# verif que les martingales ont ete bien calculees
#msg <- "Problem in the estimation of the random effects (perhaps high number of events in some clusters)"
#if (Frailty){
# if (any(MartinGale[,1]==0)){
# fit$martingale.res <- msg
# fit$martingaledeath.res <- msg
# fit$frailty.pred <- msg
# fit$linear.pred <- msg
# fit$lineardeath.pred <- msg
# }else{
# fit$martingale.res <- MartinGale[,1]
# fit$martingaledeath.res <- MartinGale[,2]
# fit$frailty.pred <- MartinGale[,3]
# fit$linear.pred <- ans$linearpred[order(ordre)]
# fit$lineardeath.pred <- ans$linearpreddc
# }
#}
fit$martingale.res <- MartinGale[,1] # seulement pour les familles
fit$martingaledeath.res <- MartinGale[,2]
fit$frailty.pred <- MartinGale[,3]
fit$frailty.fam.pred <- MartinGale[1:length(uni.cluster),5] ### AK 12/12/2016 (family frailties)
fit$linear.pred <- ans$linear.pred[order(ordre)]# a faire
fit$lineardeath.pred <- ans$lineardc.pred
#================================> For the Recurrent Event
#========================= Test de Wald
if (length(vec.factor) > 0){
Beta <- ans$b[(np - nvar + 1):np]
if (indic_alpha == 1) {
#if (indic_xi ==1) VarBeta <- diag(diag(fit$varH)[-c(1,2,3,4)])
#else VarBeta <- diag(diag(fit$varH)[-c(1,2,3)])
if(indic_xi == 1) VarBeta <- fit$varH[5:dim(fit$varH)[1],5:dim(fit$varH)[2]]
else VarBeta <- fit$varH[4:dim(fit$varH)[1],4:dim(fit$varH)[2]]
}else{
#if (indic_xi == 1) VarBeta <- diag(diag(fit$varH)[-c(1,2,3)])
#else VarBeta <- diag(diag(fit$varH)[-c(1,2)])
if(indic_xi == 1) VarBeta <- fit$varH[4:dim(fit$varH)[1],4:dim(fit$varH)[2]]
else VarBeta <- fit$varH[3:dim(fit$varH)[1],3:dim(fit$varH)[2]]
}
nfactor <- length(vec.factor)
p.wald <- rep(0,nfactor)
ntot <- nvarEnd + nvarRec
if(fit$istop == 1) fit$global_chisq <- waldtest(N=nvarRec,nfact=nfactor,place=ind.place,modality=occur,b=Beta,Varb=VarBeta,Llast=nvarEnd,Ntot=ntot)
else fit$global_chisq <- 0
fit$dof_chisq <- occur
fit$global_chisq.test <- 1
# Calcul de pvalue globale
for(i in 1:length(vec.factor)){
p.wald[i] <- signif(1 - pchisq(fit$global_chisq[i], occur[i]), 3)
}
fit$p.global_chisq <- p.wald
fit$names.factor <- vec.factor
}else{
fit$global_chisq.test <- 0
}
#================================> For the death
#========================= Test de Wald
if (length(vec.factordc) > 0){
Beta <- ans$b[(np-nvar+1):np]
if (indic_alpha == 1) {
if(indic_xi ==1) VarBeta <- diag(diag(fit$varH)[-c(1,2,3,4)])
else VarBeta <- diag(diag(fit$varH)[-c(1,2,3)])
}else{
if(indic_xi == 1) VarBeta <- diag(diag(fit$varH)[-c(1,2,3)])
else VarBeta <- diag(diag(fit$varH)[-c(1,2)])
}
nfactor <- length(vec.factordc)
p.walddc <- rep(0,nfactor)
ntot <- nvarEnd + nvarRec
if(fit$istop == 1) fit$global_chisq_d <- waldtest(N=nvarEnd,nfact=nfactor,place=ind.placedc,modality=occurdc,b=Beta,Varb=VarBeta,Lfirts=nvarRec,Ntot=ntot)
else fit$global_chisq_d <- 0
fit$dof_chisq_d <- occurdc
fit$global_chisq.test_d <- 1
# Calcul de pvalue globale
for(i in 1:length(vec.factordc)){
p.walddc[i] <- signif(1 - pchisq(fit$global_chisq_d[i], occurdc[i]), 3)
}
fit$p.global_chisq_d <- p.walddc
fit$names.factordc <- vec.factordc
}else fit$global_chisq.test_d <- 0
if(!is.null(fit$coef)){
if (indic_alpha == 1 & indic_xi == 1) seH <- sqrt(diag(fit$varH))[-c(1:4)]
else if (indic_alpha == 0 | indic_xi == 0) seH <- sqrt(diag(fit$varH))[-c(1:3)]
else if(indic_alpha == 0 & indic_xi == 0) seH <- sqrt(diag(fit$varH))[-c(1:2)]
fit$beta_p.value <- 1 - pchisq((fit$coef/seH)^2, 1)
}
if (length(Xlevels) >0) fit$Xlevels <- Xlevels
fit$contrasts <- contr.save
if (length(Xlevels2) >0) fit$Xlevels2 <- Xlevels2
fit$contrasts2 <- contr.save2
fit$formula <- formula
fit$formula.terminalEvent <- formula.terminalEvent
attr(fit,"joint")<-joint
attr(fit,"subcluster")<-TRUE
class(fit) <- "jointNestedPenal"
# ===================================================
# END NESTED JOINT MODEL
# ===================================================
}
if (print.times){
cost<-proc.time()-ptm
cat("The program took", round(cost[3],2), "seconds \n")
}
fit
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.