R/Jointlcmm.R
In CecileProust-Lima/lcmm: Extended Mixed Models Using Latent Classes and Latent Processes
Defines functions
Jointlcmm
Documented in
Jointlcmm
#' Estimation of joint latent class models for longitudinal and time-to-event
#' data
#'
#' This function fits joint latent class mixed models for a longitudinal
#' outcome and a right-censored (possibly left-truncated) time-to-event. The
#' function handles competing risks and Gaussian or non Gaussian (curvilinear)
#' longitudinal outcomes. For curvilinear longitudinal outcomes, normalizing
#' continuous functions (splines or Beta CDF) can be specified as in
#' \code{lcmm}.
#'
#'
#' A. BASELINE RISK FUNCTIONS
#'
#' For the baseline risk functions, the following parameterizations were
#' considered. Be careful, parametrisations changed in lcmm_V1.5:
#'
#' 1. With the "Weibull" function: 2 parameters are necessary w_1 and w_2 so
#' that the baseline risk function a_0(t) = w_1^2*w_2^2*(w_1^2*t)^(w_2^2-1) if
#' logscale=FALSE and a_0(t) = exp(w_1)*exp(w_2)(t)^(exp(w_2)-1) if
#' logscale=TRUE.
#'
#' 2. with the "piecewise" step function and nz nodes (y_1,...y_nz), nz-1
#' parameters are necesssary p_1,...p_nz-1 so that the baseline risk function
#' a_0(t) = p_j^2 for y_j < t =< y_j+1 if logscale=FALSE and a_0(t) = exp(p_j)
#' for y_j < t =< y_j+1 if logscale=TRUE.
#'
#' 3. with the "splines" function and nz nodes (y_1,...y_nz), nz+2 parameters
#' are necessary s_1,...s_nz+2 so that the baseline risk function a_0(t) =
#' sum_j s_j^2 M_j(t) if logscale=FALSE and a_0(t) = sum_j exp(s_j) M_j(t) if
#' logscale=TRUE where {M_j} is the basis of cubic M-splines.
#'
#' Two parametrizations of the baseline risk function are proposed
#' (logscale=TRUE or FALSE) because in some cases, especially when the
#' instantaneous risks are very close to 0, some convergence problems may
#' appear with one parameterization or the other. As a consequence, we
#' recommend to try the alternative parameterization (changing logscale option)
#' when a joint latent class model does not converge (maximum number of
#' iterations reached) where as convergence criteria based on the parameters
#' and likelihood are small.
#'
#' B. THE VECTOR OF PARAMETERS B
#'
#' The parameters in the vector of initial values \code{B} or in the vector of
#' maximum likelihood estimates \code{best} are included in the following
#' order: (1) ng-1 parameters are required for intercepts in the latent class
#' membership model, and if covariates are included in \code{classmb}, ng-1
#' parameters should be entered for each one; (2) parameters for the baseline
#' risk function: 2 parameters for each Weibull, nz-1 for each piecewise
#' constant risk and nz+2 for each splines risk; this number should be
#' multiplied by ng if specific hazard is specified; otherwise, ng-1 additional
#' proportional effects are expected if PH hazard is specified; otherwise
#' nothing is added if common hazard is specified. In the presence of competing
#' events, the number of parameters should be adapted to the number of causes
#' of event; (3) for all covariates in \code{survival}, ng parameters are
#' required if the covariate is inside a \code{mixture()}, otherwise 1
#' parameter is required. Covariates parameters should be included in the same
#' order as in \code{survival}. In the presence of cause-specific effects, the
#' number of parameters should be multiplied by the number of causes; (4) for
#' all covariates in \code{fixed}, one parameter is required if the covariate
#' is not in \code{mixture}, ng parameters are required if the covariate is
#' also in \code{mixture}. Parameters should be included in the same order as
#' in \code{fixed}; (5) the variance of each random-effect specified in
#' \code{random} (including the intercept) if \code{idiag=TRUE} and the
#' inferior triangular variance-covariance matrix of all the random-effects if
#' \code{idiag=FALSE}; (6) only if \code{nwg=TRUE}, ng-1 parameters for
#' class-specific proportional coefficients for the variance covariance matrix
#' of the random-effects; (7) the variance of the residual error.
#'
#' C. CAUTION
#'
#' Some caution should be made when using the program:
#'
#' (1) As the log-likelihood of a latent class model can have multiple maxima,
#' a careful choice of the initial values is crucial for ensuring convergence
#' toward the global maximum. The program can be run without entering the
#' vector of initial values (see point 2). However, we recommend to
#' systematically enter initial values in \code{B} and try different sets of
#' initial values.
#'
#' (2) The automatic choice of initial values that we provide requires the
#' estimation of a preliminary linear mixed model. The user should be aware
#' that first, this preliminary analysis can take time for large datatsets and
#' second, that the generated initial values can be very not likely and even
#' may converge slowly to a local maximum. This is a reason why several
#' alternatives exist. The vector of initial values can be directly specified
#' in \code{B} the initial values can be generated (automatically or randomly)
#' from a model with \code{ng=}. Finally, function \code{gridsearch} performs
#' an automatic grid search.
#'
#' (3) Convergence criteria are very strict as they are based on derivatives of
#' the log-likelihood in addition to the parameter and log-likelihood
#' stability. In some cases, the program may not converge and reach the
#' maximum number of iterations fixed at 150. In this case, the user should
#' check that parameter estimates at the last iteration are not on the
#' boundaries of the parameter space. If the parameters are on the boundaries
#' of the parameter space, the identifiability of the model is critical. This
#' may happen especially when baseline risk functions involve splines (value
#' close to the lower boundary - 0 with logscale=F -infinity with logscale=F)
#' or classmb parameters that are too high or low (perfect classification) or
#' linkfunction parameters. When identifiability of some parameters is
#' suspected, the program can be run again from the former estimates by fixing
#' the suspected parameters to their value with option posfix. This usually
#' solves the problem. An alternative is to remove the parameters of the Beta
#' of Splines link function from the inverse of the Hessian with option
#' partialH. If not, the program should be run again with other initial
#' values. Some problems of convergence may happen when the instantaneous
#' risks of event are very low and "piecewise" or "splines" baseline risk
#' functions are specified. In this case, changing the parameterization of the
#' baseline risk functions with option logscale is recommended (see paragraph A
#' for details).
#'
#' @aliases Jointlcmm jlcmm
#' @param fixed two-sided linear formula object for the fixed-effects in the
#' linear mixed model. The response outcome is on the left of \code{~} and the
#' covariates are separated by \code{+} on the right of the \code{~}. By
#' default, an intercept is included. If no intercept, \code{-1} should be the
#' first term included on the right of \code{~}.
#' @param mixture one-sided formula object for the class-specific fixed effects
#' in the linear mixed model (to specify only for a number of latent classes
#' greater than 1). Among the list of covariates included in \code{fixed}, the
#' covariates with class-specific regression parameters are entered in
#' \code{mixture} separated by \code{+}. By default, an intercept is included.
#' If no intercept, \code{-1} should be the first term included.
#' @param random optional one-sided formula for the random-effects in the
#' linear mixed model. Covariates with a random-effect are separated by
#' \code{+}. By default, an intercept is included. If no intercept, \code{-1}
#' should be the first term included.
#' @param subject name of the covariate representing the grouping structure
#' (called subject identifier) specified with ''.
#' @param classmb optional one-sided formula describing the covariates in the
#' class-membership multinomial logistic model. Covariates included are
#' separated by \code{+}. No intercept should be included in this formula.
#' @param ng optional number of latent classes considered. If \code{ng=1} (by
#' default) no \code{mixture} nor \code{classmb} should be specified. If
#' \code{ng>1}, \code{mixture} is required.
#' @param idiag optional logical for the structure of the variance-covariance
#' matrix of the random-effects. If \code{FALSE}, a non structured matrix of
#' variance-covariance is considered (by default). If \code{TRUE} a diagonal
#' matrix of variance-covariance is considered.
#' @param nwg optional logical indicating if the variance-covariance of the
#' random-effects is class-specific. If \code{FALSE} the variance-covariance
#' matrix is common over latent classes (by default). If \code{TRUE} a
#' class-specific proportional parameter multiplies the variance-covariance
#' matrix in each class (the proportional parameter in the last latent class
#' equals 1 to ensure identifiability).
#' @param survival two-sided formula object. The left side of the formula
#' corresponds to a \code{surv()} object of type "counting" for right-censored
#' and left-truncated data (example: \code{Surv(Time,EntryTime,Indicator)}) or
#' of type "right" for right-censored data (example:
#' \code{Surv(Time,Indicator)}). Multiple causes of event can be considered in
#' the Indicator (0 for censored, k for cause k of event). The right side of
#' the formula specifies the names of covariates to include in the survival
#' model with \code{mixture()} when the effect is class-specific (example:
#' \code{Surv(Time,Indicator) ~} \code{ X1 + mixture(X2)} for a class-common
#' effect of X1 and a class-specific effect of X2). In the presence of
#' competing events, covariate effects are common by default. Code
#' \code{cause(X3)} specifies a cause-specific covariate effect for X3 on each
#' cause of event while \code{cause1(X3)} (or \code{cause2(X3)}, ...) specifies
#' a cause-specific effect of X3 on the first (or second, ...) cause only.
#' @param hazard optional family of hazard function assumed for the survival
#' model. By default, "Weibull" specifies a Weibull baseline risk function.
#' Other possibilities are "piecewise" for a piecewise constant risk function
#' or "splines" for a cubic M-splines baseline risk function. For these two
#' latter families, the number of nodes and the location of the nodes should be
#' specified as well, separated by \code{-}. The number of nodes is entered
#' first followed by \code{-}, then the location is specified with "equi",
#' "quant" or "manual" for respectively equidistant nodes, nodes at quantiles
#' of the times of event distribution or interior nodes entered manually in
#' argument \code{hazardnodes}. It is followed by \code{-} and finally
#' "piecewise" or "splines" indicates the family of baseline risk function
#' considered. Examples include "5-equi-splines" for M-splines with 5
#' equidistant nodes, "6-quant-piecewise" for piecewise constant risk over 5
#' intervals and nodes defined at the quantiles of the times of events
#' distribution and "9-manual-splines" for M-splines risk function with 9
#' nodes, the vector of 7 interior nodes being entered in the argument
#' \code{hazardnodes}. In the presence of competing events, a vector of hazards
#' should be provided such as \code{hazard=c("Weibull","splines"} with 2 causes
#' of event, the first one modelled by a Weibull baseline cause-specific risk
#' function and the second one by splines.
#' @param hazardtype optional indicator for the type of baseline risk function
#' when ng>1. By default "Specific" indicates a class-specific baseline risk
#' function. Other possibilities are "PH" for a baseline risk function
#' proportional in each latent class, and "Common" for a baseline risk function
#' that is common over classes. In the presence of competing events, a vector
#' of hazardtypes should be given.
#' @param hazardnodes optional vector containing interior nodes if
#' \code{splines} or \code{piecewise} is specified for the baseline hazard
#' function in \code{hazard}.
#' @param hazardrange optional vector indicating the range of the survival times (that is
#' the minimum and maximum). By default, the range is defined according to the
#' minimum and maximum observed values of the survival times. The option should be
#' used only for piecewise constant and Splines hazard functions.
#' @param TimeDepVar optional vector containing an intermediate time
#' corresponding to a change in the risk of event. This time-dependent
#' covariate can only take the form of a time variable with the assumption that
#' there is no effect on the risk before this time and a constant effect on the
#' risk of event after this time (example: initiation of a treatment to account
#' for).
#' @param link optional family of link functions to estimate. By default,
#' "linear" option specifies a linear link function leading to a standard
#' linear mixed model (homogeneous or heterogeneous as estimated in
#' \code{hlme}). Other possibilities include "beta" for estimating a link
#' function from the family of Beta cumulative distribution functions,
#' "thresholds" for using a threshold model to describe the correspondence
#' between each level of an ordinal outcome and the underlying latent process,
#' and "Splines" for approximating the link function by I-splines. For this
#' latter case, the number of nodes and the nodes location should be also
#' specified. The number of nodes is first entered followed by \code{-}, then
#' the location is specified with "equi", "quant" or "manual" for respectively
#' equidistant nodes, nodes at quantiles of the marker distribution or interior
#' nodes entered manually in argument \code{intnodes}. It is followed by
#' \code{-} and finally "splines" is indicated. For example, "7-equi-splines"
#' means I-splines with 7 equidistant nodes, "6-quant-splines" means I-splines
#' with 6 nodes located at the quantiles of the marker distribution and
#' "9-manual-splines" means I-splines with 9 nodes, the vector of 7 interior
#' nodes being entered in the argument \code{intnodes}.
#' @param intnodes optional vector of interior nodes. This argument is only
#' required for a I-splines link function with nodes entered manually.
#' @param epsY optional definite positive real used to rescale the marker in
#' (0,1) when the beta link function is used. By default, epsY=0.5.
#' @param range optional vector indicating the range of the outcome (that is
#' the minimum and maximum). By default, the range is defined according to the
#' minimum and maximum observed values of the outcome. The option should be
#' used only for Beta and Splines transformations.
#' @param cor optional brownian motion or autoregressive process modeling the
#' correlation between the observations. "BM" or "AR" should be specified,
#' followed by the time variable between brackets. By default, no correlation
#' is added.
#' @param data optional data frame containing the variables named in
#' \code{fixed}, \code{mixture}, \code{random}, \code{classmb} and
#' \code{subject}.
#' @param B optional specification for the initial values for the parameters.
#' Three options are allowed: (1) a vector of initial values is entered (the
#' order in which the parameters are included is detailed in \code{details}
#' section). (2) nothing is specified. A preliminary analysis involving the
#' estimation of a standard linear mixed model is performed to choose initial
#' values. (3) when ng>1, a Jointlcmm object is entered. It should correspond
#' to the exact same structure of model but with ng=1. The program will
#' automatically generate initial values from this model. This specification
#' avoids the preliminary analysis indicated in (2) Note that due to possible
#' local maxima, the \code{B} vector should be specified and several different
#' starting points should be tried.
#' @param convB optional threshold for the convergence criterion based on the
#' parameter stability. By default, convB=0.0001.
#' @param convL optional threshold for the convergence criterion based on the
#' log-likelihood stability. By default, convL=0.0001.
#' @param convG optional threshold for the convergence criterion based on the
#' derivatives. By default, convG=0.0001.
#' @param maxiter optional maximum number of iterations for the Marquardt
#' iterative algorithm. By default, maxiter=150.
#' @param nsim optional number of points for the predicted survival curves and
#' predicted baseline risk curves. By default, nsim=100.
#' @param prior optional name of a covariate containing a prior information
#' about the latent class membership. The covariate should be an integer with
#' values in 0,1,...,ng. Value O indicates no prior for the subject while a
#' value in 1,...,ng indicates that the subject belongs to the corresponding
#' latent class.
#' @param pprior optional vector specifying the names of the covariates containing the
#' prior probabilities to belong to each latent class. These probabilities should be
#' between 0 and 1 and should sum up to 1 for each subject.
#' @param logscale optional boolean indicating whether an exponential
#' (logscale=TRUE) or a square (logscale=FALSE -by default) transformation is
#' used to ensure positivity of parameters in the baseline risk functions. See
#' details section
#' @param subset a specification of the rows to be used: defaults to all rows.
#' This can be any valid indexing vector for the rows of data or if that is not
#' supplied, a data frame made up of the variable used in formula.
#' @param na.action Integer indicating how NAs are managed. The default is 1
#' for 'na.omit'. The alternative is 2 for 'na.fail'. Other options such as
#' 'na.pass' or 'na.exclude' are not implemented in the current version.
#' @param posfix Optional vector specifying the indices in vector B of the
#' parameters that should not be estimated. Default to NULL, all parameters are
#' estimated.
#' @param partialH optional logical for Piecewise and Splines baseline risk
#' functions and Splines link functions only. Indicates whether the parameters of the
#' baseline risk or link functions can be dropped from the Hessian matrix to
#' define convergence criteria.
#' @param verbose logical indicating if information about computation should be
#' reported. Default to TRUE.
#' @param returndata logical indicating if data used for computation should be
#' returned. Default to FALSE, data are not returned.
#' @param var.time optional character indicating the name of the time variable.
#' @param nproc the number cores for parallel computation.
#' Default to 1 (sequential mode).
#' @param clustertype optional character indicating the type of cluster for parallel computation.
#' @return The list returned is: \item{loglik}{log-likelihood of the model}
#' \item{best}{vector of parameter estimates in the same order as specified in
#' \code{B} and detailed in section \code{details}}
#' \item{V}{if the model converged (conv=1 or 3), vector containing
#' the upper triangle matrix of variance-covariance estimates of \code{Best}
#' with exception for variance-covariance parameters of the random-effects for
#' which \code{V} contains the variance-covariance estimates of the Cholesky
#' transformed parameters displayed in \code{cholesky}.
#' If conv=2, \code{V} contains the second derivatives of the log-likelihood.}
#' \item{gconv}{vector of
#' convergence criteria: 1. on the parameters, 2. on the likelihood, 3. on the
#' derivatives} \item{conv}{status of convergence: =1 if the convergence
#' criteria were satisfied, =2 if the maximum number of iterations was reached,
#' =4 or 5 if a problem occured during optimisation} \item{call}{the matched
#' call} \item{niter}{number of Marquardt iterations} \item{pred}{table of
#' individual predictions and residuals; it includes marginal predictions
#' (pred_m), marginal residuals (resid_m), subject-specific predictions
#' (pred_ss) and subject-specific residuals (resid_ss) averaged over classes,
#' the observation (obs) and finally the class-specific marginal and
#' subject-specific predictions (with the number of the latent class:
#' pred_m_1,pred_m_2,...,pred_ss_1,pred_ss_2,...). If \code{var.time}
#' is specified, the corresponding measurement time is also included.}
#' \item{pprob}{table of
#' posterior classification and posterior individual class-membership
#' probabilities based on the longitudinal data and the time-to-event data}
#' \item{pprobY}{table of posterior classification and posterior individual
#' class-membership probabilities based only on the longitudinal data}
#' \item{predRE}{table containing individual predictions of the random-effects:
#' a column per random-effect, a line per subject} \item{cholesky}{vector
#' containing the estimates of the Cholesky transformed parameters of the
#' variance-covariance matrix of the random-effects} \item{scoretest}{Statistic
#' of the Score Test for the conditional independence assumption of the
#' longitudinal and survival data given the latent class structure. Under the
#' null hypothesis, the statistics is a Chi-square with p degrees of freedom
#' where p indicates the number of random-effects in the longitudinal mixed
#' model. See Jacqmin-Gadda and Proust-Lima (2009) for more details.}
#' \item{predSurv}{table of predictions giving for the window of times to event
#' (called "time"), the predicted baseline risk function in each latent class
#' (called "RiskFct") and the predicted cumulative baseline risk function in
#' each latent class (called "CumRiskFct").} \item{hazard}{internal information
#' about the hazard specification used in related functions} \item{data}{the
#' original data set (if returndata is TRUE)}
#' %\item{specif}{internal information used in related functions}
#' %\item{Names}{internal information used in related fnctions}
#' %\item{Names2}{internal information used in related functions}
#' @author Cecile Proust Lima, Amadou Diakite and Viviane Philipps
#'
#' \email{cecile.proust-lima@@inserm.fr}
#' @seealso \code{\link{postprob}}, \code{\link{plot.Jointlcmm}},
#' \code{\link{plot.predict}}, \code{\link{epoce}}
#' @references
#'
#' Proust-Lima C, Philipps V, Liquet B (2017). Estimation of Extended Mixed
#' Models Using Latent Classes and Latent Processes: The R Package lcmm.
#' Journal of Statistical Software, 78(2), 1-56. doi:10.18637/jss.v078.i02
#'
#'
#' Lin, H., Turnbull, B. W., McCulloch, C. E. and Slate, E. H. (2002). Latent
#' class models for joint analysis of longitudinal biomarker and event process
#' data: application to longitudinal prostate-specific antigen readings and
#' prostate cancer. Journal of the American Statistical Association 97, 53-65.
#'
#' Proust-Lima, C. and Taylor, J. (2009). Development and validation of a
#' dynamic prognostic tool for prostate cancer recurrence using repeated
#' measures of post-treatment PSA: a joint modelling approach. Biostatistics
#' 10, 535-49.
#'
#' Jacqmin-Gadda, H. and Proust-Lima, C. (2010). Score test for conditional
#' independence between longitudinal outcome and time-to-event given the
#' classes in the joint latent class model. Biometrics 66(1), 11-9
#'
#' Proust-Lima, Sene, Taylor and Jacqmin-Gadda (2014). Joint latent class
#' models of longitudinal and time-to-event data: a review. Statistical Methods
#' in Medical Research 23, 74-90.
#' @examples
#'
#' \dontrun{
#' #### Example of a joint latent class model estimated for a varying number
#' # of latent classes:
#' # The linear mixed model includes a subject- (ID) and class-specific
#' # linear trend (intercept and Time in fixed, random and mixture components)
#' # and a common effect of X1 and its interaction with time over classes
#' # (in fixed).
#' # The variance of the random intercept and slopes are assumed to be equal
#' # over classes (nwg=F).
#' # The covariate X3 predicts the class membership (in classmb).
#' # The baseline hazard function is modelled with cubic M-splines -3
#' # nodes at the quantiles- (in hazard) and a proportional hazard over
#' # classes is assumed (in hazardtype). Covariates X1 and X2 predict the
#' # risk of event (in survival) with a common effect over classes for X1
#' # and a class-specific effect of X2.
#' # !CAUTION: for illustration, only default initial values where used but
#' # other sets of initial values should be tried to ensure convergence
#' # towards the global maximum.
#'
#'
#' #### estimation with 1 latent class (ng=1): independent models for the
#' # longitudinal outcome and the time of event
#' m1 <- Jointlcmm(fixed= Ydep1~X1*Time,random=~Time,subject='ID',
#' survival = Surv(Tevent,Event)~ X1+X2 ,hazard="3-quant-splines",
#' hazardtype="PH",ng=1,data=data_lcmm)
#' summary(m1)
#' #Goodness-of-fit statistics for m1:
#' # maximum log-likelihood: -3944.77 ; AIC: 7919.54 ; BIC: 7975.09
#' }
#'
#' #### estimation with 2 latent classes (ng=2)
#' m2 <- Jointlcmm(fixed= Ydep1~Time*X1,mixture=~Time,random=~Time,
#' classmb=~X3,subject='ID',survival = Surv(Tevent,Event)~X1+mixture(X2),
#' hazard="3-quant-splines",hazardtype="PH",ng=2,data=data_lcmm,
#' B=c(0.64,-0.62,0,0,0.52,0.81,0.41,0.78,0.1,0.77,-0.05,10.43,11.3,-2.6,
#' -0.52,1.41,-0.05,0.91,0.05,0.21,1.5))
#' summary(m2)
#' #Goodness-of-fit statistics for m2:
#' # maximum log-likelihood: -3921.27; AIC: 7884.54; BIC: 7962.32
#'
#' \dontrun{
#' #### estimation with 3 latent classes (ng=3)
#' m3 <- Jointlcmm(fixed= Ydep1~Time*X1,mixture=~Time,random=~Time,
#' classmb=~X3,subject='ID',survival = Surv(Tevent,Event)~ X1+mixture(X2),
#' hazard="3-quant-splines",hazardtype="PH",ng=3,data=data_lcmm,
#' B=c(0.77,0.4,-0.82,-0.27,0,0,0,0.3,0.62,2.62,5.31,-0.03,1.36,0.82,
#' -13.5,10.17,10.24,11.51,-2.62,-0.43,-0.61,1.47,-0.04,0.85,0.04,0.26,1.5))
#' summary(m3)
#' #Goodness-of-fit statistics for m3:
#' # maximum log-likelihood: -3890.26 ; AIC: 7834.53; BIC: 7934.53
#'
#' #### estimation with 4 latent classes (ng=4)
#' m4 <- Jointlcmm(fixed= Ydep1~Time*X1,mixture=~Time,random=~Time,
#' classmb=~X3,subject='ID',survival = Surv(Tevent,Event)~ X1+mixture(X2),
#' hazard="3-quant-splines",hazardtype="PH",ng=4,data=data_lcmm,
#' B=c(0.54,-0.42,0.36,-0.94,-0.64,-0.28,0,0,0,0.34,0.59,2.6,2.56,5.26,
#' -0.1,1.27,1.34,0.7,-5.72,10.54,9.02,10.2,11.58,-2.47,-2.78,-0.28,-0.57,
#' 1.48,-0.06,0.61,-0.07,0.31,1.5))
#' summary(m4)
#' #Goodness-of-fit statistics for m4:
#' # maximum log-likelihood: -3886.93 ; AIC: 7839.86; BIC: 7962.09
#'
#'
#' ##### The model with 3 latent classes is retained according to the BIC
#' ##### and the conditional independence assumption is not rejected at
#' ##### the 5% level.
#' # posterior classification
#' plot(m3,which="postprob")
#' # Class-specific predicted baseline risk & survival functions in the
#' # 3-class model retained (for the reference value of the covariates)
#' plot(m3,which="baselinerisk",bty="l")
#' plot(m3,which="baselinerisk",ylim=c(0,5),bty="l")
#' plot(m3,which="survival",bty="l")
#' # class-specific predicted trajectories in the 3-class model retained
#' # (with characteristics of subject ID=193)
#' data <- data_lcmm[data_lcmm$ID==193,]
#' plot(predictY(m3,var.time="Time",newdata=data,bty="l"))
#' # predictive accuracy of the model evaluated with EPOCE
#' vect <- 1:15
#' cvpl <- epoce(m3,var.time="Time",pred.times=vect)
#' summary(cvpl)
#' plot(cvpl,bty="l",ylim=c(0,2))
#' ############## end of example ##############
#' }
#'
#' @export
#'
#'
Jointlcmm <- function(fixed,mixture,random,subject,classmb,ng=1,idiag=FALSE,nwg=FALSE,
survival,hazard="Weibull",hazardtype="Specific",
hazardnodes=NULL,hazardrange=NULL,TimeDepVar=NULL,
link=NULL,intnodes=NULL,epsY=0.5,range=NULL,
cor=NULL,data,B,convB=0.0001,convL=0.0001,convG=0.0001,maxiter=100,
nsim=100,prior=NULL,pprior=NULL,logscale=FALSE,subset=NULL,na.action=1,
posfix=NULL,partialH=FALSE,verbose=FALSE,returndata=FALSE,
var.time=NULL,nproc=1,clustertype=NULL)
{
ptm<-proc.time()
cl <- match.call()
if(!missing(mixture) & ng==1) stop("No mixture can be specified with ng=1")
if(missing(mixture) & ng>1) stop("The argument mixture has to be specified for ng > 1")
if(!missing(classmb) & ng==1) stop("No classmb can be specified with ng=1")
if(missing(fixed)) stop("The argument Fixed must be specified in any model")
if(missing(random)) random <- ~-1
if(missing(classmb) & ng==1) classmb <- ~-1
if(missing(classmb) & ng>1) classmb <- ~1
if(missing(mixture)) mixture <- ~-1
if(ng==1 & nwg==TRUE) stop ("The argument nwg should be FALSE for ng=1")
if(ng==1) hazardtype <- "Specific" #??
if(any(!(hazardtype %in% c("Common","Specific","PH")))) stop("'hazardtype' should be either 'Common' or 'Specific' or 'PH'")
if(!inherits(fixed,"formula")) stop("The argument fixed must be a formula")
if(!inherits(mixture,"formula")) stop("The argument mixture must be a formula")
if(!inherits(random,"formula")) stop("The argument random must be a formula")
if(!inherits(classmb,"formula")) stop("The argument classmb must be a formula")
if(missing(data)){ stop("The argument data should be specified and defined as a data.frame")}
if(nrow(data)==0) stop("Data should not be empty")
if(missing(subject)){ stop("The argument subject must be specified in any model even without random-effects")}
if(!is.numeric(data[[subject]])) stop("The argument subject must be numeric")
if(is.null(link)) link <- "NULL"
if(any(link=="thresholds")) stop("The link function thresholds is not available yet")
if(link %in% c("linear","beta") & !is.null(intnodes)) stop("Intnodes should only be specified with splines links")
nom.subject <- as.character(subject)
if(!isTRUE(nom.subject %in% colnames(data))) stop(paste("Data should contain variable",nom.subject))
nom.prior <- NULL
if(!is.null(prior))
{
nom.prior <- as.character(prior)
if(!isTRUE(nom.prior %in% colnames(data))) stop(paste("Data should contain variable",nom.prior))
}
nom.pprior <- NULL
if(!is.null(pprior))
{
if(ng==1) stop("pprior is only useful if ng>1")
if(!is.character(pprior)) stop("pprior should be a character vector")
if(length(pprior) != ng) stop("pprior should be of length ng")
nom.pprior <- as.character(pprior)
if(!all(nom.pprior %in% colnames(data))) stop(paste("Data should contain variable",paste(nom.pprior, collapse=" ")))
}
nom.timedepvar <- NULL
if(!missing(TimeDepVar))
{
if(!is.null(TimeDepVar))
{
nom.timedepvar <- as.character(TimeDepVar)
if(!isTRUE(nom.timedepvar %in% colnames(data))) stop(paste("Data should contain variable",nom.timedepvar))
}
}
if(!(na.action %in% c(1,2))) stop("only 1 for 'na.omit' or 2 for 'na.fail' are required in na.action argument")
#if(length(posfix) & missing(B)) stop("A set of initial parameters must be specified if some parameters are not estimated")
## garder data tel quel pour le renvoyer
if(returndata==TRUE)
{
datareturn <- data
}
else
{
datareturn <- NULL
}
if(!isTRUE(all.equal(as.character(cl$subset),character(0))))
{
cc <- cl
cc <- cc[c(1,which(names(cl)=="subset"))]
cc[[1]] <- as.name("model.frame")
cc$formula <- formula(paste("~",paste(colnames(data),collapse="+")))
cc$data <- data
cc$na.action <- na.pass
data <- eval(cc)
}
attributes(data)$terms <- NULL
### test de l'argument cor
ncor0 <- 0
cor.type <- cl$cor[1]
cor.time <- cl$cor[2]
cor <- paste(cor.type,cor.time,sep="-")
if (!isTRUE(all.equal(cor,character(0))))
{
if (substr(cor,1,2)=="AR") { ncor0 <- 2 }
else if (substr(cor,1,2)=="BM") { ncor0 <- 1 }
else { stop("The argument cor must be of type AR or BM") }
if(!(strsplit(cor,"-")[[1]][2] %in% colnames(data))) stop("Unable to find time variable from argument 'cor' in 'data'")
else { cor.var.time <- strsplit(cor,"-")[[1]][2] }
}
### fin test argument cor
## objet Surv
surv <- cl$survival[[2]]
if(length(surv)==3) #censure droite sans troncature gauche
{
idtrunc <- 0
Tevent <- getElement(object=data,name=as.character(surv[2]))
Event <- getElement(object=data,name=as.character(surv[3]))
Tentry <- rep(0,length(Tevent)) #si pas de troncature, Tentry=0
noms.surv <- c(as.character(surv[2]),as.character(surv[3]))
surv <- do.call("Surv",list(time=Tevent,event=Event,type="mstate"))
}
else
{
if(length(surv)==4) #censure droite et troncature
{
idtrunc <- 1
Tentry <- getElement(object=data,name=as.character(surv[2]))
Tevent <- getElement(object=data,name=as.character(surv[3]))
Event <- getElement(object=data,name=as.character(surv[4]))
noms.surv <- c(as.character(surv[2]),as.character(surv[3]),as.character(surv[4]))
surv <- do.call("Surv",list(time=Tentry,time2=Tevent,event=Event,type="mstate"))
}
}
## nombre d'evenement concurrents
nbevt <- length(attr(surv,"states"))
if(nbevt<1)
{
if(maxiter != 0) stop("No observed event in the data")
nbevt <- 1
}
##pour acces aux attributs des formules
afixed <- terms(fixed)
if(link!="NULL" & attr(afixed,"intercept")==0) stop("An intercept should appear in fixed for identifiability purposes")
amixture <- terms(mixture)
arandom <- terms(random)
aclassmb <- terms(classmb)
## pour la formule pour survivial, creer 3 formules :
## une pour les covariables en mixture, une pour les covariables avec effet specifique a la cause, et une pour les effets communs.
form.surv <- cl$survival[3]
noms.form.surv <- all.vars(attr(terms(formula(paste("~",form.surv))),"variables"))
if(length(noms.form.surv)==0)
{
form.cause <- ~-1
form.causek <- vector("list",nbevt)
for(k in 1:nbevt) form.causek[[k]] <- ~-1
form.mixture <- ~-1
form.commun <- ~-1
asurv <- terms(~-1)
}
else
{
##creer la formula pour cause
form1 <- gsub("mixture","",form.surv)
form1 <- formula(paste("~",form1))
asurv1 <- terms(form1,specials="cause")
ind.cause <- attr(asurv1,"specials")$cause
if(length(ind.cause))
{
form.cause <- paste(labels(asurv1)[ind.cause],collapse="+")
form.cause <- gsub("cause","",form.cause)
form.cause <- formula(paste("~",form.cause))
}
else
{
form.cause <- ~-1
}
## formules pour causek
form.causek <- vector("list",nbevt)
for(k in 1:nbevt)
{
formk <- gsub("mixture","",form.surv)
for(kk in 1:nbevt)
{
if(kk != k) formk <- gsub(paste("cause",kk,sep=""),"",formk)
}
asurvk <- terms(formula(paste("~",formk)),specials=paste("cause",k,sep=""))
ind.causek <- attr(asurvk,"specials")$cause
if(length(ind.causek))
{
formcausek <- paste(labels(asurvk)[ind.causek],collapse="+")
formcausek <- gsub(paste("cause",k,sep=""),"",formcausek)
formcausek <- formula(paste("~",formcausek))
form.causek[[k]] <- formcausek
}
else
{
form.causek[[k]] <- ~-1
}
}
##creer la formule pour mixture
form2 <- form.surv
for( k in 1:nbevt)
{
form2 <- gsub(paste("cause",k,sep=""),"",form2)
}
form2 <- gsub("cause","",form2)
form2 <- formula(paste("~",form2))
asurv2 <- terms(form2,specials="mixture")
ind.mixture <- attr(asurv2,"specials")$mixture
if(length(ind.mixture))
{
form.mixture <- paste(labels(asurv2)[ind.mixture],collapse="+")
form.mixture <- gsub("mixture","",form.mixture)
form.mixture <- formula(paste("~",form.mixture))
}
else
{
form.mixture <- ~-1
}
## creer la formule pour ni cause ni mixture
asurv <- terms(formula(paste("~",form.surv)),specials=c("cause","mixture",paste("cause",1:nbevt,sep="")))
ind.commun <- setdiff(1:length(labels(asurv)),unlist(attr(asurv,"specials")))
if(length(ind.commun))
{
form.commun <- paste(labels(asurv)[ind.commun],collapse="+")
form.commun <- gsub("mixture","",form.commun) #si X1*mixture(X2), alors X1:mixture(X2) dans form.commun
form.commun <- gsub("cause","",form.commun) # si X1:cause(X2)
form.commun <- formula(paste("~",form.commun))
##NB: si mixture(X1)*cause(X2), X1:X2 en commun
}
else
{
form.commun <- ~-1
}
}
##verifier si toutes les variables sont dans data
variables <- c(attr(afixed,"variables"),attr(arandom,"variables"),attr(amixture,"variables"),attr(aclassmb,"variables"),attr(asurv,"variables"))
variables <- unlist(lapply(variables,all.vars))
if(!all(variables %in% colnames(data))) stop(paste("Data should contain the variables",paste(unique(variables),collapse=" ")))
###liste des variables utilisees (sans les interactions et sans Y)
if (ncor0>0) ttesLesVar <- unique(c(variables,cor.var.time))
else ttesLesVar <- unique(c(variables))
nomY <- as.character(attr(afixed,"variables")[2])
ttesLesVar <- setdiff(ttesLesVar,nomY)
###subset de data avec les variables utilisees
newdata <- data[,unique(c(nom.subject,nomY,noms.surv,ttesLesVar,nom.prior,nom.pprior,nom.timedepvar)),drop=FALSE]
## remplacer les NA de prior par 0
if(!is.null(nom.prior))
{
prior <- newdata[,nom.prior]
newdata[which(is.na(prior)),nom.prior] <- 0
prior[which(is.na(prior))] <- 0
}
##pprior : proba d'appartenance a chaque classe
if(!is.null(pprior))
{
pprior <- newdata[,nom.pprior]
}
## remplacer les NA de TimeDepVar par Tevent
Tint <- Tevent
nvdepsurv <- 0
if(!is.null(nom.timedepvar))
{
Tint <- newdata[,nom.timedepvar]
Tint[(is.na(Tint))] <- Tevent[(is.na(Tint))]
Tint[Tint>Tevent] <- Tevent[Tint>Tevent]
Tint[Tint<Tentry] <- Tentry[Tint<Tentry]
nvdepsurv <- 1
if (length(Tint[Tint<Tevent])==0)
{
stop("TimeDepVar is always greater than Time of Event. \n")
nvdepsurv <- 0
}
if (length(Tint[Tint>Tentry])==0)
{
Tint <- Tevent
stop("TimeDepVar is always lower than Time of Entry (0 by default). \n")
nvdepsurv <- 0
}
newdata[,nom.timedepvar] <- Tint
}
##enlever les NA
linesNA <- apply(newdata,2,function(v) which(is.na(v)))
linesNA <- unique(unlist(linesNA))
if(length(linesNA))
{
if(na.action==1) newdata <- newdata[-linesNA,,drop=FALSE]
if(na.action==2) stop("Data contain missing values.")
Tentry <- Tentry[-linesNA]
Tevent <- Tevent[-linesNA]
Event <- Event[-linesNA]
Tint <- Tint[-linesNA]
prior <- prior[-linesNA]
pprior <- pprior[-linesNA,]
}
### Y0
Y0 <- newdata[,nomY]
###creation de X0 (ttes les var + interactions)
Xintercept <- model.matrix(~1,data=newdata) #pr etre sur d'avoir I en premier
Xfixed <- model.matrix(fixed[-2], data=newdata)
Xmixture <- model.matrix(mixture, data=newdata)
Xrandom <- model.matrix(random, data=newdata)
Xclassmb <- model.matrix(classmb, data=newdata)
Xsurv <- model.matrix(form.commun,data=newdata)
Xsurvmix <- model.matrix(form.mixture,data=newdata)
Xsurvcause <- model.matrix(form.cause,data=newdata)
for (k in 1:nbevt)
{
assign(paste("Xsurvcause",k,sep=""),model.matrix(form.causek[[k]],data=newdata))
}
z.fixed <- strsplit(colnames(Xfixed),split=":",fixed=TRUE)
z.fixed <- lapply(z.fixed,sort)
if(random != ~-1)
{
z.random <- strsplit(colnames(Xrandom),split=":",fixed=TRUE)
z.random <- lapply(z.random,sort)
}
else
{
z.random <- list()
}
if(mixture != ~-1)
{
z.mixture <- strsplit(colnames(Xmixture),split=":",fixed=TRUE)
z.mixture <- lapply(z.mixture,sort)
}
else
{
z.mixture <- list()
}
if(classmb != ~-1)
{
z.classmb <- strsplit(colnames(Xclassmb),split=":",fixed=TRUE)
z.classmb <- lapply(z.classmb,sort)
}
else
{
z.classmb <- list()
}
if(form.commun != ~-1)
{
z.surv <- strsplit(colnames(Xsurv),split=":",fixed=TRUE)
z.surv <- lapply(z.surv,sort)
}
else
{
z.surv <- list()
}
if(form.mixture != ~-1)
{
z.survmix <- strsplit(colnames(Xsurvmix),split=":",fixed=TRUE)
z.survmix <- lapply(z.survmix,sort)
}
else
{
z.survmix <- list()
}
if(form.cause != ~-1)
{
z.survcause <- strsplit(colnames(Xsurvcause),split=":",fixed=TRUE)
z.survcause <- lapply(z.survcause,sort)
}
else
{
z.survcause <- list()
}
for(k in 1:nbevt)
{
if(form.causek[[k]] != ~-1)
{
assign(paste("z.survcause",k,sep=""),strsplit(colnames(get(paste("Xsurvcause",k,sep=""))),split=":",fixed=TRUE))
assign(paste("z.survcause",k,sep=""),lapply(get(paste("z.survcause",k,sep="")),sort))
}
else
{
assign(paste("z.survcause",k,sep=""),list())
}
}
if(!all(z.mixture %in% z.fixed)) stop("The covariates in mixture should also be included in the argument fixed")
X0 <- cbind(Xintercept,Xfixed, Xrandom, Xclassmb,Xsurv,Xsurvmix,Xsurvcause)
for(k in 1:nbevt)
{
X0 <- cbind(X0,get(paste("Xsurvcause",k,sep="")))
}
nom.unique <- unique(colnames(X0))
X0 <- X0[,nom.unique,drop=FALSE]
if(ncor0>0)
{
if(!(cor.var.time %in% colnames(X0)))
{
X0 <- cbind(X0, newdata[,cor.var.time])
colnames(X0) <- c(nom.unique, cor.var.time)
nom.unique <- colnames(X0)
}
}
X0 <- as.matrix(X0)
###X0 fini
timeobs <- rep(0, nrow(newdata))
if(!is.null(var.time))
{
timeobs <- newdata[,var.time]
if(any(is.na(timeobs))) stop(paste("Cannot use",var.time,"as time variable because it contains missing data"))
}
###test de link
if(link %in% c("splines","Splines"))
{
link <- "5-quant-splines"
}
idlink <- 2
if(link=="linear") idlink <- 0
if(link=="beta") idlink <- 1
if(link=="NULL") idlink <- -1
if (idlink==1)
{
if (epsY<=0)
{
stop("Argument 'epsY' should be positive.")
}
}
##remplir range si pas specifie
if(!is.null(range) & length(range)!=2) stop("Vector range should be of length 2.")
if(length(range)==2)
{
if(max(Y0)>range[2] | min(Y0)<range[1]) stop("The range specified do not cover the entire range of the data")
}
if(is.null(range))
{
min1 <- min(Y0)
min2 <- round(min1,3)
if(min1<min2) min2 <- min2-0.001
max1 <- max(Y0)
max2 <- round(max1,3)
if(max1>max2) max2 <- max2+0.001
range <- c(min2,max2)
}
nbzitr <- 2 #nbzitr = nb de noeuds si splines, 2 sinon
spltype <- NULL
if(idlink==2)
{
spl <- strsplit(link,"-")
if(length(spl[[1]])!=3) stop("Invalid argument link")
nbzitr <- as.numeric(spl[[1]][1])
spltype <- spl[[1]][2]
if(!(spltype %in% c("equi","quant","manual"))) stop("The location of the nodes should be 'equi', 'quant' or 'manual'")
if(!(spl[[1]][3] %in% c("splines","Splines"))) stop("Invalid argument link")
if(!is.null(intnodes))
{
if(spltype != "manual") stop("intnodes should be NULL if the nodes are not chosen manually")
if(length(intnodes) != (nbzitr-2)) stop(paste("Vector intnodes should be of length",nbzitr-2))
}
}
##remplir zitr (contient les noeuds si splines, sinon min et max)
zitr0 <- rep(0,nbzitr)
if(idlink %in% c(0,1)) zitr0[1:2] <- range
if(idlink==2)
{
if(spltype=="manual")
{
zitr0[1] <- range[1]
zitr0[nbzitr] <- range[2]
zitr0[2:(nbzitr-1)] <- intnodes
##verifier s'il y a des obs entre les noeuds
hcounts <- hist(Y0,breaks=zitr0,plot=FALSE,include.lowest=TRUE,right=TRUE)$counts
if(any(hcounts==0)) stop("Link function can not be estimated. Please try other nodes such that there are observations in each interval.")
}
if(spltype=="equi") zitr0[1:nbzitr] <- seq(range[1],range[2],length.out=nbzitr)
if(spltype=="quant")
{
Y0bis <- Y0
if(range[1]<min(Y0)) Y0bis <- c(range[1],Y0)
if(range[2]>max(Y0)) Y0bis <- c(range[2],Y0bis)
zitr0[1:nbzitr] <- quantile(Y0bis,probs=seq(0,1,length.out=nbzitr))
}
}
###uniqueY0 et indiceY0
uniqueY0 <- 0
indiceY0 <- 0
nvalSPL0 <- 0
if(idlink == 2)
{
uniqueY0 <- sort(unique(Y0))
permut <- order(order(Y0)) # sort(y)[order(order(y))] = y
indice <- rep(1:length(uniqueY0), as.vector(table(Y0)))
indiceY0 <- indice[permut]
nvalSPL0 <- length(uniqueY0)
}
###ordonner les mesures par individu
IND <- newdata[,nom.subject]
#IDnum <- as.numeric(IND)
if(is.null(nom.prior)) prior <- rep(0,length(Y0))
if(is.null(nom.pprior)) pprior <- matrix(1,length(Y0),ng)
if(!length(indiceY0)) indiceY0 <- rep(0,length(Y0))
matYX <- cbind(IND,timeobs,prior,pprior,Y0,indiceY0,Tentry,Tevent,Event,Tint,X0)
matYXord <- matYX[order(IND),]
Y0 <- as.numeric(matYXord[,4+ng])
X0 <- apply(matYXord[,-c(1:(9+ng)),drop=FALSE],2,as.numeric)
#IDnum <- matYXord[,1]
IND <- matYXord[,1]
indiceY0 <- as.numeric(matYXord[,5+ng])
prior0 <- as.numeric(unique(matYXord[,c(1,3)])[,2])
if(length(prior0)!=length(unique(IND))) stop("Please check 'prior' argument. Subjects can not have multiple assigned classes.")
pprior0 <- t(matYXord[!duplicated(matYXord[,1]), 3+1:ng, drop=FALSE])
timeobs <- matYXord[,2]
## nombre de sujets
ns0 <- length(unique(IND))
### Tevent, Tentry et Event de dim ns
#data.surv <- unique(cbind(IND,matYX[,c(6,7,8,9)]))
#if(nrow(data.surv) != ns0) stop("Subjects cannot have several times to event.")
nmes <- as.vector(table(IND))
##data.surv <- sweep(matYXord[cumsum(nmes),ng+c(6,7,8,9), drop=FALSE], MARGIN=1, STATS=0, FUN=function(x,y){as.numeric(x)+y})
data.surv <- matYXord[cumsum(nmes),ng+c(6,7,8,9), drop=FALSE]
tsurv0 <- as.numeric(data.surv[,1])
tsurv <- as.numeric(data.surv[,2])
devt <- as.numeric(data.surv[,3])
tsurvint <- as.numeric(data.surv[,4])
ind_survint <- (tsurvint<tsurv) + 0
### test de hazard
arghaz <- hazard
hazard <- rep(hazard,length.out=nbevt)
if(any(hazard %in% c("splines","Splines")))
{
hazard[which(hazard %in% c("splines","Splines"))] <- "5-quant-splines"
}
if(any(hazard %in% c("piecewise","Piecewise")))
{
hazard[which(hazard %in% c("piecewise","Piecewise"))] <- "5-quant-piecewise"
}
haz13 <- strsplit(hazard[which(!(hazard=="Weibull"))],"-")
if(any(sapply(haz13,length)!=3)) stop("Invalid argument hazard")
## si plusieurs splines, noeuds doivent etre identiques
nbspl <- length(which(sapply(haz13,getElement,3) %in% c("splines","Splines")))
if(nbspl>1)
{
haz3 <- haz13[which(sapply(haz13,getElement,3) %in% c("splines","Splines"))]
if(length(unique(sapply(haz3,getElement,1)))>1) stop("The nodes location of all splines hazard functions must be identical")
if(length(unique(sapply(haz3,getElement,2)))>1) stop("The nodes location of all splines hazard functions must be identical")
}
nz <- rep(2,nbevt)
locnodes <- NULL
typrisq <- rep(2,nbevt)
nprisq <- rep(2,nbevt)
nbnodes <- 0 #longueur de hazardnodes
nbrange <- 0 #longueur de hazardrange
ii <- 0
dejaspl <- 0
if(any(hazard!="Weibull"))
{
for (i in 1:nbevt)
{
if(hazard[i]=="Weibull") next;
ii <- ii+1
nz[i] <- as.numeric(haz13[[ii]][1])
if(nz[i]<3) stop("At least 3 nodes are required")
typrisq[i] <- ifelse(haz13[[ii]][3] %in% c("splines","Splines"),3,1)
nprisq[i] <- ifelse(haz13[[ii]][3] %in% c("splines","Splines"),nz[i]+2,nz[i]-1)
locnodes <- c(locnodes, haz13[[ii]][2])
if(!(haz13[[ii]][3] %in% c("splines","Splines","piecewise","Piecewise"))) stop("Invalid argument hazard")
if((haz13[[ii]][2]=="manual"))
{
if(typrisq[i]==1 | dejaspl==0)
{
if(length(arghaz)>1 | i==1 )
{
nbnodes <- nbnodes + nz[i]-2
nbrange <- nbrange + 2
}
}
if(typrisq[i]==3) dejaspl <- 1
}
if(!all(locnodes %in% c("equi","quant","manual"))) stop("The location of the nodes should be 'equi', 'quant' or 'manual'")
}
if(!is.null(hazardnodes))
{
if(!any(locnodes == "manual")) stop("hazardnodes should be NULL if the nodes are not chosen manually")
if(length(hazardnodes) != nbnodes) stop(paste("Vector hazardnodes should be of length",nbnodes))
}
}
else
{
if(!is.null(hazardnodes)) stop("hazardnodes should be NULL if Weibull baseline risk functions are chosen")
}
if(nbevt>1 & length(arghaz)==1 & nbnodes>0)
{
hazardnodes <- rep(hazardnodes,length.out=nbnodes*nbevt)
}
zi <- matrix(0,nrow=max(nz),ncol=nbevt)
nb <- 0
minT1 <- 0
maxT1 <- max(tsurv)
tsurvevt <- tsurv #[which(devt!=0)]
if(idtrunc==1)
{
minT1 <- min(tsurv,tsurv0)
maxT1 <- max(tsurv,tsurv0)
}
## arrondir
minT2 <- round(minT1,3)
if(minT1<minT2) minT2 <- minT2-0.001
minT <- minT2
maxT2 <- round(maxT1,3)
if(maxT1>maxT2) maxT2 <- maxT2+0.001
maxT <- maxT2
if(!is.null(hazardrange))
{
if(length(hazardrange) != nbrange) stop(paste("Vector hazardrange should be of length", nbrange))
if(any(hazardrange[seq(2,length(hazardrange), 2)] < maxT1) | any(hazardrange[seq(1,length(hazardrange), 2)] > minT1)) stop("The hazard range specified do not cover the entire range of the data")
if(nbevt>1 & length(arghaz)==1)
{
hazardrange <- rep(hazardrange,length.out=nbrange*nbevt)
}
}
ii <- 0
for(i in 1:nbevt)
{
if(typrisq[i]==2)
{
zi[1:2,i] <- c(minT,maxT)
}
else
{
ii <- ii+1
minTii <- minT
maxTii <- maxT
if(!is.null(hazardrange))
{
minTii <- hazardrange[2*(ii-1) + 1]
maxTii <- hazardrange[2*ii]
}
if(locnodes[ii]=="manual")
{
zi[1:nz[i],i] <- c(minTii,hazardnodes[nb+1:(nz[i]-2)],maxTii)
nb <- nb + nz[i]-2
}
if(locnodes[ii]=="equi")
{
zi[1:nz[i],i] <- seq(minTii,maxTii,length.out=nz[i])
}
if(locnodes[ii]=="quant")
{
#pi <- seq(0,1,length.out=nz[i])
#pi <- pi[-length(pi)]
#pi <- pi[-1]
pi <- c(1:(nz[i]-2))/(nz[i]-1)
qi <- quantile(tsurvevt,prob=pi)
zi[1,i] <- minTii
zi[2:(nz[i]-1),i] <- qi
zi[nz[i],i] <- maxTii
}
}
}
hazardtype <- rep(hazardtype,length.out=nbevt)
risqcom <- (hazardtype=="Common") + (hazardtype=="PH")*2
nrisq <- (risqcom==1)*nprisq + (risqcom==0)*nprisq*ng + (risqcom==2)*(nprisq+ng-1)
###parametres pour fortran
ng0 <- ng
nv0 <- dim(X0)[2]
nobs0 <- length(Y0)
idiag0 <- ifelse(idiag==TRUE,1,0)
nwg0 <- ifelse(nwg==TRUE,1,0)
nrisqtot <- sum(nrisq)
loglik <- 0
ni <- 0
istop <- 0
gconv <- rep(0,3)
ppi0 <- rep(0,ns0*ng0)
ppitest0 <- rep(0,ns0*ng0)
resid_m <- rep(0,nobs0)
resid_ss <- rep(0,nobs0)
pred_m_g <- rep(0,nobs0*ng0)
pred_ss_g <- rep(0,nobs0*ng0)
Yobs <- rep(0,nobs0)
#rlindiv <- rep(0,ns0)
marker <- rep(0,nsim)
transfY <- rep(0,nsim)
#Ydiscret <- 0
#UACV <- 0
#vraisdiscret <- 0
nmes0 <- as.vector(table(IND))
logspecif <- as.numeric(logscale)
time <- seq(minT,maxT,length.out=nsim)
risq_est <- matrix(0,nrow=nsim*ng0,ncol=nbevt)
risqcum_est <- matrix(0,nrow=nsim*ng0,ncol=nbevt)
statglob <- 0
statevt <- rep(0,nbevt)
###remplir idprob, etc
z.X0 <- strsplit(nom.unique,split=":",fixed=TRUE)
z.X0 <- lapply(z.X0,sort)
idprob0 <- z.X0 %in% z.classmb + 0
idea0 <- z.X0 %in% z.random + 0
idg0 <- (z.X0 %in% z.fixed) + (z.X0 %in% z.mixture)
## indicateurs pr variables survie
idcause <- rep(0,nv0)
for(k in 1:nbevt)
{
idcause <- idcause + (z.X0 %in% get(paste("z.survcause",k,sep="")) )
}
idcom <- rep(0,nv0)
idspecif <- matrix(0,ncol=nv0,nrow=nbevt)
for(j in 1:nv0)
{
if((z.X0[j] %in% z.surv) & !(z.X0[j] %in% z.survcause) & (idcause[j]==0))
{
idcom[j] <- 1
idspecif[,j] <- 1
#print("X")
}
if((z.X0[j] %in% z.survmix) & ( !(z.X0[j] %in% z.survcause) & (idcause[j]==0)))
{
idcom[j] <- 1
idspecif[,j] <- 2
#print("mixture(X)")
}
if((z.X0[j] %in% z.survmix) & (z.X0[j] %in% z.survcause))
{
idcom[j] <- 0
idspecif[,j] <- 2
#print("cause(mixture(X))")
}
if((z.X0[j] %in% z.survcause) & (!(z.X0[j] %in% z.survmix)))
{
idcom[j] <- 0
idspecif[,j] <- 1
#print("cause(X)")
}
if(idcause[j]!=0)
{
if(z.X0[j] %in% z.survmix)
{
for(k in 1:nbevt)
{
if(z.X0[j] %in% get(paste("z.survcause",k,sep="")))
{
idcom[j] <- 0
idspecif[k,j] <- 2
#cat("causek(mixture(X)) ,k=",k,"\n")
}
}
}
else
{
for(k in 1:nbevt)
{
if(z.X0[j] %in% get(paste("z.survcause",k,sep="")))
{
idcom[j] <- 0
idspecif[k,j] <- 1
#cat("causek(X) ,k=",k,"\n")
}
}
}
}
}
## mettre des 0 pour l'intercept
idcom[1] <- 0
idspecif[,1] <- 0
## quelle variable est TimeDepVar
idtdv <- z.X0 %in% nom.timedepvar + 0
#if(any(idcom>1) & nbevt<2) stop("No event specific effect can be estimated with less than two events")
if (ncor0>0) idcor0 <- z.X0 %in% cor.var.time +0
else idcor0 <- rep(0,nv0)
## Si pas TimeDepVar dans formule survival
if(length(nom.timedepvar) & all(idtdv==0))
{
stop("Variable in 'TimeDepVar' should also appear as a covariate in the 'survival' argument")
## ou l'ajouter?
}
## si on a ignore TimeDepVar #pas utile si stop
if(length(nom.timedepvar) & nvdepsurv==0)
{
jj <- which(idtdv==1)
idtdv[jj] <- 0
idcom[jj] <- 0
idspecif[,jj] <- 0
idcause[jj] <- 0
nom.timedepvar <- NULL
}
nea0 <- sum(idea0)
predRE <- rep(0,ns0*nea0)
## nb coef pr survie
nvarxevt <- 0
nvarxevt2 <- 0 # pr valeurs initiales calculees avec Fortran
for(j in 1:nv0)
{
if(idcom[j]==1)
{
if(all(idspecif[,j]==1))
{
nvarxevt <- nvarxevt + 1
nvarxevt2 <- nvarxevt2 + 1
}
if(all(idspecif[,j]==2))
{
nvarxevt <- nvarxevt + ng
nvarxevt2 <- nvarxevt2 + 1
}
}
if(idcom[j]==0)
{
if(all(idspecif[,j]==0)) next
for(k in 1:nbevt)
{
if(idspecif[k,j]==1)
{
nvarxevt <- nvarxevt + 1
nvarxevt2 <- nvarxevt2 + 1
}
if(idspecif[k,j]==2)
{
nvarxevt <- nvarxevt + ng
nvarxevt2 <- nvarxevt2 + 1
}
}
}
}
ntrtot0 <- (idlink==-1) + 2*(idlink==0) + 4*(idlink==1) + (nbzitr+2)*(idlink==2)
nef <- sum(idg0==1) + ng0*sum(idg0==2)
if(idlink!=-1) nef <- nef-1
nvc <- ifelse(idiag0==1,nea0,nea0*(nea0+1)/2)
nw <- (ng0-1)*nwg0
nprob <- sum(idprob0)*(ng0-1)
##nombre total de parametres
NPM <- (ng0-1)*sum(idprob0) +
nrisqtot + nvarxevt +
nef + nvc + nw + ncor0 + ntrtot0
V <- rep(0, NPM*(NPM+1)/2) #pr variance des parametres
## prm fixes
fix0 <- rep(0,NPM)
if(length(posfix))
{
if(any(!(posfix %in% 1:NPM))) stop("Indexes in posfix are not correct")
fix0[posfix] <- 1
}
if(length(posfix)==NPM) stop("No parameters to estimate")
## pour H restreint
Hr0 <- as.numeric(partialH)
pbH0 <- rep(0,NPM)
if(is.logical(partialH))
{
if(partialH)
{
if(any(typrisq %in% c(1,3)))
{
for(k in 1:nbevt)
{
if(typrisq[k] %in% c(1,3))
{
pbH0[nprob+sum(nrisq[1:k])-nrisq[k]+1:nrisq[k]] <- 1
if(risqcom[k]==2)
{
pbH0[nprob+sum(nrisq[1:k])] <- 0
}
}
}
}
if(any(idlink==2))
{
pbH0[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+1:ntrtot0] <- 1
}
}
#pbH0[posfix] <- 0
if(sum(pbH0)==0 & Hr0==1) stop("No partial Hessian matrix can be defined in the absence of baseline risk function approximated by splines or piecewise linear function")
}
else
{
if(!all(Hr0 %in% 1:NPM)) stop("Indexes in partialH are not correct")
pbH0[Hr0] <- 1
#pbH0[posfix] <- 0
}
indexHr <- NULL
if(sum(pbH0)>0)
{
if(length(posfix)) pbH1 <- pbH0[-posfix] else pbH1 <- pbH0
indexHr <- which(pbH1==1)
}
## gestion de B=random(mod)
Brandom <- FALSE
tryB <- try(as.numeric(B), silent=TRUE)
if(inherits(tryB, "try-error"))
{
if(length(cl$B)==2)
{
if(!inherits(eval(cl$B[[2]], parent.env(environment())),"Jointlcmm")) stop("The model specified in B should be of class Jointlcmm")
if(as.character(cl$B[1])!="random") stop("Please use random() to specify random initial values")
Brandom <- TRUE
B <- eval(cl$B[[2]], parent.env(environment()))
if(B$conv != 1) stop("Model in argument B did not converge properly")
#if(length(posfix)) stop("Argument posfix is not compatible with random intial values")
}
}
##valeurs initiales
if(!(missing(B)))
{
if(is.vector(B))
{
if (length(B)==NPM)
{
b <- B
if(nvc>0)
{
## remplacer varcov des EA par les prm a estimer
if(idiag0==1)
{
b[nprob+nrisqtot+nvarxevt+nef+1:nvc] <- sqrt(b[nprob+nrisqtot+nvarxevt+nef+1:nvc])
}
else
{
varcov <- matrix(0,nrow=nea0,ncol=nea0)
varcov[upper.tri(varcov,diag=TRUE)] <- b[nprob+nrisqtot+nvarxevt+nef+1:nvc]
varcov <- t(varcov)
varcov[upper.tri(varcov,diag=TRUE)] <- b[nprob+nrisqtot+nvarxevt+nef+1:nvc]
ch <- chol(varcov)
b[nprob+nrisqtot+nvarxevt+nef+1:nvc] <- ch[upper.tri(ch,diag=TRUE)]
}
}
}
else
{
stop(paste("Vector B should be of length",NPM))
}
}
else
{
if(!inherits(B,"Jointlcmm")) stop("B should be either a vector or an object of class Jointlcmm")
b <- rep(0,NPM)
if(ng>1 & B$ng==1)
{
nef2 <- sum(idg0!=0)
if(idlink != -1)
{
nef2 <- nef2-1 #car intercept fixe a 0
}
NPM2 <- sum(nprisq)+nvarxevt2+nef2+nvc+ncor0+ntrtot0
if(length(B$best)!=NPM2) stop("B is not correct")
if(Brandom==FALSE)
{
## B deterministe
for(ke in 1:nbevt)
{
if(risqcom[ke]==0)
{
indb <- nprob+sum(nrisq[1:ke])-nrisq[ke]+1:nrisq[ke]
indinit <- sum(nprisq[1:ke])-nprisq[ke]+1:nprisq[ke]
if(B$conv==1)
{
b[indb] <- rep(B$best[indinit],ng0)+rep((1:ng0)-(ng0+1)/2,each=nprisq[ke])*rep(sqrt(B$V[indinit*(indinit+1)/2]),ng0)
}
else
{
b[indb] <- rep(B$best[indinit],ng0)+rep((1:ng0)-(ng0+1)/2,each=nprisq[ke])*rep(B$best[indinit],ng0)
}
}
if(risqcom[ke]==1)
{
b[nprob+sum(nrisq[1:ke])-nrisq[ke]+1:nrisq[ke]] <- B$best[sum(nprisq[1:ke])-nprisq[ke]+1:nprisq[ke]]
}
if(risqcom[ke]==2)
{
b[nprob+sum(nrisq[1:ke])-nrisq[ke]+1:nrisq[ke]] <- c(B$best[sum(nprisq[1:ke])-nprisq[ke]+1:nprisq[ke]],0.5+(0:(ng0-2))*0.5)
}
}
## pr bevt
avtj <- nprob+nrisqtot
avtj2 <- sum(nprisq)
for(j in 1:nv0)
{
if(idcom[j]==0 & all(idspecif[,j]==0)) next
if(idcom[j]==1 & all(idspecif[,j]==1))
{
b[avtj+1] <- B$best[avtj2+1]
avtj <- avtj+1
avtj2 <- avtj2+1
}
if(idcom[j]==1 & all(idspecif[,j]==2))
{
if(B$conv==1) b[avtj+1:ng0] <- abs(rep(B$best[avtj2+1],ng0)+c(1:ng0-(ng0+1)/2)*rep(sqrt(B$V[(avtj2+1)*(avtj2+2)/2]),ng0))
else b[avtj+1:ng0] <- abs(rep(B$best[avtj2+1],ng0)+c(1:ng0-(ng0+1)/2)*rep(B$best[avtj2+1],ng0))
avtj <- avtj+ng0
avtj2 <- avtj2+1
}
if(idcom[j]==0 & any(idspecif[,j]!=0))
{
for(k in 1:nbevt)
{
if(idspecif[k,j]==0) next
if(idspecif[k,j]==1)
{
b[avtj+1] <- B$best[avtj2+1]
avtj <- avtj+1
avtj2 <- avtj2+1
}
if(idspecif[k,j]==2)
{
if(B$conv==1) b[avtj+1:ng0] <- abs(rep(B$best[avtj2+1],ng0)+c(1:ng0-(ng0+1)/2)*rep(sqrt(B$V[(avtj2+1)*(avtj2+2)/2]),ng0))
else b[avtj+1:ng0] <- abs(rep(B$best[avtj2+1],ng0)+c(1:ng0-(ng0+1)/2)*rep(B$best[avtj2+1],ng0))
avtj <- avtj+ng0
avtj2 <- avtj2+1
}
}
}
}
## pr nef
avtj <- nprob+nrisqtot+nvarxevt
avtj2 <- sum(nprisq)+nvarxevt2
for(j in 1:nv0)
{
if(idg0[j]==1)
{
if(j==1 & idlink!=-1) next
else
{
if(B$conv==1) b[avtj+1] <- B$best[avtj2+1]
avtj <- avtj+1
avtj2 <- avtj2+1
}
}
if(idg0[j]==2)
{
if(j==1)
{
if(idlink!=-1)
{
b[avtj+1:(ng0-1)] <- -0.5*(1:(ng0-1))
avtj <- avtj+ng0-1
}
else
{
b[avtj+1:ng0] <- B$best[avtj2+1]+(1:ng0-(ng0+1)/2)*B$best[avtj2+1]
avtj <- avtj+ng0
avtj2 <- avtj2+1
}
}
else
{
if(B$conv==1) b[avtj+1:ng0] <- B$best[avtj2+1]+(1:ng0-(ng0+1)/2)*sqrt(B$V[(avtj2+1)*(avtj2+1+1)/2])
else b[avtj+1:ng0] <- B$best[avtj2+1]+(1:ng0-(ng0+1)/2)*B$best[avtj2+1]
avtj <- avtj+ng0
avtj2 <- avtj2+1
}
}
}
## pr nvc
if(nvc>0)
{
b[nprob+nrisqtot+nvarxevt+nef+1:nvc] <- B$cholesky
}
## cor et transfo
b[(nprob+nrisqtot+nvarxevt+nef+nvc+nw+1):NPM] <- B$best[(sum(nprisq)+nvarxevt2+nef2+nvc+1):NPM2]
}
else
{
## B random
bb <- rep(0,NPM-nprob-(ng-1)*sum(length(which(risqcom==2)))-nw) # tous les prm sauf nprob, bPH et nw
vbb <- matrix(0,length(bb),length(bb))
VB <- matrix(0,NPM2,NPM2)
VB[upper.tri(VB,diag=TRUE)] <- B$V
VB <- t(VB)
VB[upper.tri(VB,diag=TRUE)] <- B$V
copcov <- rep(0,length(bb))
copcov2 <- rep(0,length(B$best))
avt <- 0
for(ke in 1:nbevt)
{
if(risqcom[ke]==0)
{
indbb <- avt+1:nrisq[ke]
indB <- sum(nprisq[1:ke])-nprisq[ke]+1:nprisq[ke]
bb[indbb] <- rep(B$best[indB],ng)
diag(vbb[indbb,indbb]) <- rep(diag(VB[indB,indB]),ng)
avt <- avt + nrisq[ke]
}
if(risqcom[ke]==1)
{
indbb <- avt+1:nrisq[ke]
indB <- sum(nprisq[1:ke])-nprisq[ke]+1:nprisq[ke]
bb[indbb] <- B$best[indB]
copcov[indbb] <- 1
copcov2[indB] <- 1
avt <- avt + nrisq[ke]
}
if(risqcom[ke]==2)
{
indbb <- avt+1:nprisq[ke]
indB <- sum(nprisq[1:ke])-nprisq[ke]+1:nprisq[ke]
bb[indbb] <- B$best[indB]
copcov[indbb] <- 1
copcov2[indB] <- 1
avt <- avt + nprisq[ke]
}
}
avt2 <- sum(nprisq)
for(j in 1:nv0)
{
if(idcom[j]==0 & all(idspecif[,j]==0)) next
if(idcom[j]==1 & all(idspecif[,j]==1))
{
bb[avt+1] <- B$best[avt2+1]
copcov[avt+1] <- 1
copcov2[avt2+1] <- 1
avt <- avt+1
avt2 <- avt2+1
}
if(idcom[j]==1 & all(idspecif[,j]==2))
{
bb[avt+1:ng] <- rep(B$best[avt2+1],ng)
diag(vbb[avt+1:ng,avt+1:ng]) <- rep(VB[avt2+1,avt2+1],ng)
avt <- avt+ng
avt2 <- avt2+1
}
if(idcom[j]==0 & any(idspecif[,j]!=0))
{
for(k in 1:nbevt)
{
if(idspecif[k,j]==0) next
if(idspecif[k,j]==1)
{
bb[avt+1] <- B$best[avt2+1]
copcov[avt+1] <- 1
copcov2[avt2+1] <- 1
avt <- avt+1
avt2 <- avt2+1
}
if(idspecif[k,j]==2)
{
bb[avt+1:ng] <- rep(B$best[avt2+1],ng)
diag(vbb[avt+1:ng,avt+1:ng]) <- rep(VB[avt2+1,avt2+1],ng)
avt <- avt+ng
avt2 <- avt2+1
}
}
}
}
for(j in 1:nv0)
{
if(idg0[j]==1)
{
if(j==1 & idlink!=-1) next
else
{
bb[avt+1] <- B$best[avt2+1]
copcov[avt+1] <- 1
copcov2[avt2+1] <- 1
avt <- avt+1
avt2 <- avt2+1
}
}
if(idg0[j]==2)
{
if(j==1 & idlink!=-1)
{
avt <- avt+ng-1
}
else
{
bb[avt+1:ng] <- rep(B$best[avt2+1],ng)
diag(vbb[avt+1:ng,avt+1:ng]) <- rep(VB[avt2+1,avt2+1],ng)
avt <- avt+ng
avt2 <- avt2+1
}
}
}
if(nvc>0)
{
bb[avt+1:nvc] <- B$cholesky
copcov[avt+1:nvc] <- 1
copcov2[avt2+1:nvc] <- 1
avt <- avt+nvc
avt2 <- avt2+nvc
}
bb[(avt+1):length(bb)] <- B$best[(avt2+1):NPM2]
copcov[(avt+1):length(bb)] <- 1
copcov2[(avt2+1):length(B$best)] <- 1
vbb[which(copcov==1),which(copcov==1)] <- VB[which(copcov2==1),which(copcov2==1)]
if(idlink!=-1 & idg0[1]>1)
{
bb <- bb[-(sum(nrisq)-(ng-1)*sum(length(which(risqcom==2)))+nvarxevt+1:(ng-1))]
vbb <- vbb[-(sum(nrisq)-(ng-1)*sum(length(which(risqcom==2)))+nvarxevt+1:(ng-1)),-(sum(nrisq)-(ng-1)*sum(length(which(risqcom==2)))+nvarxevt+1:(ng-1))]
}
##Chol <- chol(vbb)
##Chol <- t(Chol)
bb <- rmvnorm(n=1,mean=bb,sigma = vbb) #bb + Chol %*% rnorm(length(bb))
if(nprob>0) b[1:nprob] <- 0
avt <- 0
for(ke in 1:nbevt)
{
if(risqcom[ke]==0)
{
b[nprob+sum(nrisq[1:ke])-nrisq[ke]+1:nrisq[ke]] <- bb[avt+1:nrisq[ke]]
avt <- avt + nrisq[ke]
}
if(risqcom[ke]==1)
{
b[nprob+sum(nrisq[1:ke])-nrisq[ke]+1:nrisq[ke]] <- bb[avt+1:nrisq[ke]]
avt <- avt + nrisq[ke]
}
if(risqcom[ke]==2)
{
b[nprob+sum(nrisq[1:ke])-nrisq[ke]+1:nprisq[ke]] <- bb[avt+1:nprisq[ke]]
b[nprob+sum(nrisq[1:ke])-nrisq[ke]+nprisq[ke]+1:(ng-1)] <- 1 #bPH
avt <- avt + nprisq[ke]
}
}
if(nvarxevt>0) b[nprob+nrisqtot+1:nvarxevt] <- bb[avt+1:nvarxevt]
if(idlink!=-1 & idg0[1]>1)
{
b[nprob+nrisqtot+nvarxevt+1:(ng-1)] <- 0
b[(nprob+nrisqtot+nvarxevt+ng):(nprob+nrisqtot+nvarxevt+nef)] <- bb[avt+nvarxevt+1:(nef-(ng-1))]
nefssI <- nef-ng+1
}
else
{
b[nprob+nrisqtot+nvarxevt+1:nef] <- bb[avt+nvarxevt+1:nef]
nefssI <- nef
}
if(nvc>0)
{
b[nprob+nrisqtot+nvarxevt+nef+1:nvc] <- bb[nrisqtot-((ng-1)*length(which(risqcom==2)))+nvarxevt+nefssI+1:nvc]
}
if(nw>0) b[nprob+nrisqtot+nvarxevt+nef+nvc+1:nw] <- 1
b[(nprob+nrisqtot+nvarxevt+nef+nvc+nw+1):NPM] <- bb[(nrisqtot-((ng-1)*length(which(risqcom==2)))+nvarxevt+nefssI+nvc+1):length(bb)]
}
}
}
}
else
{
if(ng>1) stop("Please specify initial values with argument 'B'")
b <- rep(0,NPM)
for(i in 1:nbevt)
{
if (typrisq[i]==2)
{
if(logspecif==1)
{
b[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- c(rep(c(log(sum(devt==i)/sum(tsurv[devt==i])),0),ifelse(risqcom[i]==0,ng,1)),rep(1,(ng-1)*(risqcom[i]==2)))
}
else
{
b[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- c(rep(c(sqrt(sum(devt==i)/sum(tsurv[devt==i])),1),ifelse(risqcom[i]==0,ng,1)),rep(1,(ng-1)*(risqcom[i]==2)))
}
}
else
{
if(logspecif==1)
{
b[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- c(rep(log(1/nprisq[i]),ifelse(risqcom[i]==0,ng*nprisq[i],nprisq[i])),rep(1,(ng-1)*(risqcom[i]==2)))
}
else
{
b[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- c(rep(sqrt(1/nprisq[i]),ifelse(risqcom[i]==0,ng*nprisq[i],nprisq[i])),rep(1,(ng-1)*(risqcom[i]==2)))
}
}
}
if (nvc>0)
{
if(idiag0==1) b[nprob+nrisqtot+nvarxevt+nef+1:nvc] <- rep(1,nvc)
if(idiag0==0)
{
init.nvc <- diag(nea0)
init.nvc <- init.nvc[upper.tri(init.nvc, diag=TRUE)]
b[nprob+nrisqtot+nvarxevt+nef+1:nvc] <- init.nvc
}
}
if(nwg0>0) b[nprob+nrisqtot+nvarxevt+nef+nvc+1:nw] <- 1
if(ncor0==1) b[nprob+nrisqtot+nvarxevt+nef+nvc+nw+1] <- 1
if(ncor0==2) b[nprob+nrisqtot+nvarxevt+nef+nvc+nw+1:2] <- c(0,1)
if(idlink==-1) b[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+1] <- 1
if(idlink==0)
{
b[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+1] <- mean(Y0)
b[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+2] <- 1
}
if(idlink==1)
{
b[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+1] <- 0
b[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+2] <- -log(2)
b[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+3] <- 0.7
b[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+4] <- 0.1
}
if(idlink==2)
{
b[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+1] <- -2
b[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+2:ntrtot0] <- 0.1
}
#cat("B inti pour ng=1 ", b,"\n")
}
#cat("B init a partir du modele pour ng=1 ", b,"\n")
### nom au vecteur best
nom.X0 <- colnames(X0)
nom.X0[which(nom.X0=="(Intercept)")] <- "intercept"
##prm classmb
if(ng0>=2 & nprob>0)
{
nom <- rep(nom.X0[idprob0==1],each=ng0-1)
nom1 <- paste(nom," class",c(1:(ng0-1)),sep="")
names(b)[1:nprob] <- nom1
}
##prm fct de risque
if(isTRUE(logscale))
{
for(i in 1:nbevt)
{
nom1 <- rep(paste("event",i,sep=""),nrisq[i])
if(typrisq[i]==2)
{
nom2 <- paste(nom1[1:2]," log(Weibull",1:2,")",sep="")
nom1[1:(2*ifelse(risqcom[i]==0,ng,1))] <- rep(nom2,ifelse(risqcom[i]==0,ng*(risqcom[i]==0),1))
if(risqcom[i]==2) nom1[2+1:(ng-1)] <- paste(nom1[2+1:(ng-1)],"SurvPH")
names(b)[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- nom1
}
if(typrisq[i]==1)
{
nom2 <- paste(nom1[1:(nz[i]-1)]," log(piecewise",1:(nz[i]-1),")",sep="")
nom1[1:((nz[i]-1)*ifelse(risqcom[i]==0,ng,1))] <- rep(nom2,ifelse(risqcom[i]==0,ng*(risqcom[i]==0),1))
if(risqcom[i]==2) nom1[nz[i]-1+1:(ng-1)] <- paste(nom1[nz[i]-1+1:(ng-1)],"SurvPH")
names(b)[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- nom1
}
if(typrisq[i]==3)
{
nom2 <- paste(nom1[1:(nz[i]-1)]," log(splines",1:(nz[i]+2),")",sep="")
nom1[1:((nz[i]+2)*ifelse(risqcom[i]==0,ng,1))] <- rep(nom2,ifelse(risqcom[i]==0,ng*(risqcom[i]==0),1))
if(risqcom[i]==2) nom1[nz[i]+2+1:(ng-1)] <- paste(nom1[nz[i]+2+1:(ng-1)],"SurvPH")
names(b)[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- nom1
}
}
}
else
{
for(i in 1:nbevt)
{
nom1 <- rep(paste("event",i,sep=""),nrisq[i])
if(typrisq[i]==2)
{
nom2 <- paste(nom1[1:2]," +/-sqrt(Weibull",1:2,")",sep="")
nom1[1:(2*ifelse(risqcom[i]==0,ng,1))] <- rep(nom2,ifelse(risqcom[i]==0,ng*(risqcom[i]==0),1))
if(risqcom[i]==2) nom1[2+1:(ng-1)] <- paste(nom1[2+1:(ng-1)],"SurvPH")
names(b)[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- nom1
}
if(typrisq[i]==1)
{
nom2 <- paste(nom1[1:(nz[i]-1)]," +/-sqrt(piecewise",1:(nz[i]-1),")",sep="")
nom1[1:((nz[i]-1)*ifelse(risqcom[i]==0,ng,1))] <- rep(nom2,ifelse(risqcom[i]==0,ng*(risqcom[i]==0),1))
if(risqcom[i]==2) nom1[nz[i]-1+1:(ng-1)] <- paste(nom1[nz[i]-1+1:(ng-1)],"SurvPH")
names(b)[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- nom1
}
if(typrisq[i]==3)
{
nom2 <- paste(nom1[1:(nz[i]-1)]," +/-sqrt(splines",1:(nz[i]+2),")",sep="")
nom1[1:((nz[i]+2)*ifelse(risqcom[i]==0,ng,1))] <- rep(nom2,ifelse(risqcom[i]==0,ng*(risqcom[i]==0),1))
if(risqcom[i]==2) nom1[nz[i]+2+1:(ng-1)] <- paste(nom1[nz[i]+2+1:(ng-1)],"SurvPH")
names(b)[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- nom1
}
}
}
if(ng>1)
{
for(i in 1:nbevt)
{
if(risqcom[i]==1) next;
if(risqcom[i]==0)
{
names(b)[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- paste(names(b)[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]],paste("class",rep(1:ng,each=nprisq[i])))
}
if(risqcom[i]==2)
{
names(b)[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- paste(names(b)[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]],c(rep("",nprisq[i]),paste(" class",1:(ng-1),sep="")),sep="")
}
}
}
##prm covariables survival
nom1 <- NULL
for(j in 1:nv0)
{
if(idcom[j]==0 & all(idspecif[,j]==0)) next
if(idcom[j]==1 & all(idspecif[,j]==1)) #X
{
if(idtdv[j]==1)
{
nom1 <- c(nom1,paste("I(t>",nom.timedepvar,")",sep=""))
}
else
{
nom1 <- c(nom1,nom.X0[j])
}
next
}
if(idcom[j]==1 & all(idspecif[,j]==2)) #mixture(X)
{
if(idtdv[j]==1)
{
nom1 <- c(nom1,paste("I(t>",nom.timedepvar,") class",1:ng0,sep=""))
}
else
{
nom1 <- c(nom1,paste(nom.X0[j],paste("class",1:ng0,sep="")))
}
next
}
if(idcom[j]==0 & all(idspecif[,j]==1)) #cause(X)
{
if(idtdv[j]==1)
{
nom1 <- c(nom1,paste("I(t>",nom.timedepvar,") event",1:nbevt,sep=""))
}
else
{
nom1 <- c(nom1,paste(nom.X0[j],paste("event",1:nbevt,sep="")))
}
next
}
if(idcom[j]==0 & all(idspecif[,j]==2)) #cause(mixture(X))
{
if(idtdv[j]==1)
{
xevt <- paste("I(t>",nom.timedepvar,") event",1:nbevt,sep="")
classg <- paste("class",1:ng0,sep="")
nom1 <- c(nom1,paste(rep(xevt,each=ng0),rep(classg,nbevt)))
}
else
{
xevt <- paste(nom.X0[j],paste("event",1:nbevt,sep=""))
classg <- paste("class",1:ng0,sep="")
nom1 <- c(nom1,paste(rep(xevt,each=ng0),rep(classg,nbevt)))
}
next
}
if(idcom[j]==0 & idcause[j]!=0) #causek
{
for(k in 1:nbevt)
{
if(idspecif[k,j]==0) next
if(idspecif[k,j]==1)
{
if(idtdv[j]==1)
{
nom1 <- c(nom1,paste("I(t>",nom.timedepvar,") event",k,sep=""))
}
else
{
nom1 <- c(nom1,paste(nom.X0[j],paste("event",k,sep="")))
}
next
}
if(idspecif[k,j]==2)
{
if(idtdv[j]==1)
{
xevtk <- paste("I(t>",nom.timedepvar,") event",k,sep="")
classg <- paste("class",1:ng0,sep="")
nom1 <- c(nom1,paste(xevtk,classg))
}
else
{
xevtk <- paste(nom.X0[j],paste("event",k,sep=""))
classg <- paste("class",1:ng0,sep="")
nom1 <- c(nom1,paste(xevtk,classg))
}
next
}
}
}
}
if(nvarxevt>0) names(b)[nprob+nrisqtot+1:nvarxevt] <- nom1
##NB : pour chaque variable, coef ranges par evenement puis par classe
##prm fixed
if(ng0==1 & nef>0)
{
if(idlink!=-1)
{
names(b)[nprob+nrisqtot+nvarxevt+1:nef] <- nom.X0[-1][idg0[-1]!=0]
}
else
{
# print(nprob+nrisqtot+nvarxevt)
# print(nef)
# print(nom.X0)
# print(idg0)
names(b)[nprob+nrisqtot+nvarxevt+1:nef] <- nom.X0[idg0!=0]
}
}
if(ng0>1)
{
nom1<- NULL
for(i in 1:nv0)
{
if(idg0[i]==2)
{
if(i==1)
{
if(idlink==-1)
{
nom <- paste("intercept class",c(1:ng0),sep="")
nom1 <- c(nom1,nom)
}
else
{
nom <- paste("intercept class",c(2:ng0),sep="")
nom1 <- c(nom1,nom)
}
}
if (i>1)
{
nom <- paste(nom.X0[i]," class",c(1:ng0),sep="")
nom1 <- c(nom1,nom)
}
}
if(idg0[i]==1 & i==1 & idlink==-1) nom1 <- c(nom1,nom.X0[i])
if(idg0[i]==1 & i>1) nom1 <- c(nom1,nom.X0[i])
}
names(b)[nprob+nrisqtot+nvarxevt+1:nef]<- nom1
}
##prm random et nwg
if(nvc!=0) names(b)[nprob+nrisqtot+nvarxevt+nef+1:nvc] <- paste("varcov",c(1:nvc))
if(nw!=0) names(b)[nprob+nrisqtot+nvarxevt+nef+nvc+1:nw] <- paste("varprop class",c(1:(ng0-1)))
##prm cor
if(ncor0>0) {names(b)[nprob+nrisqtot+nvarxevt+nef+nvc+nw+1:ncor0] <- paste("cor",1:ncor0,sep="")}
##prm link
if(idlink==-1) names(b)[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+1] <- "stderr"
if(idlink==0) names(b)[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+1:2]<- c("Linear 1","Linear 2")
if(idlink==1) names(b)[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+1:4]<- paste("Beta",1:4,sep="")
if(idlink==2) names(b)[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+1:ntrtot0]<- paste("I-splines",c(1:ntrtot0),sep="")
names_best <- names(b)
#browser()
#print(b)
#### estimation
idspecif <- as.vector(t(idspecif))
#return(b)
nfix <- sum(fix0)
bfix <- 0
if(nfix>0)
{
bfix <- b[which(fix0==1)]
b <- b[which(fix0==0)]
NPM <- NPM-nfix
}
if(maxiter==0)
{
vrais <- loglikJointlcmm(b,Y0,X0,prior0,pprior0,tsurv0,tsurv,devt,ind_survint,
idprob0,idea0,idg0,idcor0,idcom,idspecif,idtdv,
idlink,epsY,nbzitr,zitr0,uniqueY0,nvalSPL0,
indiceY0,typrisq,risqcom,nz,zi,ns0,ng0,nv0,
nobs0,nmes0,nbevt,nea0,nw,ncor0,idiag0,idtrunc,
logspecif,NPM,fix0,nfix,bfix)
out <- list(conv=2, V=rep(NA, length(b)), best=b,
ppi=rep(NA,ns0*ng0), ppitest=rep(NA,ns0*ng0),
predRE=rep(NA,ns0*nea0), Yobs=rep(NA,nobs0),
resid_m=rep(NA,nobs0), resid_ss=rep(NA,nobs0),
marker=rep(NA,nsim), transfY=rep(NA,nsim),
pred_m_g=rep(NA,nobs0*ng0), pred_ss_g=rep(NA,nobs0*ng0),
time=rep(NA,nsim), risq_est=rep(NA,nsim*ng0*nbevt),
risqcum_est=rep(NA,nsim*ng0*nbevt),
statglob=NA, statevt=rep(NA,nbevt),
gconv=rep(NA,3), niter=0, loglik=vrais)
}
else
{
res <- mla(b=b, m=length(b), fn=loglikJointlcmm,
clustertype=clustertype,.packages=NULL,
epsa=convB,epsb=convL,epsd=convG,partialH=indexHr,
digits=8,print.info=verbose,blinding=FALSE,
multipleTry=25,file="",
nproc=nproc,maxiter=maxiter,minimize=FALSE,
Y0=Y0,X0=X0,prior0=prior0,pprior0=pprior0,tentr0=tsurv0,tevt0=tsurv,devt0=devt,
ind_survint0=ind_survint,idprob0=idprob0,idea0=idea0,idg0=idg0,
idcor0=idcor0,idcom0=idcom,idspecif0=idspecif,idtdv0=idtdv,
idlink0=idlink,epsY0=epsY,nbzitr0=nbzitr,zitr0=zitr0,
uniqueY0=uniqueY0,nvalSPL0=nvalSPL0,indiceY0=indiceY0,
typrisq0=typrisq,risqcom0=risqcom,nz0=nz,zi0=zi,ns0=ns0,
ng0=ng0,nv0=nv0,nobs0=nobs0,nmes0=nmes0,nbevt0=nbevt,nea0=nea0,
nwg0=nw,ncor0=ncor0,idiag0=idiag0,idtrunc0=idtrunc,
logspecif0=logspecif,npm0=NPM,fix0=fix0,nfix0=nfix,bfix0=bfix)
out <- list(conv=res$istop, V=res$v, best=res$b,
ppi=rep(NA,ns0*ng0), ppitest=rep(NA,ns0*ng0),
predRE=rep(NA,ns0*nea0), Yobs=rep(NA,nobs0),
resid_m=rep(NA,nobs0), resid_ss=rep(NA,nobs0),
marker=rep(NA,nsim), transfY=rep(NA,nsim),
pred_m_g=rep(NA,nobs0*ng0), pred_ss_g=rep(NA,nobs0*ng0),
time=rep(NA,nsim), risq_est=rep(NA,nsim*ng0*nbevt),
risqcum_est=rep(NA,nsim*ng0*nbevt),
statglob=NA, statevt=rep(NA,nbevt),
gconv=c(res$ca, res$cb, res$rdm), niter=res$ni,
loglik=res$fn.value)
}
if(out$conv %in% c(1,2,3))
{
estim0 <- 0
ll <- 0
ppi0 <- rep(0,ns0*ng0)
ppitest0 <- rep(0,ns0*ng0)
resid_m <- rep(0,nobs0)
resid_ss <- rep(0,nobs0)
pred_m_g <- rep(0,nobs0*ng0)
pred_ss_g <- rep(0,nobs0*ng0)
predRE <- rep(0,ns0*nea0)
marker <- rep(0,nsim)
transfY <- rep(0,nsim)
Yobs <- rep(0,nobs0)
risq_est <- rep(0,nsim*ng0*nbevt)
risqcum_est <- rep(0,nsim*ng0*nbevt)
statglob <- 0
statevt <- rep(0,nbevt)
estim0 <- 0
post <- .Fortran(C_loglikjointlcmm,
as.double(Y0),
as.double(X0),
as.integer(prior0),
as.double(pprior0),
as.double(tsurv0),
as.double(tsurv),
as.integer(devt),
as.integer(ind_survint),
as.integer(idprob0),
as.integer(idea0),
as.integer(idg0),
as.integer(idcor0),
as.integer(idcom),
as.integer(idspecif),
as.integer(idtdv),
as.integer(idlink),
as.double(epsY),
as.integer(nbzitr),
as.double(zitr0),
as.double(uniqueY0),
as.integer(nvalSPL0),
as.integer(indiceY0),
as.integer(typrisq),
as.integer(risqcom),
as.integer(nz),
as.double(zi),
as.integer(ns0),
as.integer(ng0),
as.integer(nv0),
as.integer(nobs0),
as.integer(nmes0),
as.integer(nbevt),
as.integer(nea0),
as.integer(nw),
as.integer(ncor0),
as.integer(idiag0),
as.integer(idtrunc),
as.integer(logspecif),
as.integer(NPM),
best=as.double(out$best),
ppi=as.double(ppi0),
ppitest=as.double(ppitest0),
resid_m=as.double(resid_m),
resid_ss=as.double(resid_ss),
pred_m_g=as.double(pred_m_g),
pred_ss_g=as.double(pred_ss_g),
predRE=as.double(predRE),
time=as.double(time),
risq_est=as.double(risq_est),
risqcum_est=as.double(risqcum_est),
marker=as.double(marker),
transfY=as.double(transfY),
as.integer(nsim),
Yobs=as.double(Yobs),
statglob=as.double(statglob),
statevt=as.double(statevt),
as.integer(fix0),
as.integer(nfix),
as.double(bfix),
as.integer(estim0),
loglik=as.double(ll))
out$ppi <- post$ppi
out$ppitest <- post$ppitest
out$predRE <- post$predRE
out$Yobs <- post$Yobs
out$resid_m <- post$resid_m
out$resid_ss <- post$resid_ss
out$marker <- post$marker
out$transfY <- post$transfY
out$pred_m_g <- post$pred_m_g
out$pred_ss_g <- post$pred_ss_g
out$risq_est <- post$risq_est
out$risqcum_est <- post$risqcum_est
out$statglob <- post$statglob
out$statevt <- post$statevt
}
## creer best a partir de b et bfix
best <- rep(NA,length(fix0))
best[which(fix0==0)] <- out$best
best[which(fix0==1)] <- bfix
names(best) <- names_best
out$best <- best
NPM <- NPM+nfix
### mettre NA pour les variances et covariances non calculees et 0 pr les prm fixes
if(length(posfix))
{
if(out$conv==3)
{
mr <- NPM-sum(pbH0)-length(posfix)
Vr <- matrix(0,mr,mr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
Vr <- t(Vr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
V <- matrix(NA,NPM,NPM)
V[setdiff(1:NPM,c(which(pbH0==1),posfix)),setdiff(1:NPM,c(which(pbH0==1),posfix))] <- Vr
V[,posfix] <- 0
V[posfix,] <- 0
V <- V[upper.tri(V,diag=TRUE)]
}
else
{
mr <- NPM-length(posfix)
Vr <- matrix(0,mr,mr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
Vr <- t(Vr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
V <- matrix(0,NPM,NPM)
V[setdiff(1:NPM,posfix),setdiff(1:NPM,posfix)] <- Vr
V <- V[upper.tri(V,diag=TRUE)]
}
}
else
{
if(out$conv==3)
{
mr <- NPM-sum(pbH0)
Vr <- matrix(0,mr,mr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
Vr <- t(Vr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
V <- matrix(NA,NPM,NPM)
V[setdiff(1:NPM,which(pbH0==1)),setdiff(1:NPM,which(pbH0==1))] <- Vr
V <- V[upper.tri(V,diag=TRUE)]
}
else
{
V <- out$V
}
}
### remplacer cholesky par varcov dans best
if(nvc>0)
{
ch <- out$best[nprob+nrisqtot+nvarxevt+nef+1:nvc]
if(idiag0==0)
{
U <- matrix(0,nea0,nea0)
U[upper.tri(U,diag=TRUE)] <- ch
z <- t(U) %*% U
out$best[nprob+nrisqtot+nvarxevt+nef+1:nvc] <- z[upper.tri(z,diag=TRUE)]
}
if(idiag0==1)
{
out$best[nprob+nrisqtot+nvarxevt+nef+1:nvc] <- out$best[nprob+nrisqtot+nvarxevt+nef+1:nvc]**2
}
}
else
{
ch <- NA
}
nom.unique[which(nom.unique=="(Intercept)")] <- "intercept"
### predictions des effets aleatoires
if (nea0>0)
{
predRE <- matrix(out$predRE,ncol=nea0,byrow=TRUE)
predRE <- data.frame(unique(IND),predRE)
colnames(predRE) <- c(subject,nom.unique[idea0!=0])
}
### proba des classes a posteroiri et classification
if(ng0>1)
{
ppi <- matrix(out$ppi,ncol=ng0,nrow=ns0,byrow=TRUE)
ppitest <- matrix(out$ppitest,ncol=ng0,byrow=TRUE)
}
else
{
ppi <- matrix(rep(1,ns0),ncol=ng0)
ppitest <- matrix(rep(1,ns0),ncol=ng0)
}
if(!(out$conv %in% c(1,2)))
{
classif <- rep(NA,ns0)
}
else
{
chooseClass <- function(ppi)
{
res <- which.max(ppi)
if(!length(res)) res <- NA
return(res)
}
classif <- apply(ppi,1,chooseClass)
}
ppi <- data.frame(unique(IND),classif,ppi)
temp <- paste("probYT",1:ng0,sep="")
colnames(ppi) <- c(subject,"class",temp)
rownames(ppi) <- 1:ns0
## faire pareil pour ppitest
if(!(out$conv %in% c(1,2)))
{
classif <- rep(NA,ns0)
}
else
{
classif <- apply(ppitest,1,chooseClass)
if(any(!is.finite(ppitest)))
{
classif <- rep(NA,ns0)
iok <- which(is.finite(ppitest[,1]))
classif[iok] <- apply(ppitest[iok,,drop=FALSE],1,which.max)
}
}
ppitest <- data.frame(unique(IND),classif,ppitest)
temp <- paste("probY",1:ng0,sep="")
colnames(ppitest) <- c(subject,"class",temp)
rownames(ppitest) <- 1:ns0
### predictions marginales et subject-specifiques
pred_m_g <- matrix(out$pred_m_g,nrow=nobs0,ncol=ng0)
pred_ss_g <- matrix(out$pred_ss_g,nrow=nobs0,ncol=ng0)
if((out$conv %in% c(1,2)))
{
pred_m <- out$Yobs-out$resid_m
pred_ss <- out$Yobs-out$resid_ss
}
else
{
pred_m <- rep(NA,nobs0)
pred_ss <- rep(NA,nobs0)
}
pred <- data.frame(IND,pred_m,out$resid_m,pred_ss,out$resid_ss,out$Yobs,pred_m_g,pred_ss_g)
temp<-paste("pred_m",1:ng0,sep="")
temp1<-paste("pred_ss",1:ng0,sep="")
colnames(pred)<-c(subject,"pred_m","resid_m","pred_ss","resid_ss","obs",temp,temp1)
if(!is.null(var.time))
{
pred <- data.frame(IND,pred_m,out$resid_m,pred_ss,out$resid_ss,out$Yobs,pred_m_g,pred_ss_g,timeobs)
colnames(pred)<-c(subject,"pred_m","resid_m","pred_ss","resid_ss","obs",temp,temp1,var.time)
}
### risques
risqcum_est <- matrix(out$risqcum_est,nrow=nsim,ncol=ng0*nbevt)
risq_est <- matrix(out$risq_est,nrow=nsim,ncol=ng0*nbevt)
predSurv <- cbind(time,risq_est,risqcum_est)
temp <- paste(paste("event",rep(1:nbevt,each=ng0),".RiskFct",sep=""),1:ng0,sep="")
temp1 <- paste(paste("event",rep(1:nbevt,each=ng0),".CumRiskFct",sep=""),1:ng0,sep="")
colnames(predSurv) <- c("time",temp,temp1)
rownames(predSurv) <- 1:nsim
###estimlink
estimlink <- cbind(out$marker,out$transfY)
colnames(estimlink) <- c(nomY,paste("transf",nomY,sep="."))
### score test
if(out$conv!=1) stats <- rep(NA,nbevt+1)
else
{
stats <- c(out$statglob,out$statevt)
stats[which((!is.finite(stats)) | (stats==9999))] <- NA
}
if(nbevt==1) stats <- stats[1]
### N
N <- NULL
N[1] <- nprob
N[2] <- nrisqtot
N[3] <- nvarxevt
N[4] <- nef
N[5] <- nvc
N[6] <- nw
N[7] <- ncor0
N[8] <- ntrtot0
N[9] <- nobs0
nevent <- rep(0,nbevt)
for(ke in 1:nbevt)
{
nevent[ke] <- length(which(devt==ke))
}
N <- c(N,nevent)
### noms des variables
Names <- list(Xnames=nom.unique,Xnames2=ttesLesVar,Yname=nomY,
ID=nom.subject,Tnames=noms.surv,prior.name=nom.prior,
TimeDepVar.name=nom.timedepvar, pprior.name=nom.pprior)
## levels = modalites des variables dans X0 (si facteurs)
levelsdata <- vector("list", length(ttesLesVar))
levelsfixed <- vector("list", length(ttesLesVar))
levelsrandom <- vector("list", length(ttesLesVar))
levelsmixture <- vector("list", length(ttesLesVar))
levelsclassmb <- vector("list", length(ttesLesVar))
levelssurv <- vector("list", length(ttesLesVar))
names(levelsdata) <- ttesLesVar
names(levelsfixed) <- ttesLesVar
names(levelsrandom) <- ttesLesVar
names(levelsmixture) <- ttesLesVar
names(levelsclassmb) <- ttesLesVar
names(levelssurv) <- ttesLesVar
for(v in ttesLesVar)
{
if(v == "intercept") next
if(is.factor(data[,v]))
{
levelsdata[[v]] <- levels(data[,v])
}
if(length(grep(paste("factor\\(",v,"\\)",sep=""), fixed)))
{
levelsfixed[[v]] <- levels(as.factor(data[,v]))
}
if(length(grep(paste("factor\\(",v,"\\)",sep=""), random)))
{
levelsrandom[[v]] <- levels(as.factor(data[,v]))
}
if(length(grep(paste("factor\\(",v,"\\)",sep=""), mixture)))
{
levelsmixture[[v]] <- levels(as.factor(data[,v]))
}
if(length(grep(paste("factor\\(",v,"\\)",sep=""), classmb)))
{
levelsclassmb[[v]] <- levels(as.factor(data[,v]))
}
if(length(grep(paste("factor\\(",v,"\\)",sep=""), form.surv)))
{
levelssurv[[v]] <- levels(as.factor(data[,v]))
}
}
levels <- list(levelsdata=levelsdata, levelsfixed=levelsfixed, levelsrandom=levelsrandom,
levelsmixture=levelsmixture, levelsclassmb=levelsclassmb,
levelssurv=levelssurv)
cost<-proc.time()-ptm
res <-list(ns=ns0,ng=ng0,idprob=idprob0,idcom=idcom,
idspecif=idspecif,idtdv=idtdv,idg=idg0,idea=idea0,
idcor=idcor0,loglik=out$loglik,best=out$best,V=V,
gconv=out$gconv,conv=out$conv,call=cl,niter=out$niter,
N=N,idiag=idiag0,pred=pred,pprob=ppi,pprobY=ppitest,
predRE=predRE,Names=Names,cholesky=ch,logspecif=logspecif,
estimlink=estimlink,epsY=epsY,linktype=idlink,linknodes=zitr0,
predSurv=predSurv,hazard=list(typrisq,hazardtype,zi,nz),
scoretest=stats,na.action=linesNA,
AIC=2*(length(out$best)-length(posfix)-out$loglik),
BIC=(length(out$best)-length(posfix))*log(ns0)-2*out$loglik,
data=datareturn,levels=levels,var.time=var.time,runtime=cost[3])
class(res) <-c("Jointlcmm")
if(verbose==TRUE) cat("The program took", round(cost[3],2), "seconds \n")
return(res)
}
#' @rdname Jointlcmm
#' @export
jlcmm <- Jointlcmm
CecileProust-Lima/lcmm documentation built on March 3, 2024, 5:23 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions Jointlcmm
Documented in Jointlcmm
#' Estimation of joint latent class models for longitudinal and time-to-event
#' data
#'
#' This function fits joint latent class mixed models for a longitudinal
#' outcome and a right-censored (possibly left-truncated) time-to-event. The
#' function handles competing risks and Gaussian or non Gaussian (curvilinear)
#' longitudinal outcomes. For curvilinear longitudinal outcomes, normalizing
#' continuous functions (splines or Beta CDF) can be specified as in
#' \code{lcmm}.
#'
#'
#' A. BASELINE RISK FUNCTIONS
#'
#' For the baseline risk functions, the following parameterizations were
#' considered. Be careful, parametrisations changed in lcmm_V1.5:
#'
#' 1. With the "Weibull" function: 2 parameters are necessary w_1 and w_2 so
#' that the baseline risk function a_0(t) = w_1^2*w_2^2*(w_1^2*t)^(w_2^2-1) if
#' logscale=FALSE and a_0(t) = exp(w_1)*exp(w_2)(t)^(exp(w_2)-1) if
#' logscale=TRUE.
#'
#' 2. with the "piecewise" step function and nz nodes (y_1,...y_nz), nz-1
#' parameters are necesssary p_1,...p_nz-1 so that the baseline risk function
#' a_0(t) = p_j^2 for y_j < t =< y_j+1 if logscale=FALSE and a_0(t) = exp(p_j)
#' for y_j < t =< y_j+1 if logscale=TRUE.
#'
#' 3. with the "splines" function and nz nodes (y_1,...y_nz), nz+2 parameters
#' are necessary s_1,...s_nz+2 so that the baseline risk function a_0(t) =
#' sum_j s_j^2 M_j(t) if logscale=FALSE and a_0(t) = sum_j exp(s_j) M_j(t) if
#' logscale=TRUE where {M_j} is the basis of cubic M-splines.
#'
#' Two parametrizations of the baseline risk function are proposed
#' (logscale=TRUE or FALSE) because in some cases, especially when the
#' instantaneous risks are very close to 0, some convergence problems may
#' appear with one parameterization or the other. As a consequence, we
#' recommend to try the alternative parameterization (changing logscale option)
#' when a joint latent class model does not converge (maximum number of
#' iterations reached) where as convergence criteria based on the parameters
#' and likelihood are small.
#'
#' B. THE VECTOR OF PARAMETERS B
#'
#' The parameters in the vector of initial values \code{B} or in the vector of
#' maximum likelihood estimates \code{best} are included in the following
#' order: (1) ng-1 parameters are required for intercepts in the latent class
#' membership model, and if covariates are included in \code{classmb}, ng-1
#' parameters should be entered for each one; (2) parameters for the baseline
#' risk function: 2 parameters for each Weibull, nz-1 for each piecewise
#' constant risk and nz+2 for each splines risk; this number should be
#' multiplied by ng if specific hazard is specified; otherwise, ng-1 additional
#' proportional effects are expected if PH hazard is specified; otherwise
#' nothing is added if common hazard is specified. In the presence of competing
#' events, the number of parameters should be adapted to the number of causes
#' of event; (3) for all covariates in \code{survival}, ng parameters are
#' required if the covariate is inside a \code{mixture()}, otherwise 1
#' parameter is required. Covariates parameters should be included in the same
#' order as in \code{survival}. In the presence of cause-specific effects, the
#' number of parameters should be multiplied by the number of causes; (4) for
#' all covariates in \code{fixed}, one parameter is required if the covariate
#' is not in \code{mixture}, ng parameters are required if the covariate is
#' also in \code{mixture}. Parameters should be included in the same order as
#' in \code{fixed}; (5) the variance of each random-effect specified in
#' \code{random} (including the intercept) if \code{idiag=TRUE} and the
#' inferior triangular variance-covariance matrix of all the random-effects if
#' \code{idiag=FALSE}; (6) only if \code{nwg=TRUE}, ng-1 parameters for
#' class-specific proportional coefficients for the variance covariance matrix
#' of the random-effects; (7) the variance of the residual error.
#'
#' C. CAUTION
#'
#' Some caution should be made when using the program:
#'
#' (1) As the log-likelihood of a latent class model can have multiple maxima,
#' a careful choice of the initial values is crucial for ensuring convergence
#' toward the global maximum. The program can be run without entering the
#' vector of initial values (see point 2). However, we recommend to
#' systematically enter initial values in \code{B} and try different sets of
#' initial values.
#'
#' (2) The automatic choice of initial values that we provide requires the
#' estimation of a preliminary linear mixed model. The user should be aware
#' that first, this preliminary analysis can take time for large datatsets and
#' second, that the generated initial values can be very not likely and even
#' may converge slowly to a local maximum. This is a reason why several
#' alternatives exist. The vector of initial values can be directly specified
#' in \code{B} the initial values can be generated (automatically or randomly)
#' from a model with \code{ng=}. Finally, function \code{gridsearch} performs
#' an automatic grid search.
#'
#' (3) Convergence criteria are very strict as they are based on derivatives of
#' the log-likelihood in addition to the parameter and log-likelihood
#' stability. In some cases, the program may not converge and reach the
#' maximum number of iterations fixed at 150. In this case, the user should
#' check that parameter estimates at the last iteration are not on the
#' boundaries of the parameter space. If the parameters are on the boundaries
#' of the parameter space, the identifiability of the model is critical. This
#' may happen especially when baseline risk functions involve splines (value
#' close to the lower boundary - 0 with logscale=F -infinity with logscale=F)
#' or classmb parameters that are too high or low (perfect classification) or
#' linkfunction parameters. When identifiability of some parameters is
#' suspected, the program can be run again from the former estimates by fixing
#' the suspected parameters to their value with option posfix. This usually
#' solves the problem. An alternative is to remove the parameters of the Beta
#' of Splines link function from the inverse of the Hessian with option
#' partialH. If not, the program should be run again with other initial
#' values. Some problems of convergence may happen when the instantaneous
#' risks of event are very low and "piecewise" or "splines" baseline risk
#' functions are specified. In this case, changing the parameterization of the
#' baseline risk functions with option logscale is recommended (see paragraph A
#' for details).
#'
#' @aliases Jointlcmm jlcmm
#' @param fixed two-sided linear formula object for the fixed-effects in the
#' linear mixed model. The response outcome is on the left of \code{~} and the
#' covariates are separated by \code{+} on the right of the \code{~}. By
#' default, an intercept is included. If no intercept, \code{-1} should be the
#' first term included on the right of \code{~}.
#' @param mixture one-sided formula object for the class-specific fixed effects
#' in the linear mixed model (to specify only for a number of latent classes
#' greater than 1). Among the list of covariates included in \code{fixed}, the
#' covariates with class-specific regression parameters are entered in
#' \code{mixture} separated by \code{+}. By default, an intercept is included.
#' If no intercept, \code{-1} should be the first term included.
#' @param random optional one-sided formula for the random-effects in the
#' linear mixed model. Covariates with a random-effect are separated by
#' \code{+}. By default, an intercept is included. If no intercept, \code{-1}
#' should be the first term included.
#' @param subject name of the covariate representing the grouping structure
#' (called subject identifier) specified with ''.
#' @param classmb optional one-sided formula describing the covariates in the
#' class-membership multinomial logistic model. Covariates included are
#' separated by \code{+}. No intercept should be included in this formula.
#' @param ng optional number of latent classes considered. If \code{ng=1} (by
#' default) no \code{mixture} nor \code{classmb} should be specified. If
#' \code{ng>1}, \code{mixture} is required.
#' @param idiag optional logical for the structure of the variance-covariance
#' matrix of the random-effects. If \code{FALSE}, a non structured matrix of
#' variance-covariance is considered (by default). If \code{TRUE} a diagonal
#' matrix of variance-covariance is considered.
#' @param nwg optional logical indicating if the variance-covariance of the
#' random-effects is class-specific. If \code{FALSE} the variance-covariance
#' matrix is common over latent classes (by default). If \code{TRUE} a
#' class-specific proportional parameter multiplies the variance-covariance
#' matrix in each class (the proportional parameter in the last latent class
#' equals 1 to ensure identifiability).
#' @param survival two-sided formula object. The left side of the formula
#' corresponds to a \code{surv()} object of type "counting" for right-censored
#' and left-truncated data (example: \code{Surv(Time,EntryTime,Indicator)}) or
#' of type "right" for right-censored data (example:
#' \code{Surv(Time,Indicator)}). Multiple causes of event can be considered in
#' the Indicator (0 for censored, k for cause k of event). The right side of
#' the formula specifies the names of covariates to include in the survival
#' model with \code{mixture()} when the effect is class-specific (example:
#' \code{Surv(Time,Indicator) ~} \code{ X1 + mixture(X2)} for a class-common
#' effect of X1 and a class-specific effect of X2). In the presence of
#' competing events, covariate effects are common by default. Code
#' \code{cause(X3)} specifies a cause-specific covariate effect for X3 on each
#' cause of event while \code{cause1(X3)} (or \code{cause2(X3)}, ...) specifies
#' a cause-specific effect of X3 on the first (or second, ...) cause only.
#' @param hazard optional family of hazard function assumed for the survival
#' model. By default, "Weibull" specifies a Weibull baseline risk function.
#' Other possibilities are "piecewise" for a piecewise constant risk function
#' or "splines" for a cubic M-splines baseline risk function. For these two
#' latter families, the number of nodes and the location of the nodes should be
#' specified as well, separated by \code{-}. The number of nodes is entered
#' first followed by \code{-}, then the location is specified with "equi",
#' "quant" or "manual" for respectively equidistant nodes, nodes at quantiles
#' of the times of event distribution or interior nodes entered manually in
#' argument \code{hazardnodes}. It is followed by \code{-} and finally
#' "piecewise" or "splines" indicates the family of baseline risk function
#' considered. Examples include "5-equi-splines" for M-splines with 5
#' equidistant nodes, "6-quant-piecewise" for piecewise constant risk over 5
#' intervals and nodes defined at the quantiles of the times of events
#' distribution and "9-manual-splines" for M-splines risk function with 9
#' nodes, the vector of 7 interior nodes being entered in the argument
#' \code{hazardnodes}. In the presence of competing events, a vector of hazards
#' should be provided such as \code{hazard=c("Weibull","splines"} with 2 causes
#' of event, the first one modelled by a Weibull baseline cause-specific risk
#' function and the second one by splines.
#' @param hazardtype optional indicator for the type of baseline risk function
#' when ng>1. By default "Specific" indicates a class-specific baseline risk
#' function. Other possibilities are "PH" for a baseline risk function
#' proportional in each latent class, and "Common" for a baseline risk function
#' that is common over classes. In the presence of competing events, a vector
#' of hazardtypes should be given.
#' @param hazardnodes optional vector containing interior nodes if
#' \code{splines} or \code{piecewise} is specified for the baseline hazard
#' function in \code{hazard}.
#' @param hazardrange optional vector indicating the range of the survival times (that is
#' the minimum and maximum). By default, the range is defined according to the
#' minimum and maximum observed values of the survival times. The option should be
#' used only for piecewise constant and Splines hazard functions.
#' @param TimeDepVar optional vector containing an intermediate time
#' corresponding to a change in the risk of event. This time-dependent
#' covariate can only take the form of a time variable with the assumption that
#' there is no effect on the risk before this time and a constant effect on the
#' risk of event after this time (example: initiation of a treatment to account
#' for).
#' @param link optional family of link functions to estimate. By default,
#' "linear" option specifies a linear link function leading to a standard
#' linear mixed model (homogeneous or heterogeneous as estimated in
#' \code{hlme}). Other possibilities include "beta" for estimating a link
#' function from the family of Beta cumulative distribution functions,
#' "thresholds" for using a threshold model to describe the correspondence
#' between each level of an ordinal outcome and the underlying latent process,
#' and "Splines" for approximating the link function by I-splines. For this
#' latter case, the number of nodes and the nodes location should be also
#' specified. The number of nodes is first entered followed by \code{-}, then
#' the location is specified with "equi", "quant" or "manual" for respectively
#' equidistant nodes, nodes at quantiles of the marker distribution or interior
#' nodes entered manually in argument \code{intnodes}. It is followed by
#' \code{-} and finally "splines" is indicated. For example, "7-equi-splines"
#' means I-splines with 7 equidistant nodes, "6-quant-splines" means I-splines
#' with 6 nodes located at the quantiles of the marker distribution and
#' "9-manual-splines" means I-splines with 9 nodes, the vector of 7 interior
#' nodes being entered in the argument \code{intnodes}.
#' @param intnodes optional vector of interior nodes. This argument is only
#' required for a I-splines link function with nodes entered manually.
#' @param epsY optional definite positive real used to rescale the marker in
#' (0,1) when the beta link function is used. By default, epsY=0.5.
#' @param range optional vector indicating the range of the outcome (that is
#' the minimum and maximum). By default, the range is defined according to the
#' minimum and maximum observed values of the outcome. The option should be
#' used only for Beta and Splines transformations.
#' @param cor optional brownian motion or autoregressive process modeling the
#' correlation between the observations. "BM" or "AR" should be specified,
#' followed by the time variable between brackets. By default, no correlation
#' is added.
#' @param data optional data frame containing the variables named in
#' \code{fixed}, \code{mixture}, \code{random}, \code{classmb} and
#' \code{subject}.
#' @param B optional specification for the initial values for the parameters.
#' Three options are allowed: (1) a vector of initial values is entered (the
#' order in which the parameters are included is detailed in \code{details}
#' section). (2) nothing is specified. A preliminary analysis involving the
#' estimation of a standard linear mixed model is performed to choose initial
#' values. (3) when ng>1, a Jointlcmm object is entered. It should correspond
#' to the exact same structure of model but with ng=1. The program will
#' automatically generate initial values from this model. This specification
#' avoids the preliminary analysis indicated in (2) Note that due to possible
#' local maxima, the \code{B} vector should be specified and several different
#' starting points should be tried.
#' @param convB optional threshold for the convergence criterion based on the
#' parameter stability. By default, convB=0.0001.
#' @param convL optional threshold for the convergence criterion based on the
#' log-likelihood stability. By default, convL=0.0001.
#' @param convG optional threshold for the convergence criterion based on the
#' derivatives. By default, convG=0.0001.
#' @param maxiter optional maximum number of iterations for the Marquardt
#' iterative algorithm. By default, maxiter=150.
#' @param nsim optional number of points for the predicted survival curves and
#' predicted baseline risk curves. By default, nsim=100.
#' @param prior optional name of a covariate containing a prior information
#' about the latent class membership. The covariate should be an integer with
#' values in 0,1,...,ng. Value O indicates no prior for the subject while a
#' value in 1,...,ng indicates that the subject belongs to the corresponding
#' latent class.
#' @param pprior optional vector specifying the names of the covariates containing the
#' prior probabilities to belong to each latent class. These probabilities should be
#' between 0 and 1 and should sum up to 1 for each subject.
#' @param logscale optional boolean indicating whether an exponential
#' (logscale=TRUE) or a square (logscale=FALSE -by default) transformation is
#' used to ensure positivity of parameters in the baseline risk functions. See
#' details section
#' @param subset a specification of the rows to be used: defaults to all rows.
#' This can be any valid indexing vector for the rows of data or if that is not
#' supplied, a data frame made up of the variable used in formula.
#' @param na.action Integer indicating how NAs are managed. The default is 1
#' for 'na.omit'. The alternative is 2 for 'na.fail'. Other options such as
#' 'na.pass' or 'na.exclude' are not implemented in the current version.
#' @param posfix Optional vector specifying the indices in vector B of the
#' parameters that should not be estimated. Default to NULL, all parameters are
#' estimated.
#' @param partialH optional logical for Piecewise and Splines baseline risk
#' functions and Splines link functions only. Indicates whether the parameters of the
#' baseline risk or link functions can be dropped from the Hessian matrix to
#' define convergence criteria.
#' @param verbose logical indicating if information about computation should be
#' reported. Default to TRUE.
#' @param returndata logical indicating if data used for computation should be
#' returned. Default to FALSE, data are not returned.
#' @param var.time optional character indicating the name of the time variable.
#' @param nproc the number cores for parallel computation.
#' Default to 1 (sequential mode).
#' @param clustertype optional character indicating the type of cluster for parallel computation.
#' @return The list returned is: \item{loglik}{log-likelihood of the model}
#' \item{best}{vector of parameter estimates in the same order as specified in
#' \code{B} and detailed in section \code{details}}
#' \item{V}{if the model converged (conv=1 or 3), vector containing
#' the upper triangle matrix of variance-covariance estimates of \code{Best}
#' with exception for variance-covariance parameters of the random-effects for
#' which \code{V} contains the variance-covariance estimates of the Cholesky
#' transformed parameters displayed in \code{cholesky}.
#' If conv=2, \code{V} contains the second derivatives of the log-likelihood.}
#' \item{gconv}{vector of
#' convergence criteria: 1. on the parameters, 2. on the likelihood, 3. on the
#' derivatives} \item{conv}{status of convergence: =1 if the convergence
#' criteria were satisfied, =2 if the maximum number of iterations was reached,
#' =4 or 5 if a problem occured during optimisation} \item{call}{the matched
#' call} \item{niter}{number of Marquardt iterations} \item{pred}{table of
#' individual predictions and residuals; it includes marginal predictions
#' (pred_m), marginal residuals (resid_m), subject-specific predictions
#' (pred_ss) and subject-specific residuals (resid_ss) averaged over classes,
#' the observation (obs) and finally the class-specific marginal and
#' subject-specific predictions (with the number of the latent class:
#' pred_m_1,pred_m_2,...,pred_ss_1,pred_ss_2,...). If \code{var.time}
#' is specified, the corresponding measurement time is also included.}
#' \item{pprob}{table of
#' posterior classification and posterior individual class-membership
#' probabilities based on the longitudinal data and the time-to-event data}
#' \item{pprobY}{table of posterior classification and posterior individual
#' class-membership probabilities based only on the longitudinal data}
#' \item{predRE}{table containing individual predictions of the random-effects:
#' a column per random-effect, a line per subject} \item{cholesky}{vector
#' containing the estimates of the Cholesky transformed parameters of the
#' variance-covariance matrix of the random-effects} \item{scoretest}{Statistic
#' of the Score Test for the conditional independence assumption of the
#' longitudinal and survival data given the latent class structure. Under the
#' null hypothesis, the statistics is a Chi-square with p degrees of freedom
#' where p indicates the number of random-effects in the longitudinal mixed
#' model. See Jacqmin-Gadda and Proust-Lima (2009) for more details.}
#' \item{predSurv}{table of predictions giving for the window of times to event
#' (called "time"), the predicted baseline risk function in each latent class
#' (called "RiskFct") and the predicted cumulative baseline risk function in
#' each latent class (called "CumRiskFct").} \item{hazard}{internal information
#' about the hazard specification used in related functions} \item{data}{the
#' original data set (if returndata is TRUE)}
#' %\item{specif}{internal information used in related functions}
#' %\item{Names}{internal information used in related fnctions}
#' %\item{Names2}{internal information used in related functions}
#' @author Cecile Proust Lima, Amadou Diakite and Viviane Philipps
#'
#' \email{cecile.proust-lima@@inserm.fr}
#' @seealso \code{\link{postprob}}, \code{\link{plot.Jointlcmm}},
#' \code{\link{plot.predict}}, \code{\link{epoce}}
#' @references
#'
#' Proust-Lima C, Philipps V, Liquet B (2017). Estimation of Extended Mixed
#' Models Using Latent Classes and Latent Processes: The R Package lcmm.
#' Journal of Statistical Software, 78(2), 1-56. doi:10.18637/jss.v078.i02
#'
#'
#' Lin, H., Turnbull, B. W., McCulloch, C. E. and Slate, E. H. (2002). Latent
#' class models for joint analysis of longitudinal biomarker and event process
#' data: application to longitudinal prostate-specific antigen readings and
#' prostate cancer. Journal of the American Statistical Association 97, 53-65.
#'
#' Proust-Lima, C. and Taylor, J. (2009). Development and validation of a
#' dynamic prognostic tool for prostate cancer recurrence using repeated
#' measures of post-treatment PSA: a joint modelling approach. Biostatistics
#' 10, 535-49.
#'
#' Jacqmin-Gadda, H. and Proust-Lima, C. (2010). Score test for conditional
#' independence between longitudinal outcome and time-to-event given the
#' classes in the joint latent class model. Biometrics 66(1), 11-9
#'
#' Proust-Lima, Sene, Taylor and Jacqmin-Gadda (2014). Joint latent class
#' models of longitudinal and time-to-event data: a review. Statistical Methods
#' in Medical Research 23, 74-90.
#' @examples
#'
#' \dontrun{
#' #### Example of a joint latent class model estimated for a varying number
#' # of latent classes:
#' # The linear mixed model includes a subject- (ID) and class-specific
#' # linear trend (intercept and Time in fixed, random and mixture components)
#' # and a common effect of X1 and its interaction with time over classes
#' # (in fixed).
#' # The variance of the random intercept and slopes are assumed to be equal
#' # over classes (nwg=F).
#' # The covariate X3 predicts the class membership (in classmb).
#' # The baseline hazard function is modelled with cubic M-splines -3
#' # nodes at the quantiles- (in hazard) and a proportional hazard over
#' # classes is assumed (in hazardtype). Covariates X1 and X2 predict the
#' # risk of event (in survival) with a common effect over classes for X1
#' # and a class-specific effect of X2.
#' # !CAUTION: for illustration, only default initial values where used but
#' # other sets of initial values should be tried to ensure convergence
#' # towards the global maximum.
#'
#'
#' #### estimation with 1 latent class (ng=1): independent models for the
#' # longitudinal outcome and the time of event
#' m1 <- Jointlcmm(fixed= Ydep1~X1*Time,random=~Time,subject='ID',
#' survival = Surv(Tevent,Event)~ X1+X2 ,hazard="3-quant-splines",
#' hazardtype="PH",ng=1,data=data_lcmm)
#' summary(m1)
#' #Goodness-of-fit statistics for m1:
#' # maximum log-likelihood: -3944.77 ; AIC: 7919.54 ; BIC: 7975.09
#' }
#'
#' #### estimation with 2 latent classes (ng=2)
#' m2 <- Jointlcmm(fixed= Ydep1~Time*X1,mixture=~Time,random=~Time,
#' classmb=~X3,subject='ID',survival = Surv(Tevent,Event)~X1+mixture(X2),
#' hazard="3-quant-splines",hazardtype="PH",ng=2,data=data_lcmm,
#' B=c(0.64,-0.62,0,0,0.52,0.81,0.41,0.78,0.1,0.77,-0.05,10.43,11.3,-2.6,
#' -0.52,1.41,-0.05,0.91,0.05,0.21,1.5))
#' summary(m2)
#' #Goodness-of-fit statistics for m2:
#' # maximum log-likelihood: -3921.27; AIC: 7884.54; BIC: 7962.32
#'
#' \dontrun{
#' #### estimation with 3 latent classes (ng=3)
#' m3 <- Jointlcmm(fixed= Ydep1~Time*X1,mixture=~Time,random=~Time,
#' classmb=~X3,subject='ID',survival = Surv(Tevent,Event)~ X1+mixture(X2),
#' hazard="3-quant-splines",hazardtype="PH",ng=3,data=data_lcmm,
#' B=c(0.77,0.4,-0.82,-0.27,0,0,0,0.3,0.62,2.62,5.31,-0.03,1.36,0.82,
#' -13.5,10.17,10.24,11.51,-2.62,-0.43,-0.61,1.47,-0.04,0.85,0.04,0.26,1.5))
#' summary(m3)
#' #Goodness-of-fit statistics for m3:
#' # maximum log-likelihood: -3890.26 ; AIC: 7834.53; BIC: 7934.53
#'
#' #### estimation with 4 latent classes (ng=4)
#' m4 <- Jointlcmm(fixed= Ydep1~Time*X1,mixture=~Time,random=~Time,
#' classmb=~X3,subject='ID',survival = Surv(Tevent,Event)~ X1+mixture(X2),
#' hazard="3-quant-splines",hazardtype="PH",ng=4,data=data_lcmm,
#' B=c(0.54,-0.42,0.36,-0.94,-0.64,-0.28,0,0,0,0.34,0.59,2.6,2.56,5.26,
#' -0.1,1.27,1.34,0.7,-5.72,10.54,9.02,10.2,11.58,-2.47,-2.78,-0.28,-0.57,
#' 1.48,-0.06,0.61,-0.07,0.31,1.5))
#' summary(m4)
#' #Goodness-of-fit statistics for m4:
#' # maximum log-likelihood: -3886.93 ; AIC: 7839.86; BIC: 7962.09
#'
#'
#' ##### The model with 3 latent classes is retained according to the BIC
#' ##### and the conditional independence assumption is not rejected at
#' ##### the 5% level.
#' # posterior classification
#' plot(m3,which="postprob")
#' # Class-specific predicted baseline risk & survival functions in the
#' # 3-class model retained (for the reference value of the covariates)
#' plot(m3,which="baselinerisk",bty="l")
#' plot(m3,which="baselinerisk",ylim=c(0,5),bty="l")
#' plot(m3,which="survival",bty="l")
#' # class-specific predicted trajectories in the 3-class model retained
#' # (with characteristics of subject ID=193)
#' data <- data_lcmm[data_lcmm$ID==193,]
#' plot(predictY(m3,var.time="Time",newdata=data,bty="l"))
#' # predictive accuracy of the model evaluated with EPOCE
#' vect <- 1:15
#' cvpl <- epoce(m3,var.time="Time",pred.times=vect)
#' summary(cvpl)
#' plot(cvpl,bty="l",ylim=c(0,2))
#' ############## end of example ##############
#' }
#'
#' @export
#'
#'
Jointlcmm <- function(fixed,mixture,random,subject,classmb,ng=1,idiag=FALSE,nwg=FALSE,
survival,hazard="Weibull",hazardtype="Specific",
hazardnodes=NULL,hazardrange=NULL,TimeDepVar=NULL,
link=NULL,intnodes=NULL,epsY=0.5,range=NULL,
cor=NULL,data,B,convB=0.0001,convL=0.0001,convG=0.0001,maxiter=100,
nsim=100,prior=NULL,pprior=NULL,logscale=FALSE,subset=NULL,na.action=1,
posfix=NULL,partialH=FALSE,verbose=FALSE,returndata=FALSE,
var.time=NULL,nproc=1,clustertype=NULL)
{
ptm<-proc.time()
cl <- match.call()
if(!missing(mixture) & ng==1) stop("No mixture can be specified with ng=1")
if(missing(mixture) & ng>1) stop("The argument mixture has to be specified for ng > 1")
if(!missing(classmb) & ng==1) stop("No classmb can be specified with ng=1")
if(missing(fixed)) stop("The argument Fixed must be specified in any model")
if(missing(random)) random <- ~-1
if(missing(classmb) & ng==1) classmb <- ~-1
if(missing(classmb) & ng>1) classmb <- ~1
if(missing(mixture)) mixture <- ~-1
if(ng==1 & nwg==TRUE) stop ("The argument nwg should be FALSE for ng=1")
if(ng==1) hazardtype <- "Specific" #??
if(any(!(hazardtype %in% c("Common","Specific","PH")))) stop("'hazardtype' should be either 'Common' or 'Specific' or 'PH'")
if(!inherits(fixed,"formula")) stop("The argument fixed must be a formula")
if(!inherits(mixture,"formula")) stop("The argument mixture must be a formula")
if(!inherits(random,"formula")) stop("The argument random must be a formula")
if(!inherits(classmb,"formula")) stop("The argument classmb must be a formula")
if(missing(data)){ stop("The argument data should be specified and defined as a data.frame")}
if(nrow(data)==0) stop("Data should not be empty")
if(missing(subject)){ stop("The argument subject must be specified in any model even without random-effects")}
if(!is.numeric(data[[subject]])) stop("The argument subject must be numeric")
if(is.null(link)) link <- "NULL"
if(any(link=="thresholds")) stop("The link function thresholds is not available yet")
if(link %in% c("linear","beta") & !is.null(intnodes)) stop("Intnodes should only be specified with splines links")
nom.subject <- as.character(subject)
if(!isTRUE(nom.subject %in% colnames(data))) stop(paste("Data should contain variable",nom.subject))
nom.prior <- NULL
if(!is.null(prior))
{
nom.prior <- as.character(prior)
if(!isTRUE(nom.prior %in% colnames(data))) stop(paste("Data should contain variable",nom.prior))
}
nom.pprior <- NULL
if(!is.null(pprior))
{
if(ng==1) stop("pprior is only useful if ng>1")
if(!is.character(pprior)) stop("pprior should be a character vector")
if(length(pprior) != ng) stop("pprior should be of length ng")
nom.pprior <- as.character(pprior)
if(!all(nom.pprior %in% colnames(data))) stop(paste("Data should contain variable",paste(nom.pprior, collapse=" ")))
}
nom.timedepvar <- NULL
if(!missing(TimeDepVar))
{
if(!is.null(TimeDepVar))
{
nom.timedepvar <- as.character(TimeDepVar)
if(!isTRUE(nom.timedepvar %in% colnames(data))) stop(paste("Data should contain variable",nom.timedepvar))
}
}
if(!(na.action %in% c(1,2))) stop("only 1 for 'na.omit' or 2 for 'na.fail' are required in na.action argument")
#if(length(posfix) & missing(B)) stop("A set of initial parameters must be specified if some parameters are not estimated")
## garder data tel quel pour le renvoyer
if(returndata==TRUE)
{
datareturn <- data
}
else
{
datareturn <- NULL
}
if(!isTRUE(all.equal(as.character(cl$subset),character(0))))
{
cc <- cl
cc <- cc[c(1,which(names(cl)=="subset"))]
cc[[1]] <- as.name("model.frame")
cc$formula <- formula(paste("~",paste(colnames(data),collapse="+")))
cc$data <- data
cc$na.action <- na.pass
data <- eval(cc)
}
attributes(data)$terms <- NULL
### test de l'argument cor
ncor0 <- 0
cor.type <- cl$cor[1]
cor.time <- cl$cor[2]
cor <- paste(cor.type,cor.time,sep="-")
if (!isTRUE(all.equal(cor,character(0))))
{
if (substr(cor,1,2)=="AR") { ncor0 <- 2 }
else if (substr(cor,1,2)=="BM") { ncor0 <- 1 }
else { stop("The argument cor must be of type AR or BM") }
if(!(strsplit(cor,"-")[[1]][2] %in% colnames(data))) stop("Unable to find time variable from argument 'cor' in 'data'")
else { cor.var.time <- strsplit(cor,"-")[[1]][2] }
}
### fin test argument cor
## objet Surv
surv <- cl$survival[[2]]
if(length(surv)==3) #censure droite sans troncature gauche
{
idtrunc <- 0
Tevent <- getElement(object=data,name=as.character(surv[2]))
Event <- getElement(object=data,name=as.character(surv[3]))
Tentry <- rep(0,length(Tevent)) #si pas de troncature, Tentry=0
noms.surv <- c(as.character(surv[2]),as.character(surv[3]))
surv <- do.call("Surv",list(time=Tevent,event=Event,type="mstate"))
}
else
{
if(length(surv)==4) #censure droite et troncature
{
idtrunc <- 1
Tentry <- getElement(object=data,name=as.character(surv[2]))
Tevent <- getElement(object=data,name=as.character(surv[3]))
Event <- getElement(object=data,name=as.character(surv[4]))
noms.surv <- c(as.character(surv[2]),as.character(surv[3]),as.character(surv[4]))
surv <- do.call("Surv",list(time=Tentry,time2=Tevent,event=Event,type="mstate"))
}
}
## nombre d'evenement concurrents
nbevt <- length(attr(surv,"states"))
if(nbevt<1)
{
if(maxiter != 0) stop("No observed event in the data")
nbevt <- 1
}
##pour acces aux attributs des formules
afixed <- terms(fixed)
if(link!="NULL" & attr(afixed,"intercept")==0) stop("An intercept should appear in fixed for identifiability purposes")
amixture <- terms(mixture)
arandom <- terms(random)
aclassmb <- terms(classmb)
## pour la formule pour survivial, creer 3 formules :
## une pour les covariables en mixture, une pour les covariables avec effet specifique a la cause, et une pour les effets communs.
form.surv <- cl$survival[3]
noms.form.surv <- all.vars(attr(terms(formula(paste("~",form.surv))),"variables"))
if(length(noms.form.surv)==0)
{
form.cause <- ~-1
form.causek <- vector("list",nbevt)
for(k in 1:nbevt) form.causek[[k]] <- ~-1
form.mixture <- ~-1
form.commun <- ~-1
asurv <- terms(~-1)
}
else
{
##creer la formula pour cause
form1 <- gsub("mixture","",form.surv)
form1 <- formula(paste("~",form1))
asurv1 <- terms(form1,specials="cause")
ind.cause <- attr(asurv1,"specials")$cause
if(length(ind.cause))
{
form.cause <- paste(labels(asurv1)[ind.cause],collapse="+")
form.cause <- gsub("cause","",form.cause)
form.cause <- formula(paste("~",form.cause))
}
else
{
form.cause <- ~-1
}
## formules pour causek
form.causek <- vector("list",nbevt)
for(k in 1:nbevt)
{
formk <- gsub("mixture","",form.surv)
for(kk in 1:nbevt)
{
if(kk != k) formk <- gsub(paste("cause",kk,sep=""),"",formk)
}
asurvk <- terms(formula(paste("~",formk)),specials=paste("cause",k,sep=""))
ind.causek <- attr(asurvk,"specials")$cause
if(length(ind.causek))
{
formcausek <- paste(labels(asurvk)[ind.causek],collapse="+")
formcausek <- gsub(paste("cause",k,sep=""),"",formcausek)
formcausek <- formula(paste("~",formcausek))
form.causek[[k]] <- formcausek
}
else
{
form.causek[[k]] <- ~-1
}
}
##creer la formule pour mixture
form2 <- form.surv
for( k in 1:nbevt)
{
form2 <- gsub(paste("cause",k,sep=""),"",form2)
}
form2 <- gsub("cause","",form2)
form2 <- formula(paste("~",form2))
asurv2 <- terms(form2,specials="mixture")
ind.mixture <- attr(asurv2,"specials")$mixture
if(length(ind.mixture))
{
form.mixture <- paste(labels(asurv2)[ind.mixture],collapse="+")
form.mixture <- gsub("mixture","",form.mixture)
form.mixture <- formula(paste("~",form.mixture))
}
else
{
form.mixture <- ~-1
}
## creer la formule pour ni cause ni mixture
asurv <- terms(formula(paste("~",form.surv)),specials=c("cause","mixture",paste("cause",1:nbevt,sep="")))
ind.commun <- setdiff(1:length(labels(asurv)),unlist(attr(asurv,"specials")))
if(length(ind.commun))
{
form.commun <- paste(labels(asurv)[ind.commun],collapse="+")
form.commun <- gsub("mixture","",form.commun) #si X1*mixture(X2), alors X1:mixture(X2) dans form.commun
form.commun <- gsub("cause","",form.commun) # si X1:cause(X2)
form.commun <- formula(paste("~",form.commun))
##NB: si mixture(X1)*cause(X2), X1:X2 en commun
}
else
{
form.commun <- ~-1
}
}
##verifier si toutes les variables sont dans data
variables <- c(attr(afixed,"variables"),attr(arandom,"variables"),attr(amixture,"variables"),attr(aclassmb,"variables"),attr(asurv,"variables"))
variables <- unlist(lapply(variables,all.vars))
if(!all(variables %in% colnames(data))) stop(paste("Data should contain the variables",paste(unique(variables),collapse=" ")))
###liste des variables utilisees (sans les interactions et sans Y)
if (ncor0>0) ttesLesVar <- unique(c(variables,cor.var.time))
else ttesLesVar <- unique(c(variables))
nomY <- as.character(attr(afixed,"variables")[2])
ttesLesVar <- setdiff(ttesLesVar,nomY)
###subset de data avec les variables utilisees
newdata <- data[,unique(c(nom.subject,nomY,noms.surv,ttesLesVar,nom.prior,nom.pprior,nom.timedepvar)),drop=FALSE]
## remplacer les NA de prior par 0
if(!is.null(nom.prior))
{
prior <- newdata[,nom.prior]
newdata[which(is.na(prior)),nom.prior] <- 0
prior[which(is.na(prior))] <- 0
}
##pprior : proba d'appartenance a chaque classe
if(!is.null(pprior))
{
pprior <- newdata[,nom.pprior]
}
## remplacer les NA de TimeDepVar par Tevent
Tint <- Tevent
nvdepsurv <- 0
if(!is.null(nom.timedepvar))
{
Tint <- newdata[,nom.timedepvar]
Tint[(is.na(Tint))] <- Tevent[(is.na(Tint))]
Tint[Tint>Tevent] <- Tevent[Tint>Tevent]
Tint[Tint<Tentry] <- Tentry[Tint<Tentry]
nvdepsurv <- 1
if (length(Tint[Tint<Tevent])==0)
{
stop("TimeDepVar is always greater than Time of Event. \n")
nvdepsurv <- 0
}
if (length(Tint[Tint>Tentry])==0)
{
Tint <- Tevent
stop("TimeDepVar is always lower than Time of Entry (0 by default). \n")
nvdepsurv <- 0
}
newdata[,nom.timedepvar] <- Tint
}
##enlever les NA
linesNA <- apply(newdata,2,function(v) which(is.na(v)))
linesNA <- unique(unlist(linesNA))
if(length(linesNA))
{
if(na.action==1) newdata <- newdata[-linesNA,,drop=FALSE]
if(na.action==2) stop("Data contain missing values.")
Tentry <- Tentry[-linesNA]
Tevent <- Tevent[-linesNA]
Event <- Event[-linesNA]
Tint <- Tint[-linesNA]
prior <- prior[-linesNA]
pprior <- pprior[-linesNA,]
}
### Y0
Y0 <- newdata[,nomY]
###creation de X0 (ttes les var + interactions)
Xintercept <- model.matrix(~1,data=newdata) #pr etre sur d'avoir I en premier
Xfixed <- model.matrix(fixed[-2], data=newdata)
Xmixture <- model.matrix(mixture, data=newdata)
Xrandom <- model.matrix(random, data=newdata)
Xclassmb <- model.matrix(classmb, data=newdata)
Xsurv <- model.matrix(form.commun,data=newdata)
Xsurvmix <- model.matrix(form.mixture,data=newdata)
Xsurvcause <- model.matrix(form.cause,data=newdata)
for (k in 1:nbevt)
{
assign(paste("Xsurvcause",k,sep=""),model.matrix(form.causek[[k]],data=newdata))
}
z.fixed <- strsplit(colnames(Xfixed),split=":",fixed=TRUE)
z.fixed <- lapply(z.fixed,sort)
if(random != ~-1)
{
z.random <- strsplit(colnames(Xrandom),split=":",fixed=TRUE)
z.random <- lapply(z.random,sort)
}
else
{
z.random <- list()
}
if(mixture != ~-1)
{
z.mixture <- strsplit(colnames(Xmixture),split=":",fixed=TRUE)
z.mixture <- lapply(z.mixture,sort)
}
else
{
z.mixture <- list()
}
if(classmb != ~-1)
{
z.classmb <- strsplit(colnames(Xclassmb),split=":",fixed=TRUE)
z.classmb <- lapply(z.classmb,sort)
}
else
{
z.classmb <- list()
}
if(form.commun != ~-1)
{
z.surv <- strsplit(colnames(Xsurv),split=":",fixed=TRUE)
z.surv <- lapply(z.surv,sort)
}
else
{
z.surv <- list()
}
if(form.mixture != ~-1)
{
z.survmix <- strsplit(colnames(Xsurvmix),split=":",fixed=TRUE)
z.survmix <- lapply(z.survmix,sort)
}
else
{
z.survmix <- list()
}
if(form.cause != ~-1)
{
z.survcause <- strsplit(colnames(Xsurvcause),split=":",fixed=TRUE)
z.survcause <- lapply(z.survcause,sort)
}
else
{
z.survcause <- list()
}
for(k in 1:nbevt)
{
if(form.causek[[k]] != ~-1)
{
assign(paste("z.survcause",k,sep=""),strsplit(colnames(get(paste("Xsurvcause",k,sep=""))),split=":",fixed=TRUE))
assign(paste("z.survcause",k,sep=""),lapply(get(paste("z.survcause",k,sep="")),sort))
}
else
{
assign(paste("z.survcause",k,sep=""),list())
}
}
if(!all(z.mixture %in% z.fixed)) stop("The covariates in mixture should also be included in the argument fixed")
X0 <- cbind(Xintercept,Xfixed, Xrandom, Xclassmb,Xsurv,Xsurvmix,Xsurvcause)
for(k in 1:nbevt)
{
X0 <- cbind(X0,get(paste("Xsurvcause",k,sep="")))
}
nom.unique <- unique(colnames(X0))
X0 <- X0[,nom.unique,drop=FALSE]
if(ncor0>0)
{
if(!(cor.var.time %in% colnames(X0)))
{
X0 <- cbind(X0, newdata[,cor.var.time])
colnames(X0) <- c(nom.unique, cor.var.time)
nom.unique <- colnames(X0)
}
}
X0 <- as.matrix(X0)
###X0 fini
timeobs <- rep(0, nrow(newdata))
if(!is.null(var.time))
{
timeobs <- newdata[,var.time]
if(any(is.na(timeobs))) stop(paste("Cannot use",var.time,"as time variable because it contains missing data"))
}
###test de link
if(link %in% c("splines","Splines"))
{
link <- "5-quant-splines"
}
idlink <- 2
if(link=="linear") idlink <- 0
if(link=="beta") idlink <- 1
if(link=="NULL") idlink <- -1
if (idlink==1)
{
if (epsY<=0)
{
stop("Argument 'epsY' should be positive.")
}
}
##remplir range si pas specifie
if(!is.null(range) & length(range)!=2) stop("Vector range should be of length 2.")
if(length(range)==2)
{
if(max(Y0)>range[2] | min(Y0)<range[1]) stop("The range specified do not cover the entire range of the data")
}
if(is.null(range))
{
min1 <- min(Y0)
min2 <- round(min1,3)
if(min1<min2) min2 <- min2-0.001
max1 <- max(Y0)
max2 <- round(max1,3)
if(max1>max2) max2 <- max2+0.001
range <- c(min2,max2)
}
nbzitr <- 2 #nbzitr = nb de noeuds si splines, 2 sinon
spltype <- NULL
if(idlink==2)
{
spl <- strsplit(link,"-")
if(length(spl[[1]])!=3) stop("Invalid argument link")
nbzitr <- as.numeric(spl[[1]][1])
spltype <- spl[[1]][2]
if(!(spltype %in% c("equi","quant","manual"))) stop("The location of the nodes should be 'equi', 'quant' or 'manual'")
if(!(spl[[1]][3] %in% c("splines","Splines"))) stop("Invalid argument link")
if(!is.null(intnodes))
{
if(spltype != "manual") stop("intnodes should be NULL if the nodes are not chosen manually")
if(length(intnodes) != (nbzitr-2)) stop(paste("Vector intnodes should be of length",nbzitr-2))
}
}
##remplir zitr (contient les noeuds si splines, sinon min et max)
zitr0 <- rep(0,nbzitr)
if(idlink %in% c(0,1)) zitr0[1:2] <- range
if(idlink==2)
{
if(spltype=="manual")
{
zitr0[1] <- range[1]
zitr0[nbzitr] <- range[2]
zitr0[2:(nbzitr-1)] <- intnodes
##verifier s'il y a des obs entre les noeuds
hcounts <- hist(Y0,breaks=zitr0,plot=FALSE,include.lowest=TRUE,right=TRUE)$counts
if(any(hcounts==0)) stop("Link function can not be estimated. Please try other nodes such that there are observations in each interval.")
}
if(spltype=="equi") zitr0[1:nbzitr] <- seq(range[1],range[2],length.out=nbzitr)
if(spltype=="quant")
{
Y0bis <- Y0
if(range[1]<min(Y0)) Y0bis <- c(range[1],Y0)
if(range[2]>max(Y0)) Y0bis <- c(range[2],Y0bis)
zitr0[1:nbzitr] <- quantile(Y0bis,probs=seq(0,1,length.out=nbzitr))
}
}
###uniqueY0 et indiceY0
uniqueY0 <- 0
indiceY0 <- 0
nvalSPL0 <- 0
if(idlink == 2)
{
uniqueY0 <- sort(unique(Y0))
permut <- order(order(Y0)) # sort(y)[order(order(y))] = y
indice <- rep(1:length(uniqueY0), as.vector(table(Y0)))
indiceY0 <- indice[permut]
nvalSPL0 <- length(uniqueY0)
}
###ordonner les mesures par individu
IND <- newdata[,nom.subject]
#IDnum <- as.numeric(IND)
if(is.null(nom.prior)) prior <- rep(0,length(Y0))
if(is.null(nom.pprior)) pprior <- matrix(1,length(Y0),ng)
if(!length(indiceY0)) indiceY0 <- rep(0,length(Y0))
matYX <- cbind(IND,timeobs,prior,pprior,Y0,indiceY0,Tentry,Tevent,Event,Tint,X0)
matYXord <- matYX[order(IND),]
Y0 <- as.numeric(matYXord[,4+ng])
X0 <- apply(matYXord[,-c(1:(9+ng)),drop=FALSE],2,as.numeric)
#IDnum <- matYXord[,1]
IND <- matYXord[,1]
indiceY0 <- as.numeric(matYXord[,5+ng])
prior0 <- as.numeric(unique(matYXord[,c(1,3)])[,2])
if(length(prior0)!=length(unique(IND))) stop("Please check 'prior' argument. Subjects can not have multiple assigned classes.")
pprior0 <- t(matYXord[!duplicated(matYXord[,1]), 3+1:ng, drop=FALSE])
timeobs <- matYXord[,2]
## nombre de sujets
ns0 <- length(unique(IND))
### Tevent, Tentry et Event de dim ns
#data.surv <- unique(cbind(IND,matYX[,c(6,7,8,9)]))
#if(nrow(data.surv) != ns0) stop("Subjects cannot have several times to event.")
nmes <- as.vector(table(IND))
##data.surv <- sweep(matYXord[cumsum(nmes),ng+c(6,7,8,9), drop=FALSE], MARGIN=1, STATS=0, FUN=function(x,y){as.numeric(x)+y})
data.surv <- matYXord[cumsum(nmes),ng+c(6,7,8,9), drop=FALSE]
tsurv0 <- as.numeric(data.surv[,1])
tsurv <- as.numeric(data.surv[,2])
devt <- as.numeric(data.surv[,3])
tsurvint <- as.numeric(data.surv[,4])
ind_survint <- (tsurvint<tsurv) + 0
### test de hazard
arghaz <- hazard
hazard <- rep(hazard,length.out=nbevt)
if(any(hazard %in% c("splines","Splines")))
{
hazard[which(hazard %in% c("splines","Splines"))] <- "5-quant-splines"
}
if(any(hazard %in% c("piecewise","Piecewise")))
{
hazard[which(hazard %in% c("piecewise","Piecewise"))] <- "5-quant-piecewise"
}
haz13 <- strsplit(hazard[which(!(hazard=="Weibull"))],"-")
if(any(sapply(haz13,length)!=3)) stop("Invalid argument hazard")
## si plusieurs splines, noeuds doivent etre identiques
nbspl <- length(which(sapply(haz13,getElement,3) %in% c("splines","Splines")))
if(nbspl>1)
{
haz3 <- haz13[which(sapply(haz13,getElement,3) %in% c("splines","Splines"))]
if(length(unique(sapply(haz3,getElement,1)))>1) stop("The nodes location of all splines hazard functions must be identical")
if(length(unique(sapply(haz3,getElement,2)))>1) stop("The nodes location of all splines hazard functions must be identical")
}
nz <- rep(2,nbevt)
locnodes <- NULL
typrisq <- rep(2,nbevt)
nprisq <- rep(2,nbevt)
nbnodes <- 0 #longueur de hazardnodes
nbrange <- 0 #longueur de hazardrange
ii <- 0
dejaspl <- 0
if(any(hazard!="Weibull"))
{
for (i in 1:nbevt)
{
if(hazard[i]=="Weibull") next;
ii <- ii+1
nz[i] <- as.numeric(haz13[[ii]][1])
if(nz[i]<3) stop("At least 3 nodes are required")
typrisq[i] <- ifelse(haz13[[ii]][3] %in% c("splines","Splines"),3,1)
nprisq[i] <- ifelse(haz13[[ii]][3] %in% c("splines","Splines"),nz[i]+2,nz[i]-1)
locnodes <- c(locnodes, haz13[[ii]][2])
if(!(haz13[[ii]][3] %in% c("splines","Splines","piecewise","Piecewise"))) stop("Invalid argument hazard")
if((haz13[[ii]][2]=="manual"))
{
if(typrisq[i]==1 | dejaspl==0)
{
if(length(arghaz)>1 | i==1 )
{
nbnodes <- nbnodes + nz[i]-2
nbrange <- nbrange + 2
}
}
if(typrisq[i]==3) dejaspl <- 1
}
if(!all(locnodes %in% c("equi","quant","manual"))) stop("The location of the nodes should be 'equi', 'quant' or 'manual'")
}
if(!is.null(hazardnodes))
{
if(!any(locnodes == "manual")) stop("hazardnodes should be NULL if the nodes are not chosen manually")
if(length(hazardnodes) != nbnodes) stop(paste("Vector hazardnodes should be of length",nbnodes))
}
}
else
{
if(!is.null(hazardnodes)) stop("hazardnodes should be NULL if Weibull baseline risk functions are chosen")
}
if(nbevt>1 & length(arghaz)==1 & nbnodes>0)
{
hazardnodes <- rep(hazardnodes,length.out=nbnodes*nbevt)
}
zi <- matrix(0,nrow=max(nz),ncol=nbevt)
nb <- 0
minT1 <- 0
maxT1 <- max(tsurv)
tsurvevt <- tsurv #[which(devt!=0)]
if(idtrunc==1)
{
minT1 <- min(tsurv,tsurv0)
maxT1 <- max(tsurv,tsurv0)
}
## arrondir
minT2 <- round(minT1,3)
if(minT1<minT2) minT2 <- minT2-0.001
minT <- minT2
maxT2 <- round(maxT1,3)
if(maxT1>maxT2) maxT2 <- maxT2+0.001
maxT <- maxT2
if(!is.null(hazardrange))
{
if(length(hazardrange) != nbrange) stop(paste("Vector hazardrange should be of length", nbrange))
if(any(hazardrange[seq(2,length(hazardrange), 2)] < maxT1) | any(hazardrange[seq(1,length(hazardrange), 2)] > minT1)) stop("The hazard range specified do not cover the entire range of the data")
if(nbevt>1 & length(arghaz)==1)
{
hazardrange <- rep(hazardrange,length.out=nbrange*nbevt)
}
}
ii <- 0
for(i in 1:nbevt)
{
if(typrisq[i]==2)
{
zi[1:2,i] <- c(minT,maxT)
}
else
{
ii <- ii+1
minTii <- minT
maxTii <- maxT
if(!is.null(hazardrange))
{
minTii <- hazardrange[2*(ii-1) + 1]
maxTii <- hazardrange[2*ii]
}
if(locnodes[ii]=="manual")
{
zi[1:nz[i],i] <- c(minTii,hazardnodes[nb+1:(nz[i]-2)],maxTii)
nb <- nb + nz[i]-2
}
if(locnodes[ii]=="equi")
{
zi[1:nz[i],i] <- seq(minTii,maxTii,length.out=nz[i])
}
if(locnodes[ii]=="quant")
{
#pi <- seq(0,1,length.out=nz[i])
#pi <- pi[-length(pi)]
#pi <- pi[-1]
pi <- c(1:(nz[i]-2))/(nz[i]-1)
qi <- quantile(tsurvevt,prob=pi)
zi[1,i] <- minTii
zi[2:(nz[i]-1),i] <- qi
zi[nz[i],i] <- maxTii
}
}
}
hazardtype <- rep(hazardtype,length.out=nbevt)
risqcom <- (hazardtype=="Common") + (hazardtype=="PH")*2
nrisq <- (risqcom==1)*nprisq + (risqcom==0)*nprisq*ng + (risqcom==2)*(nprisq+ng-1)
###parametres pour fortran
ng0 <- ng
nv0 <- dim(X0)[2]
nobs0 <- length(Y0)
idiag0 <- ifelse(idiag==TRUE,1,0)
nwg0 <- ifelse(nwg==TRUE,1,0)
nrisqtot <- sum(nrisq)
loglik <- 0
ni <- 0
istop <- 0
gconv <- rep(0,3)
ppi0 <- rep(0,ns0*ng0)
ppitest0 <- rep(0,ns0*ng0)
resid_m <- rep(0,nobs0)
resid_ss <- rep(0,nobs0)
pred_m_g <- rep(0,nobs0*ng0)
pred_ss_g <- rep(0,nobs0*ng0)
Yobs <- rep(0,nobs0)
#rlindiv <- rep(0,ns0)
marker <- rep(0,nsim)
transfY <- rep(0,nsim)
#Ydiscret <- 0
#UACV <- 0
#vraisdiscret <- 0
nmes0 <- as.vector(table(IND))
logspecif <- as.numeric(logscale)
time <- seq(minT,maxT,length.out=nsim)
risq_est <- matrix(0,nrow=nsim*ng0,ncol=nbevt)
risqcum_est <- matrix(0,nrow=nsim*ng0,ncol=nbevt)
statglob <- 0
statevt <- rep(0,nbevt)
###remplir idprob, etc
z.X0 <- strsplit(nom.unique,split=":",fixed=TRUE)
z.X0 <- lapply(z.X0,sort)
idprob0 <- z.X0 %in% z.classmb + 0
idea0 <- z.X0 %in% z.random + 0
idg0 <- (z.X0 %in% z.fixed) + (z.X0 %in% z.mixture)
## indicateurs pr variables survie
idcause <- rep(0,nv0)
for(k in 1:nbevt)
{
idcause <- idcause + (z.X0 %in% get(paste("z.survcause",k,sep="")) )
}
idcom <- rep(0,nv0)
idspecif <- matrix(0,ncol=nv0,nrow=nbevt)
for(j in 1:nv0)
{
if((z.X0[j] %in% z.surv) & !(z.X0[j] %in% z.survcause) & (idcause[j]==0))
{
idcom[j] <- 1
idspecif[,j] <- 1
#print("X")
}
if((z.X0[j] %in% z.survmix) & ( !(z.X0[j] %in% z.survcause) & (idcause[j]==0)))
{
idcom[j] <- 1
idspecif[,j] <- 2
#print("mixture(X)")
}
if((z.X0[j] %in% z.survmix) & (z.X0[j] %in% z.survcause))
{
idcom[j] <- 0
idspecif[,j] <- 2
#print("cause(mixture(X))")
}
if((z.X0[j] %in% z.survcause) & (!(z.X0[j] %in% z.survmix)))
{
idcom[j] <- 0
idspecif[,j] <- 1
#print("cause(X)")
}
if(idcause[j]!=0)
{
if(z.X0[j] %in% z.survmix)
{
for(k in 1:nbevt)
{
if(z.X0[j] %in% get(paste("z.survcause",k,sep="")))
{
idcom[j] <- 0
idspecif[k,j] <- 2
#cat("causek(mixture(X)) ,k=",k,"\n")
}
}
}
else
{
for(k in 1:nbevt)
{
if(z.X0[j] %in% get(paste("z.survcause",k,sep="")))
{
idcom[j] <- 0
idspecif[k,j] <- 1
#cat("causek(X) ,k=",k,"\n")
}
}
}
}
}
## mettre des 0 pour l'intercept
idcom[1] <- 0
idspecif[,1] <- 0
## quelle variable est TimeDepVar
idtdv <- z.X0 %in% nom.timedepvar + 0
#if(any(idcom>1) & nbevt<2) stop("No event specific effect can be estimated with less than two events")
if (ncor0>0) idcor0 <- z.X0 %in% cor.var.time +0
else idcor0 <- rep(0,nv0)
## Si pas TimeDepVar dans formule survival
if(length(nom.timedepvar) & all(idtdv==0))
{
stop("Variable in 'TimeDepVar' should also appear as a covariate in the 'survival' argument")
## ou l'ajouter?
}
## si on a ignore TimeDepVar #pas utile si stop
if(length(nom.timedepvar) & nvdepsurv==0)
{
jj <- which(idtdv==1)
idtdv[jj] <- 0
idcom[jj] <- 0
idspecif[,jj] <- 0
idcause[jj] <- 0
nom.timedepvar <- NULL
}
nea0 <- sum(idea0)
predRE <- rep(0,ns0*nea0)
## nb coef pr survie
nvarxevt <- 0
nvarxevt2 <- 0 # pr valeurs initiales calculees avec Fortran
for(j in 1:nv0)
{
if(idcom[j]==1)
{
if(all(idspecif[,j]==1))
{
nvarxevt <- nvarxevt + 1
nvarxevt2 <- nvarxevt2 + 1
}
if(all(idspecif[,j]==2))
{
nvarxevt <- nvarxevt + ng
nvarxevt2 <- nvarxevt2 + 1
}
}
if(idcom[j]==0)
{
if(all(idspecif[,j]==0)) next
for(k in 1:nbevt)
{
if(idspecif[k,j]==1)
{
nvarxevt <- nvarxevt + 1
nvarxevt2 <- nvarxevt2 + 1
}
if(idspecif[k,j]==2)
{
nvarxevt <- nvarxevt + ng
nvarxevt2 <- nvarxevt2 + 1
}
}
}
}
ntrtot0 <- (idlink==-1) + 2*(idlink==0) + 4*(idlink==1) + (nbzitr+2)*(idlink==2)
nef <- sum(idg0==1) + ng0*sum(idg0==2)
if(idlink!=-1) nef <- nef-1
nvc <- ifelse(idiag0==1,nea0,nea0*(nea0+1)/2)
nw <- (ng0-1)*nwg0
nprob <- sum(idprob0)*(ng0-1)
##nombre total de parametres
NPM <- (ng0-1)*sum(idprob0) +
nrisqtot + nvarxevt +
nef + nvc + nw + ncor0 + ntrtot0
V <- rep(0, NPM*(NPM+1)/2) #pr variance des parametres
## prm fixes
fix0 <- rep(0,NPM)
if(length(posfix))
{
if(any(!(posfix %in% 1:NPM))) stop("Indexes in posfix are not correct")
fix0[posfix] <- 1
}
if(length(posfix)==NPM) stop("No parameters to estimate")
## pour H restreint
Hr0 <- as.numeric(partialH)
pbH0 <- rep(0,NPM)
if(is.logical(partialH))
{
if(partialH)
{
if(any(typrisq %in% c(1,3)))
{
for(k in 1:nbevt)
{
if(typrisq[k] %in% c(1,3))
{
pbH0[nprob+sum(nrisq[1:k])-nrisq[k]+1:nrisq[k]] <- 1
if(risqcom[k]==2)
{
pbH0[nprob+sum(nrisq[1:k])] <- 0
}
}
}
}
if(any(idlink==2))
{
pbH0[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+1:ntrtot0] <- 1
}
}
#pbH0[posfix] <- 0
if(sum(pbH0)==0 & Hr0==1) stop("No partial Hessian matrix can be defined in the absence of baseline risk function approximated by splines or piecewise linear function")
}
else
{
if(!all(Hr0 %in% 1:NPM)) stop("Indexes in partialH are not correct")
pbH0[Hr0] <- 1
#pbH0[posfix] <- 0
}
indexHr <- NULL
if(sum(pbH0)>0)
{
if(length(posfix)) pbH1 <- pbH0[-posfix] else pbH1 <- pbH0
indexHr <- which(pbH1==1)
}
## gestion de B=random(mod)
Brandom <- FALSE
tryB <- try(as.numeric(B), silent=TRUE)
if(inherits(tryB, "try-error"))
{
if(length(cl$B)==2)
{
if(!inherits(eval(cl$B[[2]], parent.env(environment())),"Jointlcmm")) stop("The model specified in B should be of class Jointlcmm")
if(as.character(cl$B[1])!="random") stop("Please use random() to specify random initial values")
Brandom <- TRUE
B <- eval(cl$B[[2]], parent.env(environment()))
if(B$conv != 1) stop("Model in argument B did not converge properly")
#if(length(posfix)) stop("Argument posfix is not compatible with random intial values")
}
}
##valeurs initiales
if(!(missing(B)))
{
if(is.vector(B))
{
if (length(B)==NPM)
{
b <- B
if(nvc>0)
{
## remplacer varcov des EA par les prm a estimer
if(idiag0==1)
{
b[nprob+nrisqtot+nvarxevt+nef+1:nvc] <- sqrt(b[nprob+nrisqtot+nvarxevt+nef+1:nvc])
}
else
{
varcov <- matrix(0,nrow=nea0,ncol=nea0)
varcov[upper.tri(varcov,diag=TRUE)] <- b[nprob+nrisqtot+nvarxevt+nef+1:nvc]
varcov <- t(varcov)
varcov[upper.tri(varcov,diag=TRUE)] <- b[nprob+nrisqtot+nvarxevt+nef+1:nvc]
ch <- chol(varcov)
b[nprob+nrisqtot+nvarxevt+nef+1:nvc] <- ch[upper.tri(ch,diag=TRUE)]
}
}
}
else
{
stop(paste("Vector B should be of length",NPM))
}
}
else
{
if(!inherits(B,"Jointlcmm")) stop("B should be either a vector or an object of class Jointlcmm")
b <- rep(0,NPM)
if(ng>1 & B$ng==1)
{
nef2 <- sum(idg0!=0)
if(idlink != -1)
{
nef2 <- nef2-1 #car intercept fixe a 0
}
NPM2 <- sum(nprisq)+nvarxevt2+nef2+nvc+ncor0+ntrtot0
if(length(B$best)!=NPM2) stop("B is not correct")
if(Brandom==FALSE)
{
## B deterministe
for(ke in 1:nbevt)
{
if(risqcom[ke]==0)
{
indb <- nprob+sum(nrisq[1:ke])-nrisq[ke]+1:nrisq[ke]
indinit <- sum(nprisq[1:ke])-nprisq[ke]+1:nprisq[ke]
if(B$conv==1)
{
b[indb] <- rep(B$best[indinit],ng0)+rep((1:ng0)-(ng0+1)/2,each=nprisq[ke])*rep(sqrt(B$V[indinit*(indinit+1)/2]),ng0)
}
else
{
b[indb] <- rep(B$best[indinit],ng0)+rep((1:ng0)-(ng0+1)/2,each=nprisq[ke])*rep(B$best[indinit],ng0)
}
}
if(risqcom[ke]==1)
{
b[nprob+sum(nrisq[1:ke])-nrisq[ke]+1:nrisq[ke]] <- B$best[sum(nprisq[1:ke])-nprisq[ke]+1:nprisq[ke]]
}
if(risqcom[ke]==2)
{
b[nprob+sum(nrisq[1:ke])-nrisq[ke]+1:nrisq[ke]] <- c(B$best[sum(nprisq[1:ke])-nprisq[ke]+1:nprisq[ke]],0.5+(0:(ng0-2))*0.5)
}
}
## pr bevt
avtj <- nprob+nrisqtot
avtj2 <- sum(nprisq)
for(j in 1:nv0)
{
if(idcom[j]==0 & all(idspecif[,j]==0)) next
if(idcom[j]==1 & all(idspecif[,j]==1))
{
b[avtj+1] <- B$best[avtj2+1]
avtj <- avtj+1
avtj2 <- avtj2+1
}
if(idcom[j]==1 & all(idspecif[,j]==2))
{
if(B$conv==1) b[avtj+1:ng0] <- abs(rep(B$best[avtj2+1],ng0)+c(1:ng0-(ng0+1)/2)*rep(sqrt(B$V[(avtj2+1)*(avtj2+2)/2]),ng0))
else b[avtj+1:ng0] <- abs(rep(B$best[avtj2+1],ng0)+c(1:ng0-(ng0+1)/2)*rep(B$best[avtj2+1],ng0))
avtj <- avtj+ng0
avtj2 <- avtj2+1
}
if(idcom[j]==0 & any(idspecif[,j]!=0))
{
for(k in 1:nbevt)
{
if(idspecif[k,j]==0) next
if(idspecif[k,j]==1)
{
b[avtj+1] <- B$best[avtj2+1]
avtj <- avtj+1
avtj2 <- avtj2+1
}
if(idspecif[k,j]==2)
{
if(B$conv==1) b[avtj+1:ng0] <- abs(rep(B$best[avtj2+1],ng0)+c(1:ng0-(ng0+1)/2)*rep(sqrt(B$V[(avtj2+1)*(avtj2+2)/2]),ng0))
else b[avtj+1:ng0] <- abs(rep(B$best[avtj2+1],ng0)+c(1:ng0-(ng0+1)/2)*rep(B$best[avtj2+1],ng0))
avtj <- avtj+ng0
avtj2 <- avtj2+1
}
}
}
}
## pr nef
avtj <- nprob+nrisqtot+nvarxevt
avtj2 <- sum(nprisq)+nvarxevt2
for(j in 1:nv0)
{
if(idg0[j]==1)
{
if(j==1 & idlink!=-1) next
else
{
if(B$conv==1) b[avtj+1] <- B$best[avtj2+1]
avtj <- avtj+1
avtj2 <- avtj2+1
}
}
if(idg0[j]==2)
{
if(j==1)
{
if(idlink!=-1)
{
b[avtj+1:(ng0-1)] <- -0.5*(1:(ng0-1))
avtj <- avtj+ng0-1
}
else
{
b[avtj+1:ng0] <- B$best[avtj2+1]+(1:ng0-(ng0+1)/2)*B$best[avtj2+1]
avtj <- avtj+ng0
avtj2 <- avtj2+1
}
}
else
{
if(B$conv==1) b[avtj+1:ng0] <- B$best[avtj2+1]+(1:ng0-(ng0+1)/2)*sqrt(B$V[(avtj2+1)*(avtj2+1+1)/2])
else b[avtj+1:ng0] <- B$best[avtj2+1]+(1:ng0-(ng0+1)/2)*B$best[avtj2+1]
avtj <- avtj+ng0
avtj2 <- avtj2+1
}
}
}
## pr nvc
if(nvc>0)
{
b[nprob+nrisqtot+nvarxevt+nef+1:nvc] <- B$cholesky
}
## cor et transfo
b[(nprob+nrisqtot+nvarxevt+nef+nvc+nw+1):NPM] <- B$best[(sum(nprisq)+nvarxevt2+nef2+nvc+1):NPM2]
}
else
{
## B random
bb <- rep(0,NPM-nprob-(ng-1)*sum(length(which(risqcom==2)))-nw) # tous les prm sauf nprob, bPH et nw
vbb <- matrix(0,length(bb),length(bb))
VB <- matrix(0,NPM2,NPM2)
VB[upper.tri(VB,diag=TRUE)] <- B$V
VB <- t(VB)
VB[upper.tri(VB,diag=TRUE)] <- B$V
copcov <- rep(0,length(bb))
copcov2 <- rep(0,length(B$best))
avt <- 0
for(ke in 1:nbevt)
{
if(risqcom[ke]==0)
{
indbb <- avt+1:nrisq[ke]
indB <- sum(nprisq[1:ke])-nprisq[ke]+1:nprisq[ke]
bb[indbb] <- rep(B$best[indB],ng)
diag(vbb[indbb,indbb]) <- rep(diag(VB[indB,indB]),ng)
avt <- avt + nrisq[ke]
}
if(risqcom[ke]==1)
{
indbb <- avt+1:nrisq[ke]
indB <- sum(nprisq[1:ke])-nprisq[ke]+1:nprisq[ke]
bb[indbb] <- B$best[indB]
copcov[indbb] <- 1
copcov2[indB] <- 1
avt <- avt + nrisq[ke]
}
if(risqcom[ke]==2)
{
indbb <- avt+1:nprisq[ke]
indB <- sum(nprisq[1:ke])-nprisq[ke]+1:nprisq[ke]
bb[indbb] <- B$best[indB]
copcov[indbb] <- 1
copcov2[indB] <- 1
avt <- avt + nprisq[ke]
}
}
avt2 <- sum(nprisq)
for(j in 1:nv0)
{
if(idcom[j]==0 & all(idspecif[,j]==0)) next
if(idcom[j]==1 & all(idspecif[,j]==1))
{
bb[avt+1] <- B$best[avt2+1]
copcov[avt+1] <- 1
copcov2[avt2+1] <- 1
avt <- avt+1
avt2 <- avt2+1
}
if(idcom[j]==1 & all(idspecif[,j]==2))
{
bb[avt+1:ng] <- rep(B$best[avt2+1],ng)
diag(vbb[avt+1:ng,avt+1:ng]) <- rep(VB[avt2+1,avt2+1],ng)
avt <- avt+ng
avt2 <- avt2+1
}
if(idcom[j]==0 & any(idspecif[,j]!=0))
{
for(k in 1:nbevt)
{
if(idspecif[k,j]==0) next
if(idspecif[k,j]==1)
{
bb[avt+1] <- B$best[avt2+1]
copcov[avt+1] <- 1
copcov2[avt2+1] <- 1
avt <- avt+1
avt2 <- avt2+1
}
if(idspecif[k,j]==2)
{
bb[avt+1:ng] <- rep(B$best[avt2+1],ng)
diag(vbb[avt+1:ng,avt+1:ng]) <- rep(VB[avt2+1,avt2+1],ng)
avt <- avt+ng
avt2 <- avt2+1
}
}
}
}
for(j in 1:nv0)
{
if(idg0[j]==1)
{
if(j==1 & idlink!=-1) next
else
{
bb[avt+1] <- B$best[avt2+1]
copcov[avt+1] <- 1
copcov2[avt2+1] <- 1
avt <- avt+1
avt2 <- avt2+1
}
}
if(idg0[j]==2)
{
if(j==1 & idlink!=-1)
{
avt <- avt+ng-1
}
else
{
bb[avt+1:ng] <- rep(B$best[avt2+1],ng)
diag(vbb[avt+1:ng,avt+1:ng]) <- rep(VB[avt2+1,avt2+1],ng)
avt <- avt+ng
avt2 <- avt2+1
}
}
}
if(nvc>0)
{
bb[avt+1:nvc] <- B$cholesky
copcov[avt+1:nvc] <- 1
copcov2[avt2+1:nvc] <- 1
avt <- avt+nvc
avt2 <- avt2+nvc
}
bb[(avt+1):length(bb)] <- B$best[(avt2+1):NPM2]
copcov[(avt+1):length(bb)] <- 1
copcov2[(avt2+1):length(B$best)] <- 1
vbb[which(copcov==1),which(copcov==1)] <- VB[which(copcov2==1),which(copcov2==1)]
if(idlink!=-1 & idg0[1]>1)
{
bb <- bb[-(sum(nrisq)-(ng-1)*sum(length(which(risqcom==2)))+nvarxevt+1:(ng-1))]
vbb <- vbb[-(sum(nrisq)-(ng-1)*sum(length(which(risqcom==2)))+nvarxevt+1:(ng-1)),-(sum(nrisq)-(ng-1)*sum(length(which(risqcom==2)))+nvarxevt+1:(ng-1))]
}
##Chol <- chol(vbb)
##Chol <- t(Chol)
bb <- rmvnorm(n=1,mean=bb,sigma = vbb) #bb + Chol %*% rnorm(length(bb))
if(nprob>0) b[1:nprob] <- 0
avt <- 0
for(ke in 1:nbevt)
{
if(risqcom[ke]==0)
{
b[nprob+sum(nrisq[1:ke])-nrisq[ke]+1:nrisq[ke]] <- bb[avt+1:nrisq[ke]]
avt <- avt + nrisq[ke]
}
if(risqcom[ke]==1)
{
b[nprob+sum(nrisq[1:ke])-nrisq[ke]+1:nrisq[ke]] <- bb[avt+1:nrisq[ke]]
avt <- avt + nrisq[ke]
}
if(risqcom[ke]==2)
{
b[nprob+sum(nrisq[1:ke])-nrisq[ke]+1:nprisq[ke]] <- bb[avt+1:nprisq[ke]]
b[nprob+sum(nrisq[1:ke])-nrisq[ke]+nprisq[ke]+1:(ng-1)] <- 1 #bPH
avt <- avt + nprisq[ke]
}
}
if(nvarxevt>0) b[nprob+nrisqtot+1:nvarxevt] <- bb[avt+1:nvarxevt]
if(idlink!=-1 & idg0[1]>1)
{
b[nprob+nrisqtot+nvarxevt+1:(ng-1)] <- 0
b[(nprob+nrisqtot+nvarxevt+ng):(nprob+nrisqtot+nvarxevt+nef)] <- bb[avt+nvarxevt+1:(nef-(ng-1))]
nefssI <- nef-ng+1
}
else
{
b[nprob+nrisqtot+nvarxevt+1:nef] <- bb[avt+nvarxevt+1:nef]
nefssI <- nef
}
if(nvc>0)
{
b[nprob+nrisqtot+nvarxevt+nef+1:nvc] <- bb[nrisqtot-((ng-1)*length(which(risqcom==2)))+nvarxevt+nefssI+1:nvc]
}
if(nw>0) b[nprob+nrisqtot+nvarxevt+nef+nvc+1:nw] <- 1
b[(nprob+nrisqtot+nvarxevt+nef+nvc+nw+1):NPM] <- bb[(nrisqtot-((ng-1)*length(which(risqcom==2)))+nvarxevt+nefssI+nvc+1):length(bb)]
}
}
}
}
else
{
if(ng>1) stop("Please specify initial values with argument 'B'")
b <- rep(0,NPM)
for(i in 1:nbevt)
{
if (typrisq[i]==2)
{
if(logspecif==1)
{
b[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- c(rep(c(log(sum(devt==i)/sum(tsurv[devt==i])),0),ifelse(risqcom[i]==0,ng,1)),rep(1,(ng-1)*(risqcom[i]==2)))
}
else
{
b[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- c(rep(c(sqrt(sum(devt==i)/sum(tsurv[devt==i])),1),ifelse(risqcom[i]==0,ng,1)),rep(1,(ng-1)*(risqcom[i]==2)))
}
}
else
{
if(logspecif==1)
{
b[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- c(rep(log(1/nprisq[i]),ifelse(risqcom[i]==0,ng*nprisq[i],nprisq[i])),rep(1,(ng-1)*(risqcom[i]==2)))
}
else
{
b[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- c(rep(sqrt(1/nprisq[i]),ifelse(risqcom[i]==0,ng*nprisq[i],nprisq[i])),rep(1,(ng-1)*(risqcom[i]==2)))
}
}
}
if (nvc>0)
{
if(idiag0==1) b[nprob+nrisqtot+nvarxevt+nef+1:nvc] <- rep(1,nvc)
if(idiag0==0)
{
init.nvc <- diag(nea0)
init.nvc <- init.nvc[upper.tri(init.nvc, diag=TRUE)]
b[nprob+nrisqtot+nvarxevt+nef+1:nvc] <- init.nvc
}
}
if(nwg0>0) b[nprob+nrisqtot+nvarxevt+nef+nvc+1:nw] <- 1
if(ncor0==1) b[nprob+nrisqtot+nvarxevt+nef+nvc+nw+1] <- 1
if(ncor0==2) b[nprob+nrisqtot+nvarxevt+nef+nvc+nw+1:2] <- c(0,1)
if(idlink==-1) b[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+1] <- 1
if(idlink==0)
{
b[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+1] <- mean(Y0)
b[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+2] <- 1
}
if(idlink==1)
{
b[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+1] <- 0
b[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+2] <- -log(2)
b[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+3] <- 0.7
b[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+4] <- 0.1
}
if(idlink==2)
{
b[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+1] <- -2
b[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+2:ntrtot0] <- 0.1
}
#cat("B inti pour ng=1 ", b,"\n")
}
#cat("B init a partir du modele pour ng=1 ", b,"\n")
### nom au vecteur best
nom.X0 <- colnames(X0)
nom.X0[which(nom.X0=="(Intercept)")] <- "intercept"
##prm classmb
if(ng0>=2 & nprob>0)
{
nom <- rep(nom.X0[idprob0==1],each=ng0-1)
nom1 <- paste(nom," class",c(1:(ng0-1)),sep="")
names(b)[1:nprob] <- nom1
}
##prm fct de risque
if(isTRUE(logscale))
{
for(i in 1:nbevt)
{
nom1 <- rep(paste("event",i,sep=""),nrisq[i])
if(typrisq[i]==2)
{
nom2 <- paste(nom1[1:2]," log(Weibull",1:2,")",sep="")
nom1[1:(2*ifelse(risqcom[i]==0,ng,1))] <- rep(nom2,ifelse(risqcom[i]==0,ng*(risqcom[i]==0),1))
if(risqcom[i]==2) nom1[2+1:(ng-1)] <- paste(nom1[2+1:(ng-1)],"SurvPH")
names(b)[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- nom1
}
if(typrisq[i]==1)
{
nom2 <- paste(nom1[1:(nz[i]-1)]," log(piecewise",1:(nz[i]-1),")",sep="")
nom1[1:((nz[i]-1)*ifelse(risqcom[i]==0,ng,1))] <- rep(nom2,ifelse(risqcom[i]==0,ng*(risqcom[i]==0),1))
if(risqcom[i]==2) nom1[nz[i]-1+1:(ng-1)] <- paste(nom1[nz[i]-1+1:(ng-1)],"SurvPH")
names(b)[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- nom1
}
if(typrisq[i]==3)
{
nom2 <- paste(nom1[1:(nz[i]-1)]," log(splines",1:(nz[i]+2),")",sep="")
nom1[1:((nz[i]+2)*ifelse(risqcom[i]==0,ng,1))] <- rep(nom2,ifelse(risqcom[i]==0,ng*(risqcom[i]==0),1))
if(risqcom[i]==2) nom1[nz[i]+2+1:(ng-1)] <- paste(nom1[nz[i]+2+1:(ng-1)],"SurvPH")
names(b)[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- nom1
}
}
}
else
{
for(i in 1:nbevt)
{
nom1 <- rep(paste("event",i,sep=""),nrisq[i])
if(typrisq[i]==2)
{
nom2 <- paste(nom1[1:2]," +/-sqrt(Weibull",1:2,")",sep="")
nom1[1:(2*ifelse(risqcom[i]==0,ng,1))] <- rep(nom2,ifelse(risqcom[i]==0,ng*(risqcom[i]==0),1))
if(risqcom[i]==2) nom1[2+1:(ng-1)] <- paste(nom1[2+1:(ng-1)],"SurvPH")
names(b)[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- nom1
}
if(typrisq[i]==1)
{
nom2 <- paste(nom1[1:(nz[i]-1)]," +/-sqrt(piecewise",1:(nz[i]-1),")",sep="")
nom1[1:((nz[i]-1)*ifelse(risqcom[i]==0,ng,1))] <- rep(nom2,ifelse(risqcom[i]==0,ng*(risqcom[i]==0),1))
if(risqcom[i]==2) nom1[nz[i]-1+1:(ng-1)] <- paste(nom1[nz[i]-1+1:(ng-1)],"SurvPH")
names(b)[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- nom1
}
if(typrisq[i]==3)
{
nom2 <- paste(nom1[1:(nz[i]-1)]," +/-sqrt(splines",1:(nz[i]+2),")",sep="")
nom1[1:((nz[i]+2)*ifelse(risqcom[i]==0,ng,1))] <- rep(nom2,ifelse(risqcom[i]==0,ng*(risqcom[i]==0),1))
if(risqcom[i]==2) nom1[nz[i]+2+1:(ng-1)] <- paste(nom1[nz[i]+2+1:(ng-1)],"SurvPH")
names(b)[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- nom1
}
}
}
if(ng>1)
{
for(i in 1:nbevt)
{
if(risqcom[i]==1) next;
if(risqcom[i]==0)
{
names(b)[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- paste(names(b)[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]],paste("class",rep(1:ng,each=nprisq[i])))
}
if(risqcom[i]==2)
{
names(b)[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]] <- paste(names(b)[nprob+sum(nrisq[1:i])-nrisq[i]+1:nrisq[i]],c(rep("",nprisq[i]),paste(" class",1:(ng-1),sep="")),sep="")
}
}
}
##prm covariables survival
nom1 <- NULL
for(j in 1:nv0)
{
if(idcom[j]==0 & all(idspecif[,j]==0)) next
if(idcom[j]==1 & all(idspecif[,j]==1)) #X
{
if(idtdv[j]==1)
{
nom1 <- c(nom1,paste("I(t>",nom.timedepvar,")",sep=""))
}
else
{
nom1 <- c(nom1,nom.X0[j])
}
next
}
if(idcom[j]==1 & all(idspecif[,j]==2)) #mixture(X)
{
if(idtdv[j]==1)
{
nom1 <- c(nom1,paste("I(t>",nom.timedepvar,") class",1:ng0,sep=""))
}
else
{
nom1 <- c(nom1,paste(nom.X0[j],paste("class",1:ng0,sep="")))
}
next
}
if(idcom[j]==0 & all(idspecif[,j]==1)) #cause(X)
{
if(idtdv[j]==1)
{
nom1 <- c(nom1,paste("I(t>",nom.timedepvar,") event",1:nbevt,sep=""))
}
else
{
nom1 <- c(nom1,paste(nom.X0[j],paste("event",1:nbevt,sep="")))
}
next
}
if(idcom[j]==0 & all(idspecif[,j]==2)) #cause(mixture(X))
{
if(idtdv[j]==1)
{
xevt <- paste("I(t>",nom.timedepvar,") event",1:nbevt,sep="")
classg <- paste("class",1:ng0,sep="")
nom1 <- c(nom1,paste(rep(xevt,each=ng0),rep(classg,nbevt)))
}
else
{
xevt <- paste(nom.X0[j],paste("event",1:nbevt,sep=""))
classg <- paste("class",1:ng0,sep="")
nom1 <- c(nom1,paste(rep(xevt,each=ng0),rep(classg,nbevt)))
}
next
}
if(idcom[j]==0 & idcause[j]!=0) #causek
{
for(k in 1:nbevt)
{
if(idspecif[k,j]==0) next
if(idspecif[k,j]==1)
{
if(idtdv[j]==1)
{
nom1 <- c(nom1,paste("I(t>",nom.timedepvar,") event",k,sep=""))
}
else
{
nom1 <- c(nom1,paste(nom.X0[j],paste("event",k,sep="")))
}
next
}
if(idspecif[k,j]==2)
{
if(idtdv[j]==1)
{
xevtk <- paste("I(t>",nom.timedepvar,") event",k,sep="")
classg <- paste("class",1:ng0,sep="")
nom1 <- c(nom1,paste(xevtk,classg))
}
else
{
xevtk <- paste(nom.X0[j],paste("event",k,sep=""))
classg <- paste("class",1:ng0,sep="")
nom1 <- c(nom1,paste(xevtk,classg))
}
next
}
}
}
}
if(nvarxevt>0) names(b)[nprob+nrisqtot+1:nvarxevt] <- nom1
##NB : pour chaque variable, coef ranges par evenement puis par classe
##prm fixed
if(ng0==1 & nef>0)
{
if(idlink!=-1)
{
names(b)[nprob+nrisqtot+nvarxevt+1:nef] <- nom.X0[-1][idg0[-1]!=0]
}
else
{
# print(nprob+nrisqtot+nvarxevt)
# print(nef)
# print(nom.X0)
# print(idg0)
names(b)[nprob+nrisqtot+nvarxevt+1:nef] <- nom.X0[idg0!=0]
}
}
if(ng0>1)
{
nom1<- NULL
for(i in 1:nv0)
{
if(idg0[i]==2)
{
if(i==1)
{
if(idlink==-1)
{
nom <- paste("intercept class",c(1:ng0),sep="")
nom1 <- c(nom1,nom)
}
else
{
nom <- paste("intercept class",c(2:ng0),sep="")
nom1 <- c(nom1,nom)
}
}
if (i>1)
{
nom <- paste(nom.X0[i]," class",c(1:ng0),sep="")
nom1 <- c(nom1,nom)
}
}
if(idg0[i]==1 & i==1 & idlink==-1) nom1 <- c(nom1,nom.X0[i])
if(idg0[i]==1 & i>1) nom1 <- c(nom1,nom.X0[i])
}
names(b)[nprob+nrisqtot+nvarxevt+1:nef]<- nom1
}
##prm random et nwg
if(nvc!=0) names(b)[nprob+nrisqtot+nvarxevt+nef+1:nvc] <- paste("varcov",c(1:nvc))
if(nw!=0) names(b)[nprob+nrisqtot+nvarxevt+nef+nvc+1:nw] <- paste("varprop class",c(1:(ng0-1)))
##prm cor
if(ncor0>0) {names(b)[nprob+nrisqtot+nvarxevt+nef+nvc+nw+1:ncor0] <- paste("cor",1:ncor0,sep="")}
##prm link
if(idlink==-1) names(b)[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+1] <- "stderr"
if(idlink==0) names(b)[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+1:2]<- c("Linear 1","Linear 2")
if(idlink==1) names(b)[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+1:4]<- paste("Beta",1:4,sep="")
if(idlink==2) names(b)[nprob+nrisqtot+nvarxevt+nef+nvc+nw+ncor0+1:ntrtot0]<- paste("I-splines",c(1:ntrtot0),sep="")
names_best <- names(b)
#browser()
#print(b)
#### estimation
idspecif <- as.vector(t(idspecif))
#return(b)
nfix <- sum(fix0)
bfix <- 0
if(nfix>0)
{
bfix <- b[which(fix0==1)]
b <- b[which(fix0==0)]
NPM <- NPM-nfix
}
if(maxiter==0)
{
vrais <- loglikJointlcmm(b,Y0,X0,prior0,pprior0,tsurv0,tsurv,devt,ind_survint,
idprob0,idea0,idg0,idcor0,idcom,idspecif,idtdv,
idlink,epsY,nbzitr,zitr0,uniqueY0,nvalSPL0,
indiceY0,typrisq,risqcom,nz,zi,ns0,ng0,nv0,
nobs0,nmes0,nbevt,nea0,nw,ncor0,idiag0,idtrunc,
logspecif,NPM,fix0,nfix,bfix)
out <- list(conv=2, V=rep(NA, length(b)), best=b,
ppi=rep(NA,ns0*ng0), ppitest=rep(NA,ns0*ng0),
predRE=rep(NA,ns0*nea0), Yobs=rep(NA,nobs0),
resid_m=rep(NA,nobs0), resid_ss=rep(NA,nobs0),
marker=rep(NA,nsim), transfY=rep(NA,nsim),
pred_m_g=rep(NA,nobs0*ng0), pred_ss_g=rep(NA,nobs0*ng0),
time=rep(NA,nsim), risq_est=rep(NA,nsim*ng0*nbevt),
risqcum_est=rep(NA,nsim*ng0*nbevt),
statglob=NA, statevt=rep(NA,nbevt),
gconv=rep(NA,3), niter=0, loglik=vrais)
}
else
{
res <- mla(b=b, m=length(b), fn=loglikJointlcmm,
clustertype=clustertype,.packages=NULL,
epsa=convB,epsb=convL,epsd=convG,partialH=indexHr,
digits=8,print.info=verbose,blinding=FALSE,
multipleTry=25,file="",
nproc=nproc,maxiter=maxiter,minimize=FALSE,
Y0=Y0,X0=X0,prior0=prior0,pprior0=pprior0,tentr0=tsurv0,tevt0=tsurv,devt0=devt,
ind_survint0=ind_survint,idprob0=idprob0,idea0=idea0,idg0=idg0,
idcor0=idcor0,idcom0=idcom,idspecif0=idspecif,idtdv0=idtdv,
idlink0=idlink,epsY0=epsY,nbzitr0=nbzitr,zitr0=zitr0,
uniqueY0=uniqueY0,nvalSPL0=nvalSPL0,indiceY0=indiceY0,
typrisq0=typrisq,risqcom0=risqcom,nz0=nz,zi0=zi,ns0=ns0,
ng0=ng0,nv0=nv0,nobs0=nobs0,nmes0=nmes0,nbevt0=nbevt,nea0=nea0,
nwg0=nw,ncor0=ncor0,idiag0=idiag0,idtrunc0=idtrunc,
logspecif0=logspecif,npm0=NPM,fix0=fix0,nfix0=nfix,bfix0=bfix)
out <- list(conv=res$istop, V=res$v, best=res$b,
ppi=rep(NA,ns0*ng0), ppitest=rep(NA,ns0*ng0),
predRE=rep(NA,ns0*nea0), Yobs=rep(NA,nobs0),
resid_m=rep(NA,nobs0), resid_ss=rep(NA,nobs0),
marker=rep(NA,nsim), transfY=rep(NA,nsim),
pred_m_g=rep(NA,nobs0*ng0), pred_ss_g=rep(NA,nobs0*ng0),
time=rep(NA,nsim), risq_est=rep(NA,nsim*ng0*nbevt),
risqcum_est=rep(NA,nsim*ng0*nbevt),
statglob=NA, statevt=rep(NA,nbevt),
gconv=c(res$ca, res$cb, res$rdm), niter=res$ni,
loglik=res$fn.value)
}
if(out$conv %in% c(1,2,3))
{
estim0 <- 0
ll <- 0
ppi0 <- rep(0,ns0*ng0)
ppitest0 <- rep(0,ns0*ng0)
resid_m <- rep(0,nobs0)
resid_ss <- rep(0,nobs0)
pred_m_g <- rep(0,nobs0*ng0)
pred_ss_g <- rep(0,nobs0*ng0)
predRE <- rep(0,ns0*nea0)
marker <- rep(0,nsim)
transfY <- rep(0,nsim)
Yobs <- rep(0,nobs0)
risq_est <- rep(0,nsim*ng0*nbevt)
risqcum_est <- rep(0,nsim*ng0*nbevt)
statglob <- 0
statevt <- rep(0,nbevt)
estim0 <- 0
post <- .Fortran(C_loglikjointlcmm,
as.double(Y0),
as.double(X0),
as.integer(prior0),
as.double(pprior0),
as.double(tsurv0),
as.double(tsurv),
as.integer(devt),
as.integer(ind_survint),
as.integer(idprob0),
as.integer(idea0),
as.integer(idg0),
as.integer(idcor0),
as.integer(idcom),
as.integer(idspecif),
as.integer(idtdv),
as.integer(idlink),
as.double(epsY),
as.integer(nbzitr),
as.double(zitr0),
as.double(uniqueY0),
as.integer(nvalSPL0),
as.integer(indiceY0),
as.integer(typrisq),
as.integer(risqcom),
as.integer(nz),
as.double(zi),
as.integer(ns0),
as.integer(ng0),
as.integer(nv0),
as.integer(nobs0),
as.integer(nmes0),
as.integer(nbevt),
as.integer(nea0),
as.integer(nw),
as.integer(ncor0),
as.integer(idiag0),
as.integer(idtrunc),
as.integer(logspecif),
as.integer(NPM),
best=as.double(out$best),
ppi=as.double(ppi0),
ppitest=as.double(ppitest0),
resid_m=as.double(resid_m),
resid_ss=as.double(resid_ss),
pred_m_g=as.double(pred_m_g),
pred_ss_g=as.double(pred_ss_g),
predRE=as.double(predRE),
time=as.double(time),
risq_est=as.double(risq_est),
risqcum_est=as.double(risqcum_est),
marker=as.double(marker),
transfY=as.double(transfY),
as.integer(nsim),
Yobs=as.double(Yobs),
statglob=as.double(statglob),
statevt=as.double(statevt),
as.integer(fix0),
as.integer(nfix),
as.double(bfix),
as.integer(estim0),
loglik=as.double(ll))
out$ppi <- post$ppi
out$ppitest <- post$ppitest
out$predRE <- post$predRE
out$Yobs <- post$Yobs
out$resid_m <- post$resid_m
out$resid_ss <- post$resid_ss
out$marker <- post$marker
out$transfY <- post$transfY
out$pred_m_g <- post$pred_m_g
out$pred_ss_g <- post$pred_ss_g
out$risq_est <- post$risq_est
out$risqcum_est <- post$risqcum_est
out$statglob <- post$statglob
out$statevt <- post$statevt
}
## creer best a partir de b et bfix
best <- rep(NA,length(fix0))
best[which(fix0==0)] <- out$best
best[which(fix0==1)] <- bfix
names(best) <- names_best
out$best <- best
NPM <- NPM+nfix
### mettre NA pour les variances et covariances non calculees et 0 pr les prm fixes
if(length(posfix))
{
if(out$conv==3)
{
mr <- NPM-sum(pbH0)-length(posfix)
Vr <- matrix(0,mr,mr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
Vr <- t(Vr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
V <- matrix(NA,NPM,NPM)
V[setdiff(1:NPM,c(which(pbH0==1),posfix)),setdiff(1:NPM,c(which(pbH0==1),posfix))] <- Vr
V[,posfix] <- 0
V[posfix,] <- 0
V <- V[upper.tri(V,diag=TRUE)]
}
else
{
mr <- NPM-length(posfix)
Vr <- matrix(0,mr,mr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
Vr <- t(Vr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
V <- matrix(0,NPM,NPM)
V[setdiff(1:NPM,posfix),setdiff(1:NPM,posfix)] <- Vr
V <- V[upper.tri(V,diag=TRUE)]
}
}
else
{
if(out$conv==3)
{
mr <- NPM-sum(pbH0)
Vr <- matrix(0,mr,mr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
Vr <- t(Vr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
V <- matrix(NA,NPM,NPM)
V[setdiff(1:NPM,which(pbH0==1)),setdiff(1:NPM,which(pbH0==1))] <- Vr
V <- V[upper.tri(V,diag=TRUE)]
}
else
{
V <- out$V
}
}
### remplacer cholesky par varcov dans best
if(nvc>0)
{
ch <- out$best[nprob+nrisqtot+nvarxevt+nef+1:nvc]
if(idiag0==0)
{
U <- matrix(0,nea0,nea0)
U[upper.tri(U,diag=TRUE)] <- ch
z <- t(U) %*% U
out$best[nprob+nrisqtot+nvarxevt+nef+1:nvc] <- z[upper.tri(z,diag=TRUE)]
}
if(idiag0==1)
{
out$best[nprob+nrisqtot+nvarxevt+nef+1:nvc] <- out$best[nprob+nrisqtot+nvarxevt+nef+1:nvc]**2
}
}
else
{
ch <- NA
}
nom.unique[which(nom.unique=="(Intercept)")] <- "intercept"
### predictions des effets aleatoires
if (nea0>0)
{
predRE <- matrix(out$predRE,ncol=nea0,byrow=TRUE)
predRE <- data.frame(unique(IND),predRE)
colnames(predRE) <- c(subject,nom.unique[idea0!=0])
}
### proba des classes a posteroiri et classification
if(ng0>1)
{
ppi <- matrix(out$ppi,ncol=ng0,nrow=ns0,byrow=TRUE)
ppitest <- matrix(out$ppitest,ncol=ng0,byrow=TRUE)
}
else
{
ppi <- matrix(rep(1,ns0),ncol=ng0)
ppitest <- matrix(rep(1,ns0),ncol=ng0)
}
if(!(out$conv %in% c(1,2)))
{
classif <- rep(NA,ns0)
}
else
{
chooseClass <- function(ppi)
{
res <- which.max(ppi)
if(!length(res)) res <- NA
return(res)
}
classif <- apply(ppi,1,chooseClass)
}
ppi <- data.frame(unique(IND),classif,ppi)
temp <- paste("probYT",1:ng0,sep="")
colnames(ppi) <- c(subject,"class",temp)
rownames(ppi) <- 1:ns0
## faire pareil pour ppitest
if(!(out$conv %in% c(1,2)))
{
classif <- rep(NA,ns0)
}
else
{
classif <- apply(ppitest,1,chooseClass)
if(any(!is.finite(ppitest)))
{
classif <- rep(NA,ns0)
iok <- which(is.finite(ppitest[,1]))
classif[iok] <- apply(ppitest[iok,,drop=FALSE],1,which.max)
}
}
ppitest <- data.frame(unique(IND),classif,ppitest)
temp <- paste("probY",1:ng0,sep="")
colnames(ppitest) <- c(subject,"class",temp)
rownames(ppitest) <- 1:ns0
### predictions marginales et subject-specifiques
pred_m_g <- matrix(out$pred_m_g,nrow=nobs0,ncol=ng0)
pred_ss_g <- matrix(out$pred_ss_g,nrow=nobs0,ncol=ng0)
if((out$conv %in% c(1,2)))
{
pred_m <- out$Yobs-out$resid_m
pred_ss <- out$Yobs-out$resid_ss
}
else
{
pred_m <- rep(NA,nobs0)
pred_ss <- rep(NA,nobs0)
}
pred <- data.frame(IND,pred_m,out$resid_m,pred_ss,out$resid_ss,out$Yobs,pred_m_g,pred_ss_g)
temp<-paste("pred_m",1:ng0,sep="")
temp1<-paste("pred_ss",1:ng0,sep="")
colnames(pred)<-c(subject,"pred_m","resid_m","pred_ss","resid_ss","obs",temp,temp1)
if(!is.null(var.time))
{
pred <- data.frame(IND,pred_m,out$resid_m,pred_ss,out$resid_ss,out$Yobs,pred_m_g,pred_ss_g,timeobs)
colnames(pred)<-c(subject,"pred_m","resid_m","pred_ss","resid_ss","obs",temp,temp1,var.time)
}
### risques
risqcum_est <- matrix(out$risqcum_est,nrow=nsim,ncol=ng0*nbevt)
risq_est <- matrix(out$risq_est,nrow=nsim,ncol=ng0*nbevt)
predSurv <- cbind(time,risq_est,risqcum_est)
temp <- paste(paste("event",rep(1:nbevt,each=ng0),".RiskFct",sep=""),1:ng0,sep="")
temp1 <- paste(paste("event",rep(1:nbevt,each=ng0),".CumRiskFct",sep=""),1:ng0,sep="")
colnames(predSurv) <- c("time",temp,temp1)
rownames(predSurv) <- 1:nsim
###estimlink
estimlink <- cbind(out$marker,out$transfY)
colnames(estimlink) <- c(nomY,paste("transf",nomY,sep="."))
### score test
if(out$conv!=1) stats <- rep(NA,nbevt+1)
else
{
stats <- c(out$statglob,out$statevt)
stats[which((!is.finite(stats)) | (stats==9999))] <- NA
}
if(nbevt==1) stats <- stats[1]
### N
N <- NULL
N[1] <- nprob
N[2] <- nrisqtot
N[3] <- nvarxevt
N[4] <- nef
N[5] <- nvc
N[6] <- nw
N[7] <- ncor0
N[8] <- ntrtot0
N[9] <- nobs0
nevent <- rep(0,nbevt)
for(ke in 1:nbevt)
{
nevent[ke] <- length(which(devt==ke))
}
N <- c(N,nevent)
### noms des variables
Names <- list(Xnames=nom.unique,Xnames2=ttesLesVar,Yname=nomY,
ID=nom.subject,Tnames=noms.surv,prior.name=nom.prior,
TimeDepVar.name=nom.timedepvar, pprior.name=nom.pprior)
## levels = modalites des variables dans X0 (si facteurs)
levelsdata <- vector("list", length(ttesLesVar))
levelsfixed <- vector("list", length(ttesLesVar))
levelsrandom <- vector("list", length(ttesLesVar))
levelsmixture <- vector("list", length(ttesLesVar))
levelsclassmb <- vector("list", length(ttesLesVar))
levelssurv <- vector("list", length(ttesLesVar))
names(levelsdata) <- ttesLesVar
names(levelsfixed) <- ttesLesVar
names(levelsrandom) <- ttesLesVar
names(levelsmixture) <- ttesLesVar
names(levelsclassmb) <- ttesLesVar
names(levelssurv) <- ttesLesVar
for(v in ttesLesVar)
{
if(v == "intercept") next
if(is.factor(data[,v]))
{
levelsdata[[v]] <- levels(data[,v])
}
if(length(grep(paste("factor\\(",v,"\\)",sep=""), fixed)))
{
levelsfixed[[v]] <- levels(as.factor(data[,v]))
}
if(length(grep(paste("factor\\(",v,"\\)",sep=""), random)))
{
levelsrandom[[v]] <- levels(as.factor(data[,v]))
}
if(length(grep(paste("factor\\(",v,"\\)",sep=""), mixture)))
{
levelsmixture[[v]] <- levels(as.factor(data[,v]))
}
if(length(grep(paste("factor\\(",v,"\\)",sep=""), classmb)))
{
levelsclassmb[[v]] <- levels(as.factor(data[,v]))
}
if(length(grep(paste("factor\\(",v,"\\)",sep=""), form.surv)))
{
levelssurv[[v]] <- levels(as.factor(data[,v]))
}
}
levels <- list(levelsdata=levelsdata, levelsfixed=levelsfixed, levelsrandom=levelsrandom,
levelsmixture=levelsmixture, levelsclassmb=levelsclassmb,
levelssurv=levelssurv)
cost<-proc.time()-ptm
res <-list(ns=ns0,ng=ng0,idprob=idprob0,idcom=idcom,
idspecif=idspecif,idtdv=idtdv,idg=idg0,idea=idea0,
idcor=idcor0,loglik=out$loglik,best=out$best,V=V,
gconv=out$gconv,conv=out$conv,call=cl,niter=out$niter,
N=N,idiag=idiag0,pred=pred,pprob=ppi,pprobY=ppitest,
predRE=predRE,Names=Names,cholesky=ch,logspecif=logspecif,
estimlink=estimlink,epsY=epsY,linktype=idlink,linknodes=zitr0,
predSurv=predSurv,hazard=list(typrisq,hazardtype,zi,nz),
scoretest=stats,na.action=linesNA,
AIC=2*(length(out$best)-length(posfix)-out$loglik),
BIC=(length(out$best)-length(posfix))*log(ns0)-2*out$loglik,
data=datareturn,levels=levels,var.time=var.time,runtime=cost[3])
class(res) <-c("Jointlcmm")
if(verbose==TRUE) cat("The program took", round(cost[3],2), "seconds \n")
return(res)
}
#' @rdname Jointlcmm
#' @export
jlcmm <- Jointlcmm
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.