R/externVar.R
In CecileProust-Lima/lcmm: Extended Mixed Models Using Latent Classes and Latent Processes
Defines functions
externVar
Documented in
externVar
#' Estimation of secondary regression models after the estimation of a primary latent class model
#'
#' This function fits regression models to relate a latent class structure (stemmed
#' from a latent class model estimated within \code{lcmm} package) with either an external
#' outcome or external class predictors.
#' Two inference techniques are implemented to account for the classification error:
#'
#' - a 2-stage estimation of the joint likelihood of the primary latent class model
#' and the secondary/ external regression;
#'
#' - a regression between the posterior latent class assignment and the external variable
#' which internally corrects for the assignment misclassification.
#'
#' It returns an object from one of the \code{lcmm} package classes.
#'
#' A. DATA STRUCTURE
#'
#' The \code{data} argument must follow specific structure for individual variables,
#' i.e. variables with a unique constant value for each subject. For an individual variable
#' given as external outcome, data value must be present only once per subject,
#' independently of any time variable used in the primary latent class.
#' For an individual variable given as external class predictor,
#' data values must be given for every row of every individual (as usual)
#'
#' B. VARIANCE ESTIMATION
#'
#' Not taking into account first stage variance with specifing \code{"none"} may lead to
#' underestimation of the final variance. When possible, Method \code{"Hessian"}
#' which relies on the combination of Hessians from the primary and secondary
#' model is recommended. However, it may become numerically intensive in the event
#' of very high number of parameters in the primary latent class model. As an
#' alternative, especially in situations with a complex primary model but rather
#' parcimonious secondary model, method \code{"paramBoot"} which implements a
#' parametric bootstrap can be used.
#'
#' @param model an object inheriting from class \code{hlme}, \code{lcmm},
#' \code{Jointlcmm}, \code{multlcmm} or \code{mpjlcmm} giving the primary latent
#' class model.
#' @param fixed optional two sided linear formula object for specifying the
#' fixed-effects in the secondary model with an external outcome variable.
#' 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.
#' @param mixture optional one-sided formula object for the class-specific fixed effects
#' in the model for the external outcome. Among the list of covariates included in fixed,
#' the covariates with class-specific regression parameters are entered in
#' 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 linear formula object for specifying the
#' random-effects on external outcome in the secondary model, if appropriate.
#' By default, no random effect is included.
#' @param subject name of the covariate representing the grouping structure.
#' Even in the absence of a hierarchical structure.
#' @param classmb optional one-sided formula specifying the external predictors of
#' latent class membership to be modelled in the secondary class-membership multinomial
#' logistic model. Covariates are separated by \code{+} on the right of the \code{~}.
#' @param survival optional two-sided formula specifying the external survival part
#' of the model.
#' @param hazard optional family of hazard function assumed for the survival model
#' (Weibull, piecewise or splines)
#' @param hazardtype optional indicator for the type of baseline risk function
#' (Specific, PH or Common)
#' @param hazardnodes optional vector containing interior nodes if \code{splines} or
#' \code{piecewise} is specified for the baseline hazard function in \code{hazard}
#' @param TimeDepVar optional vector specifying the name of the time-depending covariate
#' in the survival model
#' @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
#' @param idiag if appropriate, optional logical for the structure of the variance-covariance
#' matrix of the random-effects in the secondary model.
#' 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 if appropriate, optional logical indicating if the variance-covariance of the
#' random-effects in the secondary model 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 randomY optional logical for including an outcome-specific random intercept.
#' If FALSE no outcome-specific random intercept is added (default). If TRUE independent
#' outcome-specific random intercept with parameterized variance are included
#' @param link optional family of parameterized link functions for the external outcome
#' if appropriate. Defaults to NULL, corresponding to continuous Gaussian distribution
#' (hlme function).
#' @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 cor optional indicator for inclusion of an auto correlated Gaussian process
#' in the latent process linear (latent process) mixed model. Option "BM" indicates
#' a brownian motion with parameterized variance. Option "AR" specifies an
#' autoregressive process of order 1 with parameterized variance and correlation
#' intensity. Each option should be followed by the time variable in brackets as
#' code{cor=BM(time)}. By default, no autocorrelated Gaussian process is added.
#' @param nsim number of points to be used in the estimated link function. By default,
#' nsom=100.
#' @param range optional vector indicating the range of the outcomes (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 data Data frame containing the variables named in
#' \code{fixed}, \code{mixture}, \code{random}, \code{classmb} and \code{subject},
#' for both the current function arguments and the primary model arguments
#' Check \code{details} to get information on the data structure, especially with
#' external outcomes.
#' @param longitudinal only with \code{mpjlcmm} primary models and "twoStageJoint"
#' method: mandatory list containing the longitudinal submodels used in the primary
#' latent class model.
#' @param method character indicating the inference technique to be used:
#' \code{"twoStageJoint"} corresponds to 2-stage estimation. \code{"conditional"}
#' corresponds to the method based on the distribution of Y conditionally to the
#' true latent class membership.
#' @param varest optional character indicating the method to be used to compute the
#' variance of the regression estimates. \code{"none"} does not account for the
#' uncertainty in the primary latent class model, \code{"paramBoot"} computes the
#' total variance using a parametric bootstrap technique, \code{"Hessian"} computes
#' the total Hessian of the joint likelihood (implemented for \code{"twoStageJoint"}
#' method only). Default to \code{"Hessian"} for \code{"twoStageJoint"} method and
#' \code{"paramBoot"} for \code{"conditional"} method.
#' @param M option integer indicating the number of draws for the parametric boostrap
#' when \code{varest="paramBoot"}. Default to 200.
#' @param B optional vector of initial parameter values for the secondary model.
#' If external outcome, the vector has the same structure as a latent class model
#' estimated in the other functions of \code{lcmm} package for the same type of
#' outcome. If external class predictors (of size p), the vector is of length
#' (ng-1)*(1+p). If \code{B=NULL} (by default), internal initial values are selected.
#' @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 secondary model
#' estimation using Marquardt iterative algorithm. Defaults to 100
#' @param posfix optional vector specifying indices in parameter vector B the
#' secondary model that should not be estimated. Default to NULL, all the
#' parameters of the secondary regression 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
#' (can solve non convergence due to estimates at the boundary of the parameter space - usually 0).
#' @param verbose logical indicating whether information about computation should be
#' reported. Default to FALSE.
#' @param nproc the number cores for parallel computation. Default to 1 (sequential mode).
#' @return an object of class \code{externVar} and
#' \code{externSurv} for external survival outcomes,
#' \code{externX} for external class predictors, and
#' \code{hlme}, \code{lcmm}, or \code{multlcmm} for external longitudinal or cross-sectional outcomes.
#' @author Maris Dussartre, Cecile Proust-Lima and Viviane Philipps
#'
#'
#' @examples
#'
#' \dontrun{
#'
#'
#' ###### Estimation of the primary latent class model ######
#'
#' set.seed(1234)
#' PrimMod <- hlme(Ydep1~Time,random=~Time,subject='ID',ng=1,data=data_lcmm)
#' PrimMod2 <- hlme(Ydep1~Time,mixture=~Time,random=~Time,subject='ID',
#' ng=2,data=data_lcmm,B=random(PrimMod))
#'
#' ###### Example 1: Relationship between a latent class structure and #
#' # external class predictors ######
#'
#' # estimation of the secondary multinomial logistic model with total variance
#' # computed with the Hessian
#'
#' XextHess <- externVar(PrimMod2,
#' classmb = ~X1 + X2 + X3 + X4,
#' subject = "ID",
#' data = data_lcmm,
#' method = "twoStageJoint")
#' summary(XextHess)
#'
#' # estimation of a secondary multinomial logistic model with total variance
#' # computed with parametric Bootstrap (much longer). When using the bootstrap
#' # estimator, we recommend running first the analysis with option varest = "none"
#' # which is faster but which underestimates the variance. And then use these values
#' # as initial values when running the model with varest = "paramBoot" to obtain
#' # a valid variance of the parameters.
#'
#' XextNone <- externVar(PrimMod2,
#' classmb = ~X1 + X2 + X3 + X4,
#' subject = "ID",
#' data = data_lcmm,
#' varest = "none",
#' method = "twoStageJoint")
#'
#' XextBoot <- externVar(PrimMod2,
#' classmb = ~X1 + X2 + X3 + X4,
#' subject = "ID",
#' data = data_lcmm,
#' varest = "paramBoot",
#' method = "twoStageJoint",
#' B = XextNone$best)
#' summary(XextBoot)
#'
#'
#' ###### Example 2: Relationship between a latent class structure and #
#' # external outcome (repeatedly measured over time) ######
#'
#' # estimation of the secondary linear mixed model with total variance
#' # computed with the Hessian
#'
#' YextHess = externVar(PrimMod2, #primary model
#' fixed = Ydep2 ~ Time*X1, #secondary model
#' random = ~Time, #secondary model
#' mixture = ~Time, #secondary model
#' subject="ID",
#' data=data_lcmm,
#' method = "twoStageJoint")
#'
#'
#' # estimation of a secondary linear mixed model with total variance
#' # computed with parametric Bootstrap (much longer). When using the bootstrap
#' # estimator, we recommend running first the analysis with option varest = "none"
#' # which is faster but which underestimates the variance. And then use these values
#' # as initial values when running the model with varest = "paramBoot" to obtain
#' # a valid variance of the parameters.
#'
#' YextNone = externVar(PrimMod2, #primary model
#' fixed = Ydep2 ~ Time*X1, #secondary model
#' random = ~Time, #secondary model
#' mixture = ~Time, #secondary model
#' subject="ID",
#' data=data_lcmm,
#' varest = "none",
#' method = "twoStageJoint")
#'
#' YextBoot = externVar(PrimMod2, #primary model
#' fixed = Ydep2 ~ Time*X1, #secondary model
#' random = ~Time, #secondary model
#' mixture = ~Time, #secondary model
#' subject="ID",
#' data=data_lcmm,
#' method = "twoStageJoint",
#' B = YextNone$best,
#' varest= "paramBoot")
#'
#' summary(YextBoot)
#'
#'
#' ###### Example 3: Relationship between a latent class structure and #
#' # external outcome (survival) ######
#'
#' # estimation of the secondary survival model with total variance
#' # computed with the Hessian
#'
#' YextHess = externVar(PrimMod2, #primary model
#' survival = Surv(Tevent,Event)~ X1+mixture(X2), #secondary model
#' hazard="3-quant-splines", #secondary model
#' hazardtype="PH", #secondary model
#' subject="ID",
#' data=data_lcmm,
#' method = "twoStageJoint")
#' summary(YextHess)
#'
#'
#' # estimation of a secondary survival model with total variance
#' # computed with parametric Bootstrap (much longer). When using the bootstrap
#' # estimator, we recommend running first the analysis with option varest = "none"
#' # which is faster but which underestimates the variance. And then use these values
#' # as initial values when running the model with varest = "paramBoot" to obtain
#' # a valid variance of the parameters.
#'
#' YextNone = externVar(PrimMod2, #primary model
#' survival = Surv(Tevent,Event)~ X1+mixture(X2), #secondary model
#' hazard="3-quant-splines", #secondary model
#' hazardtype="PH", #secondary model
#' subject="ID",
#' data=data_lcmm,
#' varest = "none",
#' method = "twoStageJoint")
#'
#' YextBoot = externVar(PrimMod2, #primary model
#' survival = Surv(Tevent,Event)~ X1+mixture(X2), #secondary model
#' hazard="3-quant-splines", #secondary model
#' hazardtype="PH", #secondary model
#' subject="ID",
#' data=data_lcmm,
#' method = "twoStageJoint",
#' B = YextNone$best,
#' varest= "paramBoot")
#'
#' summary(YextBoot)
#'
#' }
#'
#'
#'
#'
#' @export
#'
#'
#'
externVar = function(model,
fixed,
mixture,
random,
subject,
classmb,
survival,
hazard = "Weibull",
hazardtype = "Specific",
hazardnodes = NULL,
TimeDepVar = NULL,
logscale = FALSE,
idiag = FALSE,
nwg = FALSE,
randomY = NULL,
link = NULL,
intnodes = NULL,
epsY = NULL,
cor = NULL,
nsim = NULL,
range = NULL,
data,
longitudinal,
method,
varest,
M = 200,
B,
convB = 0.0001,
convL = 0.0001,
convG = 0.0001,
maxiter = 100,
posfix,
partialH = FALSE,
verbose = FALSE,
nproc = 1){
ptm <- proc.time()
if(missing(model)) stop("model argument must be given")
if(!inherits(model, c("hlme", "lcmm", "multlcmm", "Jointlcmm", "mpjlcmm"))) stop('primary model class must be either "hlme", "lcmm", "multlcmm", "Jointlcmm" or "mpjlcmm"')
if(model$conv == 2) warning("primary model did not fully converge")
if(sum(c(!missing(fixed), !missing(classmb), !missing(survival))) != 1) stop("One and only one in survival, fixed or classmb must be given")
if(missing(method) | !method %in% c("twoStageJoint", "conditional")) stop('Method must be either "twoStageJoint" or "conditional"')
if(model$ng == 1) stop("Primary model does not have latent class structure (ng=1)")
if(method == "twoStageJoint" & missing(varest)) varest = "Hessian"
if(method == "conditional" & missing(varest)) varest = "paramBoot"
if(!varest %in% c("none", "paramBoot", "Hessian")) stop('Variance estimation method "varest" must be either "none", "paramBoot" or "Hessian"')
if(!is.null(link) & missing(fixed)) stop("The argument link is not to be used with external class predictor")
if(missing(posfix)) posfix = c()
cl = match.call()
##Informations about primary model
argumentsIn = as.list(model$call)
funIn = as.character(argumentsIn[[1]])
if(funIn == "jlcmm") funIn <- "Jointlcmm"
if(funIn == "mlcmm") funIn <- "multlcmm"
argumentsIn[[1]] = NULL
ng = model$ng
nIn = length(model$best)
if(is.null(argumentsIn[["classmb"]])){
oldclassmb = ~ 1
} else {
oldclassmb = formula(argumentsIn[["classmb"]])
}
if(!missing(classmb) & oldclassmb != ~1) stop("Primary model already has class predictor")
##number of MB parameters in primary model
nInMB = ncol(model.matrix(oldclassmb, data))*(ng-1)
##Get subject
if(missing(subject)){
if (model$call$subject %in% colnames(data)){
subject = model$call$subject
} else {
stop("The argument subject must be specified if different from the subject argument used in the primary model")
}
}
##Get longitudinal
if(funIn == "mpjlcmm"){
if(missing(longitudinal)) stop("The argument longitudinal is mandatory with a mpjlcmm primary model")
longCall = substitute(longitudinal)
K = length(longitudinal)
for(k in 1:K){
cl_long = as.list(longitudinal[[k]]$call)
if(inherits(longitudinal[[k]], "lcmm")) cl_long[["computeDiscrete"]] = FALSE
longitudinal[[k]]$call = as.call(cl_long)
assign(as.character(longCall[[k+1]]), longitudinal[[k]])
}
}
##nVCIn
if(funIn == "mpjlcmm"){
nVCIn = model$Nprm[3+2*model$K+(1:model$K)]
iVCIn = c()
prev = 0
for(k in 1:model$K){
iVCIn = c(iVCIn, sum(model$Nprm[c(1:3, 3:4*model$K-model$K+k-1)])+prev+1:nVCIn[k])
prev = prev+model$npmK[k]
}
} else if(funIn == "Jointlcmm"){
nVCIn = model$N[5]
iVCIn = sum(model$N[1:4]) + 1:nVCIn
} else if(funIn == "multlcmm"){
nVCIn = model$N[4]
iVCIn = sum(model$N[3]) + 1:nVCIn
} else {
nVCIn = model$N[3]
iVCIn = sum(model$N[1:2]) + 1:nVCIn
}
##pprob with new data
argumentsInEdit = argumentsIn
argumentsInEdit[["data"]] = data
argumentsInEdit[["maxiter"]] = 0
argumentsInEdit[["B"]] = model$best
argumentsInEdit[["verbose"]] = FALSE
modelEdit = do.call(funIn, argumentsInEdit)
pprob = modelEdit$pprob
arguments = list()
##Finding out the number of parameters is needed for all survival external outcome
if(!missing(survival)){
##Informations about secondary outcome model
##number of survival parameters to estimate
##number of events
surv <- cl$survival[[2]]
if(length(surv)==3) #censure droite sans troncature gauche
{
idtrunc <- 0
nom.Tevent <- as.character(surv[2])
nom.Event <- as.character(surv[3])
nom.Tentry <- NULL #si pas de troncature, Tentry=0
noms.surv <- c(nom.Tevent,nom.Event)
}
if(length(surv)==4) #censure droite et troncature
{
idtrunc <- 1
nom.Tentry <- as.character(surv[2])
nom.Tevent <- as.character(surv[3])
nom.Event <- as.character(surv[4])
noms.surv <- c(nom.Tentry,nom.Tevent,nom.Event)
}
Tevent <- getElement(object=data,name=nom.Tevent)
Event <- getElement(object=data,name=nom.Event)
nbevt <- length(attr(do.call("Surv",list(time=Tevent,event=Event,type="mstate")),"states"))
if(nbevt<1) nbevt <- 1
##get number of parameters for baseline functions
hazard <- rep(hazard,length.out=nbevt)
hazardtype <- rep(hazardtype,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"
}
hazWhat = hazard
##when not weibull, keep only the last word of hazard
hazWhat[hazard != 'Weibull'] = sapply(strsplit(hazard[hazard != 'Weibull'], "-"), `[`, 3)
hazN = rep(2, nbevt)
hazN[hazard != 'Weibull'] = sapply(strsplit(hazard[hazard != 'Weibull'], "-"), `[`, 1)
hazN = as.integer(hazN)
hazN = hazN +
(hazWhat == "Weibull")*0 +
(hazWhat == "piecewise")*(-1) +
(hazWhat == "splines")*(2)
hazN = hazN + hazN *
(hazardtype == "Specific")*(ng-1)
##we extract the number of base function parameter with constraints (ie, not PH parameters)
nSurvConstraint = hazN
hazN = hazN +
(hazardtype == "PH")*(ng-1)
##we now have all base functions parameters :
nEstY = sum(hazN)
##get number of parameters for survival covariates
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
}
}
##I extract number of variable in each formula (through model.matrix) excluding intercept
ncols = sapply(c(form.commun, form.cause, form.mixture, form.causek), function(x, data){
mm = model.matrix(x, data)
mm = mm[,-1]
if(is.null(ncol(mm))) {return(1)}
else {return(ncol(mm))}
}, data = data)
##I also need how many parameters each kind of covariate makes
nparam = c(1, nbevt, ng, nbevt*ng)
##we now have all the survival covariates parameters :
nEstX = sum(ncols*nparam)
nEst = nEstY + nEstX
}
##A model structure is needed for all longitudinal external outcome
if(!missing(fixed)){
##Manage inputs
if(missing(mixture)) mixture = ~1
if(missing(random)) random = ~-1
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(length(fixed[[2]]) != 1){
argfunctionStrMod = "multlcmm"
} else if(is.null(link)){
argfunctionStrMod = "hlme"
} else {
argfunctionStrMod = "lcmm"
}
##let's create structure for secondary model
argumentsStrMod = list()
argumentsStrMod[["fixed"]] = fixed
argumentsStrMod[["random"]] = random
argumentsStrMod[["subject"]] = subject
argumentsStrMod[["ng"]] = 1
argumentsStrMod[["idiag"]] = idiag
argumentsStrMod[["randomY"]] = randomY
if(argfunctionStrMod %in% c("lcmm", "multlcmm", "Jointlcmm")){
argumentsStrMod[["link"]] = link
argumentsStrMod[["intnodes"]] = intnodes
}
argumentsStrMod[["epsY"]] = epsY
argumentsStrMod[["cor"]] = substitute(cor)
argumentsStrMod[["nsim"]] = nsim
argumentsStrMod[["range"]] = range
argumentsStrMod[["data"]] = data
argumentsStrMod[["maxiter"]] = 0
argumentsStrMod[["verbose"]] = FALSE
strMod = do.call(argfunctionStrMod, c(argumentsStrMod))
argumentsStrMod[["mixture"]] = mixture
argumentsStrMod[["classmb"]] = ~1
argumentsStrMod[["ng"]] = ng
argumentsStrMod[["nwg"]] = nwg
argumentsStrMod[["B"]] = as.name("strMod")
strMod = do.call(argfunctionStrMod, c(argumentsStrMod))
strMod$best[1:ng+(ng-1)] = mean(strMod$best[1:ng+(ng-1)])# pourquoi?? tous les intercepts a la meme valeur. A revoir
}
##Change general model arguments
if(method == "twoStageJoint"){
##Common argument in twoStageJoint
arguments[["data"]] = data
arguments[["ng"]] = ng
arguments[["subject"]] = subject
##technical options
arguments[["maxiter"]] = maxiter
arguments[["verbose"]] = verbose
arguments[["nproc"]] = nproc
arguments[["convB"]] = convB
arguments[["convL"]] = convL
arguments[["convG"]] = convG
arguments[["partialH"]] = partialH
##primary survival
arguments[["survival"]] = argumentsIn[["survival"]]
arguments[["hazard"]] = argumentsIn[["hazard"]]
arguments[["hazardtype"]] = argumentsIn[["hazardtype"]]
arguments[["hazardnodes"]] = argumentsIn[["hazardnodes"]]
arguments[["hazardrange"]] = argumentsIn[["hazardrange"]]
arguments[["TimeDepVar"]] = argumentsIn[["TimeDepVar"]]
arguments[["logscale"]] = argumentsIn[["logscale"]]
##Primary Jointlcmm needs transformation into lcmm to be put into longitudinal.
if(funIn == "Jointlcmm"){
if(is.null(argumentsIn[["link"]])){
argfunJoint = "hlme"
} else {
argfunJoint = "lcmm"
}
argumentsJoint = argumentsIn
argumentsJoint[["survival"]] = NULL
argumentsJoint[["hazard"]] = NULL
argumentsJoint[["hazardtype"]] = NULL
argumentsJoint[["hazardnodes"]] = NULL
argumentsJoint[["hazardrange"]] = NULL
argumentsJoint[["TimeDepVar"]] = NULL
argumentsJoint[["logscale"]] = NULL
argumentsJoint[["maxiter"]] = 0
argumentsJoint[["mixture"]] = NULL
argumentsJoint[["classmb"]] = NULL
argumentsJoint[["ng"]] = 1
argumentsJoint[["nwg"]] = FALSE
argumentsJoint[["B"]] = NULL
argumentsJoint[["verbose"]] = FALSE
argumentsJoint[["data"]] = data
modNoSurv = do.call(argfunJoint, argumentsJoint)
argumentsJoint[["mixture"]] = argumentsIn[["mixture"]]
argumentsJoint[["classmb"]]= ~1
argumentsJoint[["ng"]] = argumentsIn[["ng"]]
argumentsJoint[["nwg"]] = argumentsIn[["nwg"]]
argumentsJoint[["B"]] = modNoSurv
modNoSurv = do.call(argfunJoint, argumentsJoint)
}
##Yextern survival
if(!missing(survival)){
##manage inputs
if(!is.null(argumentsIn[["survival"]])) stop('secondary survival model is not supported with "twoStageJoint" method if primary model already includes survival')
funOut = "mpjlcmm"
##nOut : nuber of total final parameters
nOut = nIn + nEst
##index
iKeepOut = c(1:nInMB, nInMB+nEst+1:(nIn-nInMB))
iKeepIn = 1:nIn
iEst = nInMB+1:nEst
##input VC, only for keep betaa (iKeepIn)
iVCKeep = iVCIn
##Id of varcov estimates
iVCOut = c()
##list of arguments
arguments[["survival"]] = survival
arguments[["hazard"]] = hazard
arguments[["hazardtype"]] = hazardtype
arguments[["hazardnodes"]] = hazardnodes
arguments[["TimeDepVar"]] = TimeDepVar
arguments[["logscale"]] = logscale
arguments[["posfix"]] = unique(c(iKeepOut, posfix))
##what is in longitudinal ?
if(funIn == "mpjlcmm"){
arguments[["longitudinal"]] = longitudinal
} else {
arguments[["longitudinal"]] = list(model)
}
##initial values
if(missing(B)){
arguments[["B"]][iEst] = rep(0.1, nEst)
} else {
if(length(B) != length(iEst)) stop("B should be of length ", length(iEst))
arguments[["B"]][iEst] = B
}
}
##Yextern longitudinal
if(!missing(fixed)){
funOut = "mpjlcmm"
##Let's change strMod's saved call
strMod$call$data = substitute(data)
## Now join the primary and secondary model
##Informations about secondary outcome model
##number of classmb parameters to remove
nMB = ng-1
##number of remaining parameters to estimate
nStr = length(strMod$best)
nEst = nStr - nMB
##nOut : nuber of total final parameters
nOut = nIn + nEst
##index
iKeepOut = 1:nIn
iKeepIn = iKeepOut
iEst = nIn+1:nEst
##input VC, only for keep betaa (iKeepIn)
iVCKeep = iVCIn
##Id of varcov estimates
if(inherits(strMod, "multlcmm")){
nVCStr = strMod$N[4]
iVCStr = sum(strMod$N[3]) + 1:nVCStr
} else {
nVCStr = strMod$N[3]
iVCStr = sum(strMod$N[1:2]) + 1:nVCStr
}
iVCOut = nIn + iVCStr - nMB
##Liste des arguments
##on fixe nos parametres
arguments[["posfix"]] = unique(c(iKeepOut, posfix))
##On donne les modeles
if(funIn == "mpjlcmm"){
arguments[["longitudinal"]] = c(longitudinal, list(strMod))
} else if(funIn == "Jointlcmm"){
arguments[["longitudinal"]] = list(modNoSurv, strMod)
} else {
arguments[["longitudinal"]] = list(model, strMod)
} ## ici : B= prm modele entre puis prm modele ajoute
##initial values
if(missing(B)){
arguments[["B"]][iEst] = strMod$best[(nMB+1):nStr]
} else {
if(length(B) != length(iEst)) stop("B should be of length ", length(iEst))
arguments[["B"]][iEst] = B
}
}
##X extern
if(!missing(classmb)){
if(!inherits(classmb,"formula")) stop("The argument classmb must be a formula")
funOut = "mpjlcmm"
##nEst : number of MB parameters in output model
nEst1G = ncol(model.matrix(classmb, data))
nEst = nEst1G*(ng-1)
##nOut : nuber of total final parameters
nOut = nIn - nInMB + nEst
##index
iKeepOut = (nEst+1):nOut
iKeepIn = (nInMB+1):nIn
iEst = 1:nEst
##input VC, only for keep betas (iKeepIn)
iVCKeep = iVCIn-nInMB
##Id of varcov estimates (none)
iVCOut = c()
##On recree tous nos arguments
if(funIn == "mpjlcmm"){
arguments[["longitudinal"]] = longitudinal
} else if(funIn == "Jointlcmm"){
arguments[["longitudinal"]] = list(modNoSurv)
} else {
arguments[["longitudinal"]] = list(model)
}
arguments[["classmb"]] = classmb
##on ajoute des valeurs de base pour nos nouveaux estimateurs
arguments[["B"]] = rep(0, nOut)
##initial values
if(missing(B)){
arguments[["B"]][1:ng-1] = model$best[1:ng-1]
} else {
if(length(B) != length(iEst)) stop("B should be of length ", length(iEst))
arguments[["B"]][iEst] = B
}
##on fixe nos parametres
arguments[["posfix"]] = unique(c(posfix, iKeepOut))
}
}
if(method == "conditional"){
##Yextern survival
if(!missing(survival)){
## we need it all in a function in order to be able to use parametric bootstrap later on
##remaining parameters to estimate
iEst = 1:nEst
iKeepIn = 1:nIn
iKeepOut = 1:nIn+nEst+2
nOut = nEst+2
##Id of varcov estimates to keep in bootstrap
iVCKeep = iVCIn
iVCOut = c()
nVCIn = 0
##We need what is inside of longitudinal to still exist in the worker
argumentsIn[["longitudinal"]] = eval(argumentsIn[["longitudinal"]])
conditionalS = function(model,
data,
survival,
hazard,
hazardtype,
hazardnodes = NULL,
TimeDepVar = NULL,
logscale,
subject,
ng,
B,
link,
iEst,
maxiter,
verbose,
argumentsIn,
funIn,
nproc,
convB,
convL,
convG){
argumentsInEdit = argumentsIn
argumentsInEdit[["B"]] = B[iKeepOut]
argumentsInEdit[["maxiter"]] = 0
argumentsInEdit[["verbose"]] = F
argumentsInEdit[["nproc"]] = nproc
argumentsInEdit[["data"]] = data
model = do.call(funIn, argumentsInEdit)
B = B[iEst]
predCl = predictClass(model, data)
##First : let's compute P(C|\tilde C) (\tilde C : A)
pAlY = sapply(1:ng, function(g){
return(as.numeric(predCl[,2] == g))
})
pClY = as.matrix(predCl[,3:(2+ng)])
if(any(is.nan(pClY))) stop("NaN in posterior classification probability")
pA = apply(pAlY, 2, mean)
pClA = t(pAlY)%*%pClY/(model$ns*pA)
if(det(pClA) == 0 | is.na(det(pClA))) stop("Computed error matrix is singular. One class might be empty")
##Then : let's add it to dataset for each individual
indivProb = pAlY%*%pClA
indivProb = cbind(predCl[,1], indivProb)
colnames(indivProb)[1] = subject
data = merge(data, indivProb, by = subject)
##We need dummy Y for the model to run
data$dummyY = 1
##With some type of input model, prob does not have the same name
if(!"prob1" %in% colnames(data)){
for(i in 1:ng){
data[[paste0("prob", i)]] = data[[paste0("probYT", i)]]
}
}
##Finally : we need to build the model for ppriors !
arguments = list()
arguments[["data"]] = data
arguments[["fixed"]] = dummyY~1
arguments[["mixture"]] = ~-1
arguments[["survival"]] = survival
arguments[["hazard"]] = hazard
arguments[["hazardtype"]] = hazardtype
arguments[["hazardnodes"]] = hazardnodes
arguments[["TimeDepVar"]] = TimeDepVar
arguments[["logscale"]] = logscale
arguments[["subject"]] = subject
arguments[["classmb"]] = ~-1
arguments[["ng"]] = ng
arguments[["pprior"]] = c("prob1", "prob2")
arguments[["posfix"]] = 1:2+length(iEst)
##technical options
arguments[["maxiter"]] = maxiter
arguments[["verbose"]] = verbose
arguments[["nproc"]] = nproc
arguments[["convB"]] = convB
arguments[["convL"]] = convL
arguments[["convG"]] = convG
arguments[["B"]] = c(B[iEst], 1, 0.000001)
res = do.call("Jointlcmm", arguments)
res$call = match.call()
return(res)
}
##we need to build the model
arguments[["data"]] = data
arguments[["survival"]] = survival
arguments[["hazard"]] = hazard
arguments[["hazardtype"]] = hazardtype
arguments[["hazardnodes"]] = hazardnodes
arguments[["TimeDepVar"]] = TimeDepVar
arguments[["logscale"]] = logscale
arguments[["model"]] = model
arguments[["subject"]] = subject
arguments[["ng"]] = ng
arguments[["link"]] = link
arguments[["iEst"]] = iEst
##technical options
arguments[["maxiter"]] = maxiter
arguments[["verbose"]] = verbose
arguments[["argumentsIn"]] = argumentsIn
arguments[["funIn"]] = funIn
arguments[["nproc"]] = nproc
arguments[["convB"]] = convB
arguments[["convL"]] = convL
arguments[["convG"]] = convG
##on ajoute des valeurs de base pour nos nouveaux estimateurs
arguments[["B"]] = rep(0.1, nIn+nOut)
##initial values
if(!missing(B)){
if(length(B) != length(iEst)) stop("B should be of length ", length(iEst))
arguments[["B"]][iEst] = B
}
funOut = "conditionalS"
}
##Yextern longitudinal
if(!missing(fixed)){
## we need it all in a function in order to be able to use parametric bootstrap later on
##number of classmb parameters to remove
nMB = ng-1
##number of remaining parameters to estimate
nStr = length(strMod$best)
nEst = nStr - nMB
iEst = 1:nEst
iKeepIn = 1:nIn
iKeepOut = 1:nIn+nEst
nOut = nEst
##Id of varcov estimates to keep in bootstrap
iVCKeep = iVCIn
##Id of varcov estimates
if(inherits(strMod, "multlcmm")){
nVCStr = strMod$N[4]
iVCStr = sum(strMod$N[3]) + 1:nVCStr
} else {
nVCStr = strMod$N[3]
iVCStr = sum(strMod$N[1:2]) + 1:nVCStr
}
iVCOut = iVCStr - nMB
nVCIn = 0
##We need what is inside of longitudinal to still exist in the worker
argumentsIn[["longitudinal"]] = eval(argumentsIn[["longitudinal"]])
conditional = function(model,
data,
fixed,
random,
idiag,
nwg,
randomY = NULL,
link,
intnodes = NULL,
epsY = NULL,
cor = NULL,
nsim = NULL,
range = NULL,
subject,
mixture,
ng,
B,
iEst,
maxiter,
verbose,
argumentsIn,
funIn,
nproc,
convB,
convL,
convG){
argumentsInEdit = argumentsIn
argumentsInEdit[["B"]] = B[iKeepOut] ## B : prm model a estimer, prm model entree
argumentsInEdit[["maxiter"]] = 0
argumentsInEdit[["verbose"]] = F
argumentsInEdit[["nproc"]] = nproc
argumentsInEdit[["data"]] = data
model = do.call(funIn, argumentsInEdit)
B = B[iEst]
##predCl = predictClass(model, data) # pas utile # chgmt Viviane
predCl <- model$pprob
##First : let's compute P(C|\tilde C) (\tilde C : A)
pAlY = sapply(1:ng, function(g){
return(as.numeric(predCl[,2] == g))
})
pClY = as.matrix(predCl[,3:(2+ng)])
if(any(is.nan(pClY))) stop("NaN in posterior classification probability")
pA = apply(pAlY, 2, mean)
pClA = t(pAlY)%*%pClY/(model$ns*pA)
if(det(pClA) == 0 | is.na(det(pClA))) stop("Computed error matrix is singular. One class might be empty")
##Then : let's add it to dataset for each individual
indivProb = pAlY%*%pClA
indivProb = cbind(predCl[,1], indivProb)
colnames(indivProb)[1] = subject
data = merge(data, indivProb, by = subject)
##Finally : we need to build the model for ppriors !
arguments = list()
if(length(fixed[[2]]) != 1){
funOut = "multlcmm"
arguments[["link"]] = link
arguments[["intnodes"]] = intnodes
} else if(missing(link)){
funOut = "hlme"
arguments[["link"]] = NULL
} else {
funOut = "lcmm"
arguments[["link"]] = link
arguments[["intnodes"]] = intnodes
}
##With some type of input model, prob does not have the same name
if(!"prob1" %in% colnames(data)){
for(i in 1:ng){
data[[paste0("prob", i)]] = data[[paste0("probYT", i)]]
}
}
arguments[["data"]] = data
arguments[["fixed"]] = fixed
arguments[["random"]] = random
arguments[["idiag"]] = idiag
arguments[["nwg"]] = nwg
arguments[["randomY"]] = randomY
arguments[["epsY"]] = epsY
arguments[["cor"]] = substitute(cor)
arguments[["nsim"]] = nsim
arguments[["range"]] = range
arguments[["subject"]] = subject
arguments[["mixture"]] = mixture
arguments[["classmb"]] = ~-1
arguments[["ng"]] = ng
arguments[["pprior"]] = c("prob1", "prob2")
##technical options
arguments[["maxiter"]] = maxiter
arguments[["verbose"]] = verbose
arguments[["nproc"]] = nproc
arguments[["convB"]] = convB
arguments[["convL"]] = convL
arguments[["convG"]] = convG
arguments[["B"]] = B ## ici B = vi du modele secondaire
res = do.call(funOut, arguments)
res$call = match.call()
return(res)
}
##we need to build the model
arguments[["data"]] = data
arguments[["model"]] = model
arguments[["fixed"]] = fixed
arguments[["random"]] = random
arguments[["idiag"]] = idiag
arguments[["nwg"]] = nwg
arguments[["randomY"]] = randomY
arguments[["link"]] = link
arguments[["intnodes"]] = intnodes
arguments[["epsY"]] = epsY
arguments[["cor"]] = substitute(cor)
arguments[["nsim"]] = nsim
arguments[["range"]] = range
arguments[["subject"]] = subject
arguments[["mixture"]] = mixture
arguments[["ng"]] = ng
arguments[["iEst"]] = iEst
##technical options
arguments[["maxiter"]] = maxiter
arguments[["verbose"]] = verbose
arguments[["argumentsIn"]] = argumentsIn
arguments[["funIn"]] = funIn
arguments[["nproc"]] = nproc
arguments[["convB"]] = convB
arguments[["convL"]] = convL
arguments[["convG"]] = convG
##on ajoute des valeurs de base pour nos nouveaux estimateurs
arguments[["B"]] = rep(0, nIn+nOut) ## ici B = valeurs pour les 2 modeles
##initial values
if(missing(B)){
arguments[["B"]][iEst] = strMod$best[(nMB+1):nStr]
} else {
if(length(B) != length(iEst)) stop("B should be of length ", length(iEst))
arguments[["B"]][iEst] = B
}
funOut = "conditional"
}
##Xextern
if(!missing(classmb)){
##we need it all in a function in order to be able to use parametric bootstrap later
##nEst : number of MB parameters in output model
nEst1G = ncol(model.matrix(classmb, data))
nEst = nEst1G*(ng-1)
iEst = 1:nEst
iKeepIn = 1:nIn
iKeepOut = 1:nIn+nEst
nOut = nEst
##Id of varcov estimates to keep in bootstrap
iVCKeep = iVCIn
##Id of varcov estimates (none)
iVCOut = c()
nVCIn = 0
##We need what is inside of longitudinal to still exist in the worker
argumentsIn[["longitudinal"]] = eval(argumentsIn[["longitudinal"]])
conditionalX = function(classmb,
data,
ng,
B,
iKeepIn,
iEst,
argumentsIn,
funIn,
nproc,
maxiter,
verbose){
argumentsInEdit = argumentsIn
argumentsInEdit[["B"]] = B[iKeepOut]
argumentsInEdit[["maxiter"]] = 0
argumentsInEdit[["verbose"]] = F
argumentsInEdit[["nproc"]] = nproc
argumentsInEdit[["data"]] = data
model = do.call(funIn, argumentsInEdit)
B = B[iEst]
predCl = predictClass(model, data)
##First : let's compute P(\tilde C|C) (\tilde C : A)
pAlY = sapply(1:ng, function(g){
return(as.numeric(predCl[,2] == g))
})
pClY = as.matrix(predCl[,3:(2+ng)])
if(any(is.nan(pClY))) stop("NaN in posterior classification probability")
betas = c(model$best[1:(ng-1)], 0)
pC = sapply(betas, function(b, betas){
exp(b)/sum(exp(betas))
}, betas=betas)
pAlC = t(pClY)%*%pAlY/(model$ns*pC)
if(det(pAlC) == 0 | is.na(det(pAlC))) stop("Computed error matrix is singular. One class might be empty")
##Then : let's add the classification to the dataset for each individual
indivProb = pAlY%*%t(pAlC)
indivProb = cbind(predCl[,1], indivProb)
colnames(indivProb) = c(subject, paste0("class", 1:ng))
data = merge(data, indivProb, by = subject)
for(id in unique(data[[subject]])){
if(sum(data[[subject]] == id) > 1) data = data[-which(data[[subject]] == id)[-1],]
}
##negative log likelihood function
nLL = function(beta, y, X) {
beta = matrix(ncol = ncol(y)-1, byrow = T, beta)
denom <- apply(exp(X%*%beta),1,sum) + 1
num_mat <- y*cbind(exp(X%*%beta), 1)
num <- apply(num_mat,1,sum)
vrais <- sum(log(num/denom))
return(-vrais)
}
##frame
argmf = list(
formula = classmb,
data = data
)
mf = do.call(model.frame, argmf)
ns = nrow(mf)
y = data[, paste0("class", 1:ng)]
y = as.matrix(sapply(y, as.numeric))
nBy = ncol(y)-1
X = model.matrix(classmb, data)
nEst = ncol(X)*nBy
opt = mla(b=B, fn=nLL, y=y, X=X, print.info = verbose, nproc = nproc, maxiter = maxiter)
namesX = c("intercept", colnames(X)[colnames(X) != "(Intercept)"])
names(opt$b) = c(sapply(namesX, FUN = function(i){
return(paste0(i, " ", colnames(y)[-ng]))
}))
Names = list(Xnsnames = namesX,
ID = subject)
N = c(nEst)
res = list(best = opt$b,
V = opt$v,
conv = opt$istop,
loglik = -opt$fn.value,
ns = length(unique(data[[subject]])),
ng = ng,
idprob = rep(1, ncol(X)),
nv2 = ncol(X),
gconv = c(opt$ca, opt$cb, opt$rdm),
pprob = NULL,
Names = Names,
N = N,
call = match.call())
return(res)
}
##Finally : we need to build the model arguments
arguments[["classmb"]] = classmb
arguments[["data"]] = data
arguments[["ng"]] = ng
arguments[["iKeepIn"]] = iKeepIn
arguments[["iEst"]] = iEst
arguments[["argumentsIn"]] = argumentsIn
arguments[["funIn"]] = funIn
arguments[["nproc"]] = nproc
arguments[["maxiter"]] = maxiter
arguments[["verbose"]] = verbose
##on ajoute des valeurs de base pour nos nouveaux estimateurs
arguments[["B"]] = rep(0, nIn+nOut)
##initial values
if(!missing(B)){
if(length(B) != length(iEst)) stop("B should be of length ", length(iEst))
arguments[["B"]][iEst] = B
}
funOut = "conditionalX"
}
}
if(varest != "paramBoot"){
arguments[["B"]][iKeepOut] = model$best[iKeepIn]
if(verbose){cat("Model estimation...\n\n")}
##Model Estimation
modOut = do.call(funOut, c(arguments))
}
if(varest == "Hessian"){
if(method != "twoStageJoint") stop("Hessian variance estimation method only avaliable for 'twoStageJoint' method")
if(verbose){cat("Variance estimation...\n\n")}
nb11 = length(model$best)
V11 = matrix(0, nb11, nb11)
V11[upper.tri(V11, diag=T)] = model$V
V11[lower.tri(V11, diag=F)] = t(V11)[lower.tri(V11, diag=F)]
n1 = model$ns
V11 = V11[iKeepIn, iKeepIn]
nb22 = length(modOut$best)
V22 = matrix(0, nb22, nb22)
V22[upper.tri(V22, diag=T)] = modOut$V
V22[lower.tri(V22, diag=F)] = t(V22)[lower.tri(V22, diag=F)]
saveV22 = V22
V22 = V22[iEst, iEst]
n2 = modOut$ns
modOut
if(nproc == 1){
I12 = -hessienne(modOut)
} else {
I12 = -hessienne(modOut, method = "deriva", nproc = nproc)
}
I12 = I12[iEst, iKeepOut]
V = V22*n2 + (V22*n2) %*% (I12/n2) %*% ((n2/n1)*(V11*n1)) %*% t(I12/n2) %*% (V22*n2)
V = V/n2
saveV22[iEst, iEst] = V
V = saveV22
modOut$V = V[upper.tri(V, diag = TRUE)]
}
##Get Bootstrap Models
if(varest == "paramBoot"){
if(verbose){cat("Bootstrap estimation...\n\n")}
est = estimates(model)
Vin = matrix(0, length(est), length(est))
Vin[upper.tri(Vin, diag = T)] = model$V
Vin[lower.tri(Vin, diag = F)] = t(Vin)[lower.tri(t(Vin), diag=F)]
est = est[iKeepIn]
Vin = Vin[iKeepIn,iKeepIn]
coefss = rmvnorm(M, est, Vin)
coefss = as.data.frame(coefss)
colnames(coefss) = names(model$best)[iKeepIn]
##we just need to build back varcov into the coefs instead of cholesky matrix
coefss = apply(coefss, 1, function(coefs, model, data, iVCKeep){
if(funIn == "mpjlcmm"){
varcovMods = longitudinal
} else {
varcovMods = list(model)
}
chols = coefs[iVCKeep]
model$cholesky = chols
varcov = c()
countChol = 0
for(varcovMod in varcovMods){
ncolRandMod = ncol(model.matrix(formula(varcovMod$call$random), data))
if(varcovMod$idiag){
nChol = ncolRandMod-as.integer(funIn == "multlcmm")
vc = chols[1:nChol+countChol]
varcov = c(varcov, (vc^2))
} else {
nChol = ncolRandMod*(ncolRandMod+1)/2-as.integer(funIn == "multlcmm")
##cholMatrix
cholMatrix = matrix(0, ncolRandMod, ncolRandMod)
cholsToMatrix = chols[1:nChol+countChol]
if(funIn == "multlcmm") cholsToMatrix = c(1, cholsToMatrix)
cholMatrix[upper.tri(cholMatrix, diag = T)] = cholsToMatrix
vc = t(cholMatrix)%*%cholMatrix
varcov = c(varcov, vc[upper.tri(vc, diag = T)])
if(funIn == "multlcmm") varcov = varcov[-1]
}
countChol = countChol + nChol
}
coefs[iVCKeep] = varcov
return(coefs)
}, model = model, data = data, iVCKeep = iVCKeep)
if(nproc > 1)
{
clust <- parallel::makeCluster(nproc)
##load all loaded packages
packages = loadedNamespaces()
for(pack in packages){
clusterExport(clust, "pack", environment())
clusterEvalQ(clust, require(pack, character.only = T))
}
survivalMissing = missing(survival)
fixedMissing = missing(fixed)
modOuts <- parApply(clust, coefss, 2, function(coefs, arguments, iKeepOut, funOut, iEst, survivalMissing, fixedMissing, logscale){
arguments[["B"]][iKeepOut] = coefs
arguments[["nproc"]] = 1
##Model Estimation
modOut = do.call(funOut, c(arguments))
######## chgmt Viviane ########
## if(!fixedMissing){
## ##cholesky not varcov as output in best
## modOut$best = estimates(modOut)
## ## abs value of link prm :
## n = length(modOut$best)
## if(inherits(modOut, "lcmm")){ #For lcmm
## nLink <- n - sum(modOut$N[c(1:4,6)])
## if(modOut$linktype == 2){ #with spline link
## ##nSpl = n-sum(modOut$N[1:4]) #count number of link function parameters (the only one not in N)
## ##modOut$best[n-2:nSpl+1+1] = abs(modOut$best[n-2:nSpl+1+1]) #all but the first
## modOut$best[sum(modOut$N[1:4])+2:nLink] = abs(modOut$best[sum(modOut$N[1:4])+2:nLink]) #all but the first
## } else { #rest of lcmm
## modOut$best[sum(modOut$N[1:4])+nLink] = abs(modOut$best[sum(modOut$N[1:4])+nLink]) #only the last one
## }
## } else if (inherits(modOut, "hlme")){ #hlme
## modOut$best[n] = abs(modOut$best[n]) #only the last one
## } else if(inherits(modOut, "multlcmm")){
## nPreLink = sum(modOut$N[3:8]) #number of parameters before the one for the link function
## numSPL = 0
## for (ny in 1:modOut$N[8]){
## if (modOut$linktype[ny]==0) nLink = 2
## if (modOut$linktype[ny]==1) nLink = 4
## if (modOut$linktype[ny]==2){
## numSPL <- numSPL+1
## nLink = modOut$nbnodes[numSPL]+2
## modOut$best[nPreLink+2:nLink] = abs(modOut$best[nPreLink+2:nLink])
## } else {
## modOut$best[nPreLink+nLink] = abs(modOut$best[nPreLink+nLink])
## }
## nPreLink = nPreLink+nLink
## }
## } else if(inherits(modOut, "mpjlcmm")){
## #Residual Error : need to be the same sign across bootstrap iterations
## nPre = sum(modOut$N[1:(2+modOut$nbevt)])
## sumny = 0
## for(k in 1:modOut$K){
## if(modOut$contrainte[k] == 2){ ## attention : length(nrisq) = nbevt!!
## nPre = nPre + sum(modOut$Nprm[2+modOut$nbevt+1:7*modOut$K-modOut$K+k]) #add number of parameters before the one for the link function
## } else {
## nPre = nPre + sum(modOut$Nprm[2+modOut$nbevt+1:8*modOut$K-modOut$K+k]) #add the number of parameters for this K
## }
## for(y in 1:modOut$ny[k]){
## sumny = sumny+1
## if(modOut$contrainte[k] == 1 & modOut$linktype[sumny] == 2){ #for lcmm with spl link
## nSpl = sum(modOut$Nprm[2+modOut$nbevt+7*modOut$K+sumny]) #count number of link function parameters (the only one not in N)
## modOut$best[nPre-2:nSpl+1+1] = abs(modOut$best[nPre-2:nSpl+1+1]) #all but the first
## } else if (modOut$contrainte[k] == 0 | modOut$contrainte[k] == 1){ #hlme ou le reste de lcmm
## modOut$best[nPre] = abs(modOut$best[nPre]) #only the last one
## } else if (modOut$contrainte[k] == 2) {
## numSPL = 0
## if (modOut$linktype[sumny]==0) nLink = 2
## if (modOut$linktype[sumny]==1) nLink = 4
## if (modOut$linktype[sumny]==2){
## numSPL <- numSPL+1
## nLink = modOut$call$longitudinal[[k]]$nbnodes[numSPL]+2
## modOut$best[nPre+2:nLink] = abs(modOut$best[nPre+2:nLink])
## } else {
## modOut$best[nPre+nLink] = abs(modOut$best[nPre+nLink])
## }
## nPre = nPre+nLink
## }
## }
## }
## }
## }
## ##nw & cor & alea
## idnw <- NULL
## idalea <- NULL
## idcor <- NULL
## if(inherits(modOut, "multlcmm")){#with multlcmm
## if(modOut$N[5] != 0) idnw = sum(modOut$N[3:4])+1:modOut$N[5]
## if(modOut$N[6] != 0) idalea = sum(modOut$N[3:5])+1:modOut$N[6]
## if(modOut$N[7] != 0) idcor = sum(modOut$N[3:6])+modOut$N[7]
## } else if(inherits(modOut, "mpjlcmm")){#with mpjlcmm
## nPre = sum(modOut$N[1:3])
## idnw = c()
## idcor = c()
## idalea = c()
## for(k in 1:modOut$K){
## nprenw = nPre+sum(modOut$Nprm[c(3+1:3*modOut$K-modOut$K+k)])
## if(modOut$Nprm[3+4*modOut$K-modOut$K+k] != 0){
## idnw = c(idnw, nprenw+1:modOut$Nprm[3+4*modOut$K-modOut$K+k])
## }
## nprecor = nPre+sum(modOut$Nprm[c(3+1:4*modOut$K-modOut$K+k)])
## if(modOut$Nprm[3+5*modOut$K-modOut$K+k] != 0){
## idcor = c(idcor, nprecor+modOut$Nprm[3+5*modOut$K-modOut$K+k])
## }
## nprealea = nPre+sum(modOut$Nprm[c(3+1:6*modOut$K-modOut$K+k)])
## if(modOut$Nprm[3+7*modOut$K-modOut$K+k] != 0){
## idalea = c(idalea, nprealea+1:modOut$Nprm[3+7*modOut$K-modOut$K+k])
## }
## nPre = nPre+sum(modOut$Nprm[c(3+1:8*modOut$K-modOut$K+k)])
## }
## } else {#with lcmm & hlme
## if(modOut$N[4] != 0) idnw = sum(modOut$N[1:3])+1:modOut$N[4]
## if(modOut$N[5] != 0) idcor = sum(modOut$N[1:4])+modOut$N[5]
## }
## modOut$best[idnw] = abs(modOut$best[idnw])
## modOut$best[idcor] = abs(modOut$best[idcor])
## modOut$best[idalea] = abs(modOut$best[idalea])
## #Survival Base Function : need to be the same sign across bootstrap iterations
## if(!survivalMissing & !logscale){
## iSurvConstraint = 1:nSurvConstraint+modOut$N[1]
## modOut$best[iSurvConstraint] = abs(modOut$best[iSurvConstraint])
## }
######## fin chgmt Viviane ########
if(inherits(modOut, c("hlme", "lcmm", "multlcmm", "Jointlcmm", "mpjlcmm"))) modOut <- absprm(modOut) ## remplace ce que j'ai mis en commentaire
## si model multinomial, pas besoin de valeur absolue dans les prm
return(modOut)
}, arguments = arguments, iKeepOut = iKeepOut, funOut = funOut, iEst = iEst, survivalMissing = survivalMissing, fixedMissing = fixedMissing, logscale = logscale)
parallel::stopCluster(clust)
##format output
bests = as.data.frame(matrix(NA, nrow = nEst, ncol = M))
Vs = list()
conv = c()
modOut = modOuts[[M]]
for (i in 1:M){
##output V and betas
bests[,i] = modOut$best[iEst]
V = matrix(NA, nOut, nOut)
V[upper.tri(V, diag = T)] = modOut$V
V[lower.tri(V, diag = F)] = t(V)[lower.tri(t(V), diag = F)]
Vs = c(Vs, list(V[iEst, iEst]))
conv = c(conv, modOut$conv)
}
} else {
##estimate final models, extract usefull information
bests = as.data.frame(matrix(NA, nrow = nEst, ncol = M))
Vs = list()
conv = c()
for(i in 1:M){
if(verbose) cat("==================== Bootstrap Iteration", i, "====================\n")
arguments[["B"]][iKeepOut] = coefss[,i]
##Model Estimation
if(verbose){
modOut = do.call(funOut, c(arguments))
} else {
captured_log = capture.output({modOut = do.call(funOut, c(arguments))})
}
######## chgmt Viviane ########
## if(!missing(fixed)){
## #cholesky not varcov as output in best
## modOut$best = estimates(modOut)
## n = length(modOut$best)
## if(inherits(modOut, "lcmm")){ #For lcmm
## if(modOut$linktype == 2){ #with spline link
## nSpl = n-sum(modOut$N[1:4]) #count number of link function parameters (the only one not in N)
## modOut$best[n-2:nSpl+1+1] = abs(modOut$best[n-2:nSpl+1+1]) #all but the first
## } else { #rest of lcmm
## modOut$best[n] = abs(modOut$best[n]) #only the last one
## }
## } else if (inherits(modOut, "hlme")){ #hlme
## modOut$best[n] = abs(modOut$best[n]) #only the last one
## } else if(inherits(modOut, "multlcmm")){
## nPreLink = modOut$N[1:5] #number of parameters before the one for the link function
## numSPL = 0
## for (ny in 1:modOut$N[8]){
## if (modOut$linktype[ny]==0) nLink = 2
## if (modOut$linktype[ny]==1) nLink = 4
## if (modOut$linktype[ny]==2){
## numSPL <- numSPL+1
## nLink = modOut$nbnodes[numSPL]+2
## modOut$best[nPreLink+2:nLink] = abs(modOut$best[nPreLink+2:nLink])
## } else {
## modOut$best[nPreLink+nLink] = abs(modOut$best[nPreLink+nLink])
## }
## nPreLink = nPreLink+nLink
## }
## } else if(inherits(modOut, "mpjlcmm")){
## #Residual Error : need to be the same sign across bootstrap iterations
## nPre = sum(modOut$N[1:3])
## sumny = 0
## for(k in 1:modOut$K){
## if(modOut$contrainte[k] == 2){
## nPre = nPre + sum(modOut$Nprm[3+1:7*modOut$K-modOut$K+k]) #add number of parameters before the one for the link function
## } else {
## nPre = nPre + sum(modOut$Nprm[3+1:8*modOut$K-modOut$K+k]) #add the number of parameters for this K
## }
## for(y in 1:modOut$ny[k]){
## sumny = sumny+1
## if(modOut$contrainte[k] == 1 & modOut$linktype[sumny] == 2){ #for lcmm with spl link
## nSpl = sum(modOut$Nprm[3+7*modOut$K+1*modOut$K-modOut$K+k]) #count number of link function parameters (the only one not in N)
## modOut$best[nPre-2:nSpl+1+1] = abs(modOut$best[nPre-2:nSpl+1+1]) #all but the first
## } else if (modOut$contrainte[k] == 0 | modOut$contrainte[k] == 1){ #hlme ou le reste de lcmm
## modOut$best[nPre] = abs(modOut$best[nPre]) #only the last one
## } else if (modOut$contrainte[k] == 2) {
## numSPL = 0
## if (modOut$linktype[sumny]==0) nLink = 2
## if (modOut$linktype[sumny]==1) nLink = 4
## if (modOut$linktype[sumny]==2){
## numSPL <- numSPL+1
## nLink = modOut$call$longitudinal[[k]]$nbnodes[numSPL]+2
## modOut$best[nPre+2:nLink] = abs(modOut$best[nPre+2:nLink])
## } else {
## modOut$best[nPre+nLink] = abs(modOut$best[nPre+nLink])
## }
## nPre = nPre+nLink
## }
## }
## }
## }
## }
## ##nw & cor & alea
## idnw <- NULL
## idalea <- NULL
## idcor <- NULL
## if(inherits(modOut, "multlcmm")){#with multlcmm
## if(modOut$N[5] != 0) idnw = sum(modOut$N[3:4])+1:modOut$N[5]
## if(modOut$N[6] != 0) idalea = sum(modOut$N[3:5])+1:modOut$N[6]
## if(modOut$N[7] != 0) idcor = sum(modOut$N[3:6])+modOut$N[7]
## } else if(inherits(modOut, "mpjlcmm")){#with mpjlcmm
## nPre = sum(modOut$N[1:3])
## idnw = c()
## idcor = c()
## idalea = c()
## for(k in 1:modOut$K){
## nprenw = nPre+sum(modOut$Nprm[c(3+1:3*modOut$K-modOut$K+k)])
## if(modOut$Nprm[3+4*modOut$K-modOut$K+k] != 0){
## idnw = c(idnw, nprenw+1:modOut$Nprm[3+4*modOut$K-modOut$K+k])
## }
## nprecor = nPre+sum(modOut$Nprm[c(3+1:4*modOut$K-modOut$K+k)])
## if(modOut$Nprm[3+5*modOut$K-modOut$K+k] != 0){
## idcor = c(idcor, nprecor+modOut$Nprm[3+5*modOut$K-modOut$K+k])
## }
## nprealea = nPre+sum(modOut$Nprm[c(3+1:6*modOut$K-modOut$K+k)])
## if(modOut$Nprm[3+7*modOut$K-modOut$K+k] != 0){
## idalea = c(idalea, nprealea+1:modOut$Nprm[3+7*modOut$K-modOut$K+k])
## }
## nPre = nPre+sum(modOut$Nprm[c(3+1:8*modOut$K-modOut$K+k)])
## }
## } else {#with lcmm
## if(inherits(modOut, "lcmm")){
## if(modOut$N[4] != 0) idnw = sum(modOut$N[1:3])+1:modOut$N[4]
## if(modOut$N[6] != 0) idcor = sum(modOut$N[1:4])+ntr+modOut$N[5]
## }
## else { #hlme
## if(inherits(modOut, "hlme")){
## if(modOut$N[4] != 0) idnw = sum(modOut$N[1:3])+1:modOut$N[4]
## if(modOut$N[5] != 0) idcor = sum(modOut$N[1:4])+modOut$N[5]
## }
## }
## }
## modOut$best[idnw] = abs(modOut$best[idnw])
## modOut$best[idcor] = abs(modOut$best[idcor])
## modOut$best[idalea] = abs(modOut$best[idalea])
## #Survival Base Function constraint : need to be the same sign across bootstrap iterations
## if(!missing(survival) & !logscale){
## iSurvConstraint = 1:nSurvConstraint+modOut$N[1]
## modOut$best[iSurvConstraint] = abs(modOut$best[iSurvConstraint])
## }
######## chgmt Viviane ########
if(inherits(modOut, c("hlme", "lcmm", "multlcmm", "Jointlcmm", "mpjlcmm"))) modOut <- absprm(modOut) ## remplace ce que j'ai mis en commentaire
## si model multinomial, pas besoin de valeur absolue dans les prm
##output V and betas
bests[,i] = modOut$best[iEst]
V = matrix(NA, nOut, nOut)
V[upper.tri(V, diag = T)] = modOut$V
V[lower.tri(V, diag = F)] = t(V)[lower.tri(t(V), diag = F)]
Vs = c(Vs, list(V[iEst, iEst]))
conv = c(conv, modOut$conv)
}
}
Mconv = sum(conv %in% c(1,3))
if(Mconv <= 1){
stop("No parametric boostrap iteration could converge.")
}
##compute variance
bests = bests[,conv %in% c(1,3)]
mb = apply(bests, 1, mean, na.rm = T)
Vs = Vs[conv %in% c(1,3)]
V2 = Reduce("+", Vs)/Mconv
V1 = (Mconv+1)/((Mconv-1)*Mconv) * Reduce("+", lapply(bests, function(best){
(best-mb)%*%t(best-mb)
}))
covar = V2 + V1
V = matrix(0, nOut, nOut)
V[iEst, iEst] = covar
##Create output model, based on last estimated bootstrap model structure
modOut$best[iKeepOut] = model$best[iKeepIn]
modOut$best[iEst] = mb
modOut$V = V[upper.tri(V, diag = T)]
modOut$Mconv = Mconv
##replace chol for varcov in best
if(!missing(fixed)){
modOut$cholesky[-(1:sum(nVCIn))*(nVCIn != 0)] = modOut$best[iVCOut]
if(idiag & inherits(modOut, "mpjlcmm")){ #idiag but not multlcmm because diff $chol structure
if(sum(nVCIn) == 0){
modOut$best[iVCOut] = modOut$cholesky^2
} else {
modOut$best[iVCOut] = modOut$cholesky[-(1:sum(nVCIn))]^2
}
} else if(inherits(modOut, "mpjlcmm")) { #only mpj (twoStage), because $contrainte is used to get multlcmm
isMult = as.integer(modOut$contrainte[modOut$K] == 2)
NVC = sqrt(2*(length(modOut$cholesky)+isMult-sum(nVCIn))+1/4)-1/2
cholMatrix = matrix(0, NVC, NVC)
chols = modOut$best[iVCOut]
if(modOut$contrainte[modOut$K] == 2) chols = c(1, chols) #mpjlcmm & mult : no 1 in $chol
cholMatrix[upper.tri(cholMatrix, diag = T)] = chols
vc = t(cholMatrix)%*%cholMatrix
vc = vc[upper.tri(vc, diag = T)]
if(modOut$contrainte[modOut$K] == 2) vc = vc[-1]
modOut$best[iVCOut] = vc
} else { #for conditional, but mult can be idiag. Because $chol structure is full for mult even with idiag
NVC = sqrt(2*(length(modOut$cholesky)-nVCIn)+1/4)-1/2
cholMatrix = matrix(0, NVC, NVC)
chols = modOut$best[iVCOut]
if(inherits(modOut, "multlcmm")) chols = c(1, chols)
if(idiag){
chols_idiag = c()
for(i in 1:NVC){
chols_idiag = c(chols_idiag, rep(0, i-1), chols[i])
}
chols = chols_idiag
}
cholMatrix[upper.tri(cholMatrix, diag = T)] = chols
vc = t(cholMatrix)%*%cholMatrix
if(idiag){
vc = diag(vc)
} else {
vc = vc[upper.tri(vc, diag = T)]
}
if(inherits(modOut, "multlcmm")) vc = vc[-1]
modOut$best[iVCOut] = vc
}
}
##Computations on final model
z = modOut$call
z$posfix = NULL
z$B = modOut$best
z$maxiter = 0
z$verbose = FALSE
modelEdit = eval(z)
if(method == "twoStageJoint") modOut$pprob = modelEdit$pprob
if(!missing(survival)) modOut$predSurv = modelEdit$predSurv
modOut$loglik = modelEdit$loglik
if(method == "conditional") modOut$best = modOut$best[iEst]
}
modOut$call$data = substitute(data)
##On remet longicall comme considere dans environnement present, pas dans celui de mpjlcmm
if(inherits(modOut, "mpjlcmm")){
longicall = vector("list", modOut$K)
for(k in 1:modOut$K){
longicall[[k]] = eval(modOut$call$longitudinal)[[k]]$call
}
modOut$longicall = longicall
}
cost = proc.time()-ptm
##Object Class Creation and transformation from mpjlcmm
if(!missing(survival)){
if(method == "twoStageJoint"){
best = modOut$best[iEst]
V = matrix(NA, nOut, nOut)
V[upper.tri(V, diag = T)] = modOut$V
V = V[iEst, iEst]
V = V[upper.tri(V, diag = T)]
##select output
N = modOut$N[c(2:3, 11+modOut$K+modOut$nbevt)]
Nprm = modOut$Nprm[2:(2+modOut$nbevt)]
Names = modOut$Names[c("Xnsnames", "ID", "Tnames", "TimeDepVar.name")]
Names[["Xnsnames"]] = Names[["Xnsnames"]][as.logical(modOut$idspecif)]
levels = modOut$levels[c("levelsdata", "levelssurv")]
modOut = list(nbevt = modOut$nbevt, ng = modOut$ng, ns = modOut$ns, idcom = modOut$idcom,
idspecif = modOut$idspecif, idtdv = modOut$idtdv, loglik = modOut$loglik,
best = best, V = V, gconv = modOut$gconv, conv = modOut$conv, call = cl,
niter = modOut$niter, N = N, Nprm = Nprm, pprob = pprob, Names = Names,
logspecif = modOut$logspecif, predSurv = modOut$predSurv, typrisq = modOut$typrisq,
hazardtype = modOut$hazardtype, hazardnodes = modOut$hazardnodes, nz = modOut$nz,
scoretest = modOut$scoretest, na.action = modOut$na.action, levels = levels,
AIC = 2*(length(best)-length(posfix)-modOut$loglik),
BIC = (length(best)-length(posfix))*log(modOut$ns)-2*modOut$loglik,
varest = varest, runtime = cost[3])
if(varest == "paramBoot"){
modOut$Mconv = round(Mconv/M*100, 1)
modOut$conv = as.integer(modOut$Mconv > 90)
}
}
if(method == "conditional"){
data$dummyY = 1
subLoglik = hlme(dummyY~1,
data = data,
subject = subject,
maxiter = 0,
B = c(1, 0.000001))$loglik
loglik = modOut$loglik - subLoglik
best = modOut$best[iEst]
V = matrix(NA, nOut, nOut)
V[upper.tri(V, diag = T)] = modOut$V
V = V[iEst, iEst]
V = V[upper.tri(V, diag = T)]
N = modOut$N[c(2:3, 10)]
Nprm = modOut$N[2:3]
Names = modOut$Names[c("Xnames2", "ID", "Tnames", "TimeDepVar.name")]
names(Names)[names(Names) == "Xnames2"] = "Xnsnames"
levels = modOut$levels[c("levelsdata", "levelssurv")]
modOut = list(nbevt = 1, ng = modOut$ng, ns = modOut$ns, idcom = modOut$idcom,
idspecif = modOut$idspecif, idtdv = modOut$idtdv, loglik = loglik,
best = best, V = V, gconv = modOut$gconv, conv = modOut$conv, call = cl,
niter = modOut$niter, N = N, Nprm = Nprm, pprob = pprob, Names = Names,
logspecif = modOut$logspecif, predSurv = modOut$predSurv, typrisq = modOut$hazard[[1]],
hazardtype = modOut$hazard[[2]], hazardnodes = modOut$hazard[[3]], nz = modOut$hazard[[4]],
scoretest = modOut$scoretest, na.action = modOut$na.action, levels = levels,
AIC = 2*(length(best)-length(posfix)-loglik),
BIC = (length(best)-length(posfix))*log(modOut$ns)-2*loglik,
varest = varest, runtime = cost[3])
if(varest == "paramBoot"){
modOut$Mconv = round(Mconv/M*100, 1)
modOut$conv = as.integer(modOut$Mconv > 90)
}
}
class(modOut) = c("externSurv", "externVar")
}
if(!missing(fixed)){
if(method == "twoStageJoint"){
##Get info
gconv = modOut$gconv
conv = modOut$conv
niter = modOut$niter
##With exclusively longitudinal external outcome
modUpdate = update(modOut)
modOut = modUpdate[[length(modUpdate)]]
modOut$gconv = gconv
modOut$conv = conv
modOut$niter = niter
modOut$pprob = pprob
if(inherits(modOut, "multlcmm")) modOut$N[3] = modOut$N[3]-modOut$N[1]
modOut$N[1] = 0
modOut$idprob0 = rep(0, length(modOut$idprob0))
modOut$call = cl
modOut$varest = varest
modOut$runtime = cost[3]
V = matrix(NA, length(modOut$best), length(modOut$best))
V[upper.tri(V, diag = T)] = modOut$V
V[lower.tri(V, diag = F)] = t(V)[lower.tri(t(V), diag = F)]
V = V[-c(1:(modOut$ng-1)), -c(1:(modOut$ng-1))]
V = V[upper.tri(V, diag = T)]
modOut$V = V
modOut$best = modOut$best[-c(1:nMB)]
if(varest == "paramBoot"){
modOut$Mconv = round(Mconv/M*100, 1)
modOut$conv = as.integer(modOut$Mconv > 90)
}
class(modOut) = c(class(modOut), "externVar")
}
if(method == "conditional"){
##Get info
modOut$pprob = pprob
if(inherits(modOut, "multlcmm")) modOut$N[3] = modOut$N[3]-modOut$N[1]
modOut$call = cl
modOut$varest = varest
modOut$runtime = cost[3]
if(varest == "paramBoot"){
modOut$Mconv = round(Mconv/M*100, 1)
modOut$conv = as.integer(modOut$Mconv > 90)
}
class(modOut) = c(class(modOut), "externVar")
}
}
if(!missing(classmb)){
V = matrix(NA, nOut, nOut)
V[upper.tri(V, diag = T)] = modOut$V
V[lower.tri(V, diag = F)] = t(V)[lower.tri(t(V), diag = F)]
V = V[iEst, iEst]
V = V[upper.tri(V, diag = T)]
best = modOut$best[iEst]
##Select output
N = modOut$N[1]
Names = modOut$Names[c("Xnsnames", "ID")]
levels = modOut$levels[c("levelsdata", "levelsclassmb")]
modOut = list(ng = modOut$ng, ns = modOut$ns, idprob = modOut$idprob, nv2 = modOut$nv2,
loglik = modOut$loglik, best = best, V = V, gconv = modOut$gconv,
conv = modOut$conv, call = cl, niter = modOut$niter, N = N,
pprob = modOut$pprob, Names = Names, na.action = modOut$na.action,
AIC = 2*(length(best)-length(posfix)-modOut$loglik),
BIC = (length(best)-length(posfix))*log(modOut$ns)-2*modOut$loglik,
varest = varest, method = method, runtime = cost[3])
if(varest == "paramBoot"){
modOut$Mconv = round(Mconv/M*100, 1)
modOut$conv = as.integer(modOut$Mconv > 90)
}
class(modOut) = c("externX", "externVar")
}
if(verbose){cat(paste("The externVar program took", round(cost[3],2), "seconds\n"))}
return(modOut)
}
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 externVar
Documented in externVar
#' Estimation of secondary regression models after the estimation of a primary latent class model
#'
#' This function fits regression models to relate a latent class structure (stemmed
#' from a latent class model estimated within \code{lcmm} package) with either an external
#' outcome or external class predictors.
#' Two inference techniques are implemented to account for the classification error:
#'
#' - a 2-stage estimation of the joint likelihood of the primary latent class model
#' and the secondary/ external regression;
#'
#' - a regression between the posterior latent class assignment and the external variable
#' which internally corrects for the assignment misclassification.
#'
#' It returns an object from one of the \code{lcmm} package classes.
#'
#' A. DATA STRUCTURE
#'
#' The \code{data} argument must follow specific structure for individual variables,
#' i.e. variables with a unique constant value for each subject. For an individual variable
#' given as external outcome, data value must be present only once per subject,
#' independently of any time variable used in the primary latent class.
#' For an individual variable given as external class predictor,
#' data values must be given for every row of every individual (as usual)
#'
#' B. VARIANCE ESTIMATION
#'
#' Not taking into account first stage variance with specifing \code{"none"} may lead to
#' underestimation of the final variance. When possible, Method \code{"Hessian"}
#' which relies on the combination of Hessians from the primary and secondary
#' model is recommended. However, it may become numerically intensive in the event
#' of very high number of parameters in the primary latent class model. As an
#' alternative, especially in situations with a complex primary model but rather
#' parcimonious secondary model, method \code{"paramBoot"} which implements a
#' parametric bootstrap can be used.
#'
#' @param model an object inheriting from class \code{hlme}, \code{lcmm},
#' \code{Jointlcmm}, \code{multlcmm} or \code{mpjlcmm} giving the primary latent
#' class model.
#' @param fixed optional two sided linear formula object for specifying the
#' fixed-effects in the secondary model with an external outcome variable.
#' 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.
#' @param mixture optional one-sided formula object for the class-specific fixed effects
#' in the model for the external outcome. Among the list of covariates included in fixed,
#' the covariates with class-specific regression parameters are entered in
#' 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 linear formula object for specifying the
#' random-effects on external outcome in the secondary model, if appropriate.
#' By default, no random effect is included.
#' @param subject name of the covariate representing the grouping structure.
#' Even in the absence of a hierarchical structure.
#' @param classmb optional one-sided formula specifying the external predictors of
#' latent class membership to be modelled in the secondary class-membership multinomial
#' logistic model. Covariates are separated by \code{+} on the right of the \code{~}.
#' @param survival optional two-sided formula specifying the external survival part
#' of the model.
#' @param hazard optional family of hazard function assumed for the survival model
#' (Weibull, piecewise or splines)
#' @param hazardtype optional indicator for the type of baseline risk function
#' (Specific, PH or Common)
#' @param hazardnodes optional vector containing interior nodes if \code{splines} or
#' \code{piecewise} is specified for the baseline hazard function in \code{hazard}
#' @param TimeDepVar optional vector specifying the name of the time-depending covariate
#' in the survival model
#' @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
#' @param idiag if appropriate, optional logical for the structure of the variance-covariance
#' matrix of the random-effects in the secondary model.
#' 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 if appropriate, optional logical indicating if the variance-covariance of the
#' random-effects in the secondary model 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 randomY optional logical for including an outcome-specific random intercept.
#' If FALSE no outcome-specific random intercept is added (default). If TRUE independent
#' outcome-specific random intercept with parameterized variance are included
#' @param link optional family of parameterized link functions for the external outcome
#' if appropriate. Defaults to NULL, corresponding to continuous Gaussian distribution
#' (hlme function).
#' @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 cor optional indicator for inclusion of an auto correlated Gaussian process
#' in the latent process linear (latent process) mixed model. Option "BM" indicates
#' a brownian motion with parameterized variance. Option "AR" specifies an
#' autoregressive process of order 1 with parameterized variance and correlation
#' intensity. Each option should be followed by the time variable in brackets as
#' code{cor=BM(time)}. By default, no autocorrelated Gaussian process is added.
#' @param nsim number of points to be used in the estimated link function. By default,
#' nsom=100.
#' @param range optional vector indicating the range of the outcomes (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 data Data frame containing the variables named in
#' \code{fixed}, \code{mixture}, \code{random}, \code{classmb} and \code{subject},
#' for both the current function arguments and the primary model arguments
#' Check \code{details} to get information on the data structure, especially with
#' external outcomes.
#' @param longitudinal only with \code{mpjlcmm} primary models and "twoStageJoint"
#' method: mandatory list containing the longitudinal submodels used in the primary
#' latent class model.
#' @param method character indicating the inference technique to be used:
#' \code{"twoStageJoint"} corresponds to 2-stage estimation. \code{"conditional"}
#' corresponds to the method based on the distribution of Y conditionally to the
#' true latent class membership.
#' @param varest optional character indicating the method to be used to compute the
#' variance of the regression estimates. \code{"none"} does not account for the
#' uncertainty in the primary latent class model, \code{"paramBoot"} computes the
#' total variance using a parametric bootstrap technique, \code{"Hessian"} computes
#' the total Hessian of the joint likelihood (implemented for \code{"twoStageJoint"}
#' method only). Default to \code{"Hessian"} for \code{"twoStageJoint"} method and
#' \code{"paramBoot"} for \code{"conditional"} method.
#' @param M option integer indicating the number of draws for the parametric boostrap
#' when \code{varest="paramBoot"}. Default to 200.
#' @param B optional vector of initial parameter values for the secondary model.
#' If external outcome, the vector has the same structure as a latent class model
#' estimated in the other functions of \code{lcmm} package for the same type of
#' outcome. If external class predictors (of size p), the vector is of length
#' (ng-1)*(1+p). If \code{B=NULL} (by default), internal initial values are selected.
#' @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 secondary model
#' estimation using Marquardt iterative algorithm. Defaults to 100
#' @param posfix optional vector specifying indices in parameter vector B the
#' secondary model that should not be estimated. Default to NULL, all the
#' parameters of the secondary regression 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
#' (can solve non convergence due to estimates at the boundary of the parameter space - usually 0).
#' @param verbose logical indicating whether information about computation should be
#' reported. Default to FALSE.
#' @param nproc the number cores for parallel computation. Default to 1 (sequential mode).
#' @return an object of class \code{externVar} and
#' \code{externSurv} for external survival outcomes,
#' \code{externX} for external class predictors, and
#' \code{hlme}, \code{lcmm}, or \code{multlcmm} for external longitudinal or cross-sectional outcomes.
#' @author Maris Dussartre, Cecile Proust-Lima and Viviane Philipps
#'
#'
#' @examples
#'
#' \dontrun{
#'
#'
#' ###### Estimation of the primary latent class model ######
#'
#' set.seed(1234)
#' PrimMod <- hlme(Ydep1~Time,random=~Time,subject='ID',ng=1,data=data_lcmm)
#' PrimMod2 <- hlme(Ydep1~Time,mixture=~Time,random=~Time,subject='ID',
#' ng=2,data=data_lcmm,B=random(PrimMod))
#'
#' ###### Example 1: Relationship between a latent class structure and #
#' # external class predictors ######
#'
#' # estimation of the secondary multinomial logistic model with total variance
#' # computed with the Hessian
#'
#' XextHess <- externVar(PrimMod2,
#' classmb = ~X1 + X2 + X3 + X4,
#' subject = "ID",
#' data = data_lcmm,
#' method = "twoStageJoint")
#' summary(XextHess)
#'
#' # estimation of a secondary multinomial logistic model with total variance
#' # computed with parametric Bootstrap (much longer). When using the bootstrap
#' # estimator, we recommend running first the analysis with option varest = "none"
#' # which is faster but which underestimates the variance. And then use these values
#' # as initial values when running the model with varest = "paramBoot" to obtain
#' # a valid variance of the parameters.
#'
#' XextNone <- externVar(PrimMod2,
#' classmb = ~X1 + X2 + X3 + X4,
#' subject = "ID",
#' data = data_lcmm,
#' varest = "none",
#' method = "twoStageJoint")
#'
#' XextBoot <- externVar(PrimMod2,
#' classmb = ~X1 + X2 + X3 + X4,
#' subject = "ID",
#' data = data_lcmm,
#' varest = "paramBoot",
#' method = "twoStageJoint",
#' B = XextNone$best)
#' summary(XextBoot)
#'
#'
#' ###### Example 2: Relationship between a latent class structure and #
#' # external outcome (repeatedly measured over time) ######
#'
#' # estimation of the secondary linear mixed model with total variance
#' # computed with the Hessian
#'
#' YextHess = externVar(PrimMod2, #primary model
#' fixed = Ydep2 ~ Time*X1, #secondary model
#' random = ~Time, #secondary model
#' mixture = ~Time, #secondary model
#' subject="ID",
#' data=data_lcmm,
#' method = "twoStageJoint")
#'
#'
#' # estimation of a secondary linear mixed model with total variance
#' # computed with parametric Bootstrap (much longer). When using the bootstrap
#' # estimator, we recommend running first the analysis with option varest = "none"
#' # which is faster but which underestimates the variance. And then use these values
#' # as initial values when running the model with varest = "paramBoot" to obtain
#' # a valid variance of the parameters.
#'
#' YextNone = externVar(PrimMod2, #primary model
#' fixed = Ydep2 ~ Time*X1, #secondary model
#' random = ~Time, #secondary model
#' mixture = ~Time, #secondary model
#' subject="ID",
#' data=data_lcmm,
#' varest = "none",
#' method = "twoStageJoint")
#'
#' YextBoot = externVar(PrimMod2, #primary model
#' fixed = Ydep2 ~ Time*X1, #secondary model
#' random = ~Time, #secondary model
#' mixture = ~Time, #secondary model
#' subject="ID",
#' data=data_lcmm,
#' method = "twoStageJoint",
#' B = YextNone$best,
#' varest= "paramBoot")
#'
#' summary(YextBoot)
#'
#'
#' ###### Example 3: Relationship between a latent class structure and #
#' # external outcome (survival) ######
#'
#' # estimation of the secondary survival model with total variance
#' # computed with the Hessian
#'
#' YextHess = externVar(PrimMod2, #primary model
#' survival = Surv(Tevent,Event)~ X1+mixture(X2), #secondary model
#' hazard="3-quant-splines", #secondary model
#' hazardtype="PH", #secondary model
#' subject="ID",
#' data=data_lcmm,
#' method = "twoStageJoint")
#' summary(YextHess)
#'
#'
#' # estimation of a secondary survival model with total variance
#' # computed with parametric Bootstrap (much longer). When using the bootstrap
#' # estimator, we recommend running first the analysis with option varest = "none"
#' # which is faster but which underestimates the variance. And then use these values
#' # as initial values when running the model with varest = "paramBoot" to obtain
#' # a valid variance of the parameters.
#'
#' YextNone = externVar(PrimMod2, #primary model
#' survival = Surv(Tevent,Event)~ X1+mixture(X2), #secondary model
#' hazard="3-quant-splines", #secondary model
#' hazardtype="PH", #secondary model
#' subject="ID",
#' data=data_lcmm,
#' varest = "none",
#' method = "twoStageJoint")
#'
#' YextBoot = externVar(PrimMod2, #primary model
#' survival = Surv(Tevent,Event)~ X1+mixture(X2), #secondary model
#' hazard="3-quant-splines", #secondary model
#' hazardtype="PH", #secondary model
#' subject="ID",
#' data=data_lcmm,
#' method = "twoStageJoint",
#' B = YextNone$best,
#' varest= "paramBoot")
#'
#' summary(YextBoot)
#'
#' }
#'
#'
#'
#'
#' @export
#'
#'
#'
externVar = function(model,
fixed,
mixture,
random,
subject,
classmb,
survival,
hazard = "Weibull",
hazardtype = "Specific",
hazardnodes = NULL,
TimeDepVar = NULL,
logscale = FALSE,
idiag = FALSE,
nwg = FALSE,
randomY = NULL,
link = NULL,
intnodes = NULL,
epsY = NULL,
cor = NULL,
nsim = NULL,
range = NULL,
data,
longitudinal,
method,
varest,
M = 200,
B,
convB = 0.0001,
convL = 0.0001,
convG = 0.0001,
maxiter = 100,
posfix,
partialH = FALSE,
verbose = FALSE,
nproc = 1){
ptm <- proc.time()
if(missing(model)) stop("model argument must be given")
if(!inherits(model, c("hlme", "lcmm", "multlcmm", "Jointlcmm", "mpjlcmm"))) stop('primary model class must be either "hlme", "lcmm", "multlcmm", "Jointlcmm" or "mpjlcmm"')
if(model$conv == 2) warning("primary model did not fully converge")
if(sum(c(!missing(fixed), !missing(classmb), !missing(survival))) != 1) stop("One and only one in survival, fixed or classmb must be given")
if(missing(method) | !method %in% c("twoStageJoint", "conditional")) stop('Method must be either "twoStageJoint" or "conditional"')
if(model$ng == 1) stop("Primary model does not have latent class structure (ng=1)")
if(method == "twoStageJoint" & missing(varest)) varest = "Hessian"
if(method == "conditional" & missing(varest)) varest = "paramBoot"
if(!varest %in% c("none", "paramBoot", "Hessian")) stop('Variance estimation method "varest" must be either "none", "paramBoot" or "Hessian"')
if(!is.null(link) & missing(fixed)) stop("The argument link is not to be used with external class predictor")
if(missing(posfix)) posfix = c()
cl = match.call()
##Informations about primary model
argumentsIn = as.list(model$call)
funIn = as.character(argumentsIn[[1]])
if(funIn == "jlcmm") funIn <- "Jointlcmm"
if(funIn == "mlcmm") funIn <- "multlcmm"
argumentsIn[[1]] = NULL
ng = model$ng
nIn = length(model$best)
if(is.null(argumentsIn[["classmb"]])){
oldclassmb = ~ 1
} else {
oldclassmb = formula(argumentsIn[["classmb"]])
}
if(!missing(classmb) & oldclassmb != ~1) stop("Primary model already has class predictor")
##number of MB parameters in primary model
nInMB = ncol(model.matrix(oldclassmb, data))*(ng-1)
##Get subject
if(missing(subject)){
if (model$call$subject %in% colnames(data)){
subject = model$call$subject
} else {
stop("The argument subject must be specified if different from the subject argument used in the primary model")
}
}
##Get longitudinal
if(funIn == "mpjlcmm"){
if(missing(longitudinal)) stop("The argument longitudinal is mandatory with a mpjlcmm primary model")
longCall = substitute(longitudinal)
K = length(longitudinal)
for(k in 1:K){
cl_long = as.list(longitudinal[[k]]$call)
if(inherits(longitudinal[[k]], "lcmm")) cl_long[["computeDiscrete"]] = FALSE
longitudinal[[k]]$call = as.call(cl_long)
assign(as.character(longCall[[k+1]]), longitudinal[[k]])
}
}
##nVCIn
if(funIn == "mpjlcmm"){
nVCIn = model$Nprm[3+2*model$K+(1:model$K)]
iVCIn = c()
prev = 0
for(k in 1:model$K){
iVCIn = c(iVCIn, sum(model$Nprm[c(1:3, 3:4*model$K-model$K+k-1)])+prev+1:nVCIn[k])
prev = prev+model$npmK[k]
}
} else if(funIn == "Jointlcmm"){
nVCIn = model$N[5]
iVCIn = sum(model$N[1:4]) + 1:nVCIn
} else if(funIn == "multlcmm"){
nVCIn = model$N[4]
iVCIn = sum(model$N[3]) + 1:nVCIn
} else {
nVCIn = model$N[3]
iVCIn = sum(model$N[1:2]) + 1:nVCIn
}
##pprob with new data
argumentsInEdit = argumentsIn
argumentsInEdit[["data"]] = data
argumentsInEdit[["maxiter"]] = 0
argumentsInEdit[["B"]] = model$best
argumentsInEdit[["verbose"]] = FALSE
modelEdit = do.call(funIn, argumentsInEdit)
pprob = modelEdit$pprob
arguments = list()
##Finding out the number of parameters is needed for all survival external outcome
if(!missing(survival)){
##Informations about secondary outcome model
##number of survival parameters to estimate
##number of events
surv <- cl$survival[[2]]
if(length(surv)==3) #censure droite sans troncature gauche
{
idtrunc <- 0
nom.Tevent <- as.character(surv[2])
nom.Event <- as.character(surv[3])
nom.Tentry <- NULL #si pas de troncature, Tentry=0
noms.surv <- c(nom.Tevent,nom.Event)
}
if(length(surv)==4) #censure droite et troncature
{
idtrunc <- 1
nom.Tentry <- as.character(surv[2])
nom.Tevent <- as.character(surv[3])
nom.Event <- as.character(surv[4])
noms.surv <- c(nom.Tentry,nom.Tevent,nom.Event)
}
Tevent <- getElement(object=data,name=nom.Tevent)
Event <- getElement(object=data,name=nom.Event)
nbevt <- length(attr(do.call("Surv",list(time=Tevent,event=Event,type="mstate")),"states"))
if(nbevt<1) nbevt <- 1
##get number of parameters for baseline functions
hazard <- rep(hazard,length.out=nbevt)
hazardtype <- rep(hazardtype,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"
}
hazWhat = hazard
##when not weibull, keep only the last word of hazard
hazWhat[hazard != 'Weibull'] = sapply(strsplit(hazard[hazard != 'Weibull'], "-"), `[`, 3)
hazN = rep(2, nbevt)
hazN[hazard != 'Weibull'] = sapply(strsplit(hazard[hazard != 'Weibull'], "-"), `[`, 1)
hazN = as.integer(hazN)
hazN = hazN +
(hazWhat == "Weibull")*0 +
(hazWhat == "piecewise")*(-1) +
(hazWhat == "splines")*(2)
hazN = hazN + hazN *
(hazardtype == "Specific")*(ng-1)
##we extract the number of base function parameter with constraints (ie, not PH parameters)
nSurvConstraint = hazN
hazN = hazN +
(hazardtype == "PH")*(ng-1)
##we now have all base functions parameters :
nEstY = sum(hazN)
##get number of parameters for survival covariates
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
}
}
##I extract number of variable in each formula (through model.matrix) excluding intercept
ncols = sapply(c(form.commun, form.cause, form.mixture, form.causek), function(x, data){
mm = model.matrix(x, data)
mm = mm[,-1]
if(is.null(ncol(mm))) {return(1)}
else {return(ncol(mm))}
}, data = data)
##I also need how many parameters each kind of covariate makes
nparam = c(1, nbevt, ng, nbevt*ng)
##we now have all the survival covariates parameters :
nEstX = sum(ncols*nparam)
nEst = nEstY + nEstX
}
##A model structure is needed for all longitudinal external outcome
if(!missing(fixed)){
##Manage inputs
if(missing(mixture)) mixture = ~1
if(missing(random)) random = ~-1
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(length(fixed[[2]]) != 1){
argfunctionStrMod = "multlcmm"
} else if(is.null(link)){
argfunctionStrMod = "hlme"
} else {
argfunctionStrMod = "lcmm"
}
##let's create structure for secondary model
argumentsStrMod = list()
argumentsStrMod[["fixed"]] = fixed
argumentsStrMod[["random"]] = random
argumentsStrMod[["subject"]] = subject
argumentsStrMod[["ng"]] = 1
argumentsStrMod[["idiag"]] = idiag
argumentsStrMod[["randomY"]] = randomY
if(argfunctionStrMod %in% c("lcmm", "multlcmm", "Jointlcmm")){
argumentsStrMod[["link"]] = link
argumentsStrMod[["intnodes"]] = intnodes
}
argumentsStrMod[["epsY"]] = epsY
argumentsStrMod[["cor"]] = substitute(cor)
argumentsStrMod[["nsim"]] = nsim
argumentsStrMod[["range"]] = range
argumentsStrMod[["data"]] = data
argumentsStrMod[["maxiter"]] = 0
argumentsStrMod[["verbose"]] = FALSE
strMod = do.call(argfunctionStrMod, c(argumentsStrMod))
argumentsStrMod[["mixture"]] = mixture
argumentsStrMod[["classmb"]] = ~1
argumentsStrMod[["ng"]] = ng
argumentsStrMod[["nwg"]] = nwg
argumentsStrMod[["B"]] = as.name("strMod")
strMod = do.call(argfunctionStrMod, c(argumentsStrMod))
strMod$best[1:ng+(ng-1)] = mean(strMod$best[1:ng+(ng-1)])# pourquoi?? tous les intercepts a la meme valeur. A revoir
}
##Change general model arguments
if(method == "twoStageJoint"){
##Common argument in twoStageJoint
arguments[["data"]] = data
arguments[["ng"]] = ng
arguments[["subject"]] = subject
##technical options
arguments[["maxiter"]] = maxiter
arguments[["verbose"]] = verbose
arguments[["nproc"]] = nproc
arguments[["convB"]] = convB
arguments[["convL"]] = convL
arguments[["convG"]] = convG
arguments[["partialH"]] = partialH
##primary survival
arguments[["survival"]] = argumentsIn[["survival"]]
arguments[["hazard"]] = argumentsIn[["hazard"]]
arguments[["hazardtype"]] = argumentsIn[["hazardtype"]]
arguments[["hazardnodes"]] = argumentsIn[["hazardnodes"]]
arguments[["hazardrange"]] = argumentsIn[["hazardrange"]]
arguments[["TimeDepVar"]] = argumentsIn[["TimeDepVar"]]
arguments[["logscale"]] = argumentsIn[["logscale"]]
##Primary Jointlcmm needs transformation into lcmm to be put into longitudinal.
if(funIn == "Jointlcmm"){
if(is.null(argumentsIn[["link"]])){
argfunJoint = "hlme"
} else {
argfunJoint = "lcmm"
}
argumentsJoint = argumentsIn
argumentsJoint[["survival"]] = NULL
argumentsJoint[["hazard"]] = NULL
argumentsJoint[["hazardtype"]] = NULL
argumentsJoint[["hazardnodes"]] = NULL
argumentsJoint[["hazardrange"]] = NULL
argumentsJoint[["TimeDepVar"]] = NULL
argumentsJoint[["logscale"]] = NULL
argumentsJoint[["maxiter"]] = 0
argumentsJoint[["mixture"]] = NULL
argumentsJoint[["classmb"]] = NULL
argumentsJoint[["ng"]] = 1
argumentsJoint[["nwg"]] = FALSE
argumentsJoint[["B"]] = NULL
argumentsJoint[["verbose"]] = FALSE
argumentsJoint[["data"]] = data
modNoSurv = do.call(argfunJoint, argumentsJoint)
argumentsJoint[["mixture"]] = argumentsIn[["mixture"]]
argumentsJoint[["classmb"]]= ~1
argumentsJoint[["ng"]] = argumentsIn[["ng"]]
argumentsJoint[["nwg"]] = argumentsIn[["nwg"]]
argumentsJoint[["B"]] = modNoSurv
modNoSurv = do.call(argfunJoint, argumentsJoint)
}
##Yextern survival
if(!missing(survival)){
##manage inputs
if(!is.null(argumentsIn[["survival"]])) stop('secondary survival model is not supported with "twoStageJoint" method if primary model already includes survival')
funOut = "mpjlcmm"
##nOut : nuber of total final parameters
nOut = nIn + nEst
##index
iKeepOut = c(1:nInMB, nInMB+nEst+1:(nIn-nInMB))
iKeepIn = 1:nIn
iEst = nInMB+1:nEst
##input VC, only for keep betaa (iKeepIn)
iVCKeep = iVCIn
##Id of varcov estimates
iVCOut = c()
##list of arguments
arguments[["survival"]] = survival
arguments[["hazard"]] = hazard
arguments[["hazardtype"]] = hazardtype
arguments[["hazardnodes"]] = hazardnodes
arguments[["TimeDepVar"]] = TimeDepVar
arguments[["logscale"]] = logscale
arguments[["posfix"]] = unique(c(iKeepOut, posfix))
##what is in longitudinal ?
if(funIn == "mpjlcmm"){
arguments[["longitudinal"]] = longitudinal
} else {
arguments[["longitudinal"]] = list(model)
}
##initial values
if(missing(B)){
arguments[["B"]][iEst] = rep(0.1, nEst)
} else {
if(length(B) != length(iEst)) stop("B should be of length ", length(iEst))
arguments[["B"]][iEst] = B
}
}
##Yextern longitudinal
if(!missing(fixed)){
funOut = "mpjlcmm"
##Let's change strMod's saved call
strMod$call$data = substitute(data)
## Now join the primary and secondary model
##Informations about secondary outcome model
##number of classmb parameters to remove
nMB = ng-1
##number of remaining parameters to estimate
nStr = length(strMod$best)
nEst = nStr - nMB
##nOut : nuber of total final parameters
nOut = nIn + nEst
##index
iKeepOut = 1:nIn
iKeepIn = iKeepOut
iEst = nIn+1:nEst
##input VC, only for keep betaa (iKeepIn)
iVCKeep = iVCIn
##Id of varcov estimates
if(inherits(strMod, "multlcmm")){
nVCStr = strMod$N[4]
iVCStr = sum(strMod$N[3]) + 1:nVCStr
} else {
nVCStr = strMod$N[3]
iVCStr = sum(strMod$N[1:2]) + 1:nVCStr
}
iVCOut = nIn + iVCStr - nMB
##Liste des arguments
##on fixe nos parametres
arguments[["posfix"]] = unique(c(iKeepOut, posfix))
##On donne les modeles
if(funIn == "mpjlcmm"){
arguments[["longitudinal"]] = c(longitudinal, list(strMod))
} else if(funIn == "Jointlcmm"){
arguments[["longitudinal"]] = list(modNoSurv, strMod)
} else {
arguments[["longitudinal"]] = list(model, strMod)
} ## ici : B= prm modele entre puis prm modele ajoute
##initial values
if(missing(B)){
arguments[["B"]][iEst] = strMod$best[(nMB+1):nStr]
} else {
if(length(B) != length(iEst)) stop("B should be of length ", length(iEst))
arguments[["B"]][iEst] = B
}
}
##X extern
if(!missing(classmb)){
if(!inherits(classmb,"formula")) stop("The argument classmb must be a formula")
funOut = "mpjlcmm"
##nEst : number of MB parameters in output model
nEst1G = ncol(model.matrix(classmb, data))
nEst = nEst1G*(ng-1)
##nOut : nuber of total final parameters
nOut = nIn - nInMB + nEst
##index
iKeepOut = (nEst+1):nOut
iKeepIn = (nInMB+1):nIn
iEst = 1:nEst
##input VC, only for keep betas (iKeepIn)
iVCKeep = iVCIn-nInMB
##Id of varcov estimates (none)
iVCOut = c()
##On recree tous nos arguments
if(funIn == "mpjlcmm"){
arguments[["longitudinal"]] = longitudinal
} else if(funIn == "Jointlcmm"){
arguments[["longitudinal"]] = list(modNoSurv)
} else {
arguments[["longitudinal"]] = list(model)
}
arguments[["classmb"]] = classmb
##on ajoute des valeurs de base pour nos nouveaux estimateurs
arguments[["B"]] = rep(0, nOut)
##initial values
if(missing(B)){
arguments[["B"]][1:ng-1] = model$best[1:ng-1]
} else {
if(length(B) != length(iEst)) stop("B should be of length ", length(iEst))
arguments[["B"]][iEst] = B
}
##on fixe nos parametres
arguments[["posfix"]] = unique(c(posfix, iKeepOut))
}
}
if(method == "conditional"){
##Yextern survival
if(!missing(survival)){
## we need it all in a function in order to be able to use parametric bootstrap later on
##remaining parameters to estimate
iEst = 1:nEst
iKeepIn = 1:nIn
iKeepOut = 1:nIn+nEst+2
nOut = nEst+2
##Id of varcov estimates to keep in bootstrap
iVCKeep = iVCIn
iVCOut = c()
nVCIn = 0
##We need what is inside of longitudinal to still exist in the worker
argumentsIn[["longitudinal"]] = eval(argumentsIn[["longitudinal"]])
conditionalS = function(model,
data,
survival,
hazard,
hazardtype,
hazardnodes = NULL,
TimeDepVar = NULL,
logscale,
subject,
ng,
B,
link,
iEst,
maxiter,
verbose,
argumentsIn,
funIn,
nproc,
convB,
convL,
convG){
argumentsInEdit = argumentsIn
argumentsInEdit[["B"]] = B[iKeepOut]
argumentsInEdit[["maxiter"]] = 0
argumentsInEdit[["verbose"]] = F
argumentsInEdit[["nproc"]] = nproc
argumentsInEdit[["data"]] = data
model = do.call(funIn, argumentsInEdit)
B = B[iEst]
predCl = predictClass(model, data)
##First : let's compute P(C|\tilde C) (\tilde C : A)
pAlY = sapply(1:ng, function(g){
return(as.numeric(predCl[,2] == g))
})
pClY = as.matrix(predCl[,3:(2+ng)])
if(any(is.nan(pClY))) stop("NaN in posterior classification probability")
pA = apply(pAlY, 2, mean)
pClA = t(pAlY)%*%pClY/(model$ns*pA)
if(det(pClA) == 0 | is.na(det(pClA))) stop("Computed error matrix is singular. One class might be empty")
##Then : let's add it to dataset for each individual
indivProb = pAlY%*%pClA
indivProb = cbind(predCl[,1], indivProb)
colnames(indivProb)[1] = subject
data = merge(data, indivProb, by = subject)
##We need dummy Y for the model to run
data$dummyY = 1
##With some type of input model, prob does not have the same name
if(!"prob1" %in% colnames(data)){
for(i in 1:ng){
data[[paste0("prob", i)]] = data[[paste0("probYT", i)]]
}
}
##Finally : we need to build the model for ppriors !
arguments = list()
arguments[["data"]] = data
arguments[["fixed"]] = dummyY~1
arguments[["mixture"]] = ~-1
arguments[["survival"]] = survival
arguments[["hazard"]] = hazard
arguments[["hazardtype"]] = hazardtype
arguments[["hazardnodes"]] = hazardnodes
arguments[["TimeDepVar"]] = TimeDepVar
arguments[["logscale"]] = logscale
arguments[["subject"]] = subject
arguments[["classmb"]] = ~-1
arguments[["ng"]] = ng
arguments[["pprior"]] = c("prob1", "prob2")
arguments[["posfix"]] = 1:2+length(iEst)
##technical options
arguments[["maxiter"]] = maxiter
arguments[["verbose"]] = verbose
arguments[["nproc"]] = nproc
arguments[["convB"]] = convB
arguments[["convL"]] = convL
arguments[["convG"]] = convG
arguments[["B"]] = c(B[iEst], 1, 0.000001)
res = do.call("Jointlcmm", arguments)
res$call = match.call()
return(res)
}
##we need to build the model
arguments[["data"]] = data
arguments[["survival"]] = survival
arguments[["hazard"]] = hazard
arguments[["hazardtype"]] = hazardtype
arguments[["hazardnodes"]] = hazardnodes
arguments[["TimeDepVar"]] = TimeDepVar
arguments[["logscale"]] = logscale
arguments[["model"]] = model
arguments[["subject"]] = subject
arguments[["ng"]] = ng
arguments[["link"]] = link
arguments[["iEst"]] = iEst
##technical options
arguments[["maxiter"]] = maxiter
arguments[["verbose"]] = verbose
arguments[["argumentsIn"]] = argumentsIn
arguments[["funIn"]] = funIn
arguments[["nproc"]] = nproc
arguments[["convB"]] = convB
arguments[["convL"]] = convL
arguments[["convG"]] = convG
##on ajoute des valeurs de base pour nos nouveaux estimateurs
arguments[["B"]] = rep(0.1, nIn+nOut)
##initial values
if(!missing(B)){
if(length(B) != length(iEst)) stop("B should be of length ", length(iEst))
arguments[["B"]][iEst] = B
}
funOut = "conditionalS"
}
##Yextern longitudinal
if(!missing(fixed)){
## we need it all in a function in order to be able to use parametric bootstrap later on
##number of classmb parameters to remove
nMB = ng-1
##number of remaining parameters to estimate
nStr = length(strMod$best)
nEst = nStr - nMB
iEst = 1:nEst
iKeepIn = 1:nIn
iKeepOut = 1:nIn+nEst
nOut = nEst
##Id of varcov estimates to keep in bootstrap
iVCKeep = iVCIn
##Id of varcov estimates
if(inherits(strMod, "multlcmm")){
nVCStr = strMod$N[4]
iVCStr = sum(strMod$N[3]) + 1:nVCStr
} else {
nVCStr = strMod$N[3]
iVCStr = sum(strMod$N[1:2]) + 1:nVCStr
}
iVCOut = iVCStr - nMB
nVCIn = 0
##We need what is inside of longitudinal to still exist in the worker
argumentsIn[["longitudinal"]] = eval(argumentsIn[["longitudinal"]])
conditional = function(model,
data,
fixed,
random,
idiag,
nwg,
randomY = NULL,
link,
intnodes = NULL,
epsY = NULL,
cor = NULL,
nsim = NULL,
range = NULL,
subject,
mixture,
ng,
B,
iEst,
maxiter,
verbose,
argumentsIn,
funIn,
nproc,
convB,
convL,
convG){
argumentsInEdit = argumentsIn
argumentsInEdit[["B"]] = B[iKeepOut] ## B : prm model a estimer, prm model entree
argumentsInEdit[["maxiter"]] = 0
argumentsInEdit[["verbose"]] = F
argumentsInEdit[["nproc"]] = nproc
argumentsInEdit[["data"]] = data
model = do.call(funIn, argumentsInEdit)
B = B[iEst]
##predCl = predictClass(model, data) # pas utile # chgmt Viviane
predCl <- model$pprob
##First : let's compute P(C|\tilde C) (\tilde C : A)
pAlY = sapply(1:ng, function(g){
return(as.numeric(predCl[,2] == g))
})
pClY = as.matrix(predCl[,3:(2+ng)])
if(any(is.nan(pClY))) stop("NaN in posterior classification probability")
pA = apply(pAlY, 2, mean)
pClA = t(pAlY)%*%pClY/(model$ns*pA)
if(det(pClA) == 0 | is.na(det(pClA))) stop("Computed error matrix is singular. One class might be empty")
##Then : let's add it to dataset for each individual
indivProb = pAlY%*%pClA
indivProb = cbind(predCl[,1], indivProb)
colnames(indivProb)[1] = subject
data = merge(data, indivProb, by = subject)
##Finally : we need to build the model for ppriors !
arguments = list()
if(length(fixed[[2]]) != 1){
funOut = "multlcmm"
arguments[["link"]] = link
arguments[["intnodes"]] = intnodes
} else if(missing(link)){
funOut = "hlme"
arguments[["link"]] = NULL
} else {
funOut = "lcmm"
arguments[["link"]] = link
arguments[["intnodes"]] = intnodes
}
##With some type of input model, prob does not have the same name
if(!"prob1" %in% colnames(data)){
for(i in 1:ng){
data[[paste0("prob", i)]] = data[[paste0("probYT", i)]]
}
}
arguments[["data"]] = data
arguments[["fixed"]] = fixed
arguments[["random"]] = random
arguments[["idiag"]] = idiag
arguments[["nwg"]] = nwg
arguments[["randomY"]] = randomY
arguments[["epsY"]] = epsY
arguments[["cor"]] = substitute(cor)
arguments[["nsim"]] = nsim
arguments[["range"]] = range
arguments[["subject"]] = subject
arguments[["mixture"]] = mixture
arguments[["classmb"]] = ~-1
arguments[["ng"]] = ng
arguments[["pprior"]] = c("prob1", "prob2")
##technical options
arguments[["maxiter"]] = maxiter
arguments[["verbose"]] = verbose
arguments[["nproc"]] = nproc
arguments[["convB"]] = convB
arguments[["convL"]] = convL
arguments[["convG"]] = convG
arguments[["B"]] = B ## ici B = vi du modele secondaire
res = do.call(funOut, arguments)
res$call = match.call()
return(res)
}
##we need to build the model
arguments[["data"]] = data
arguments[["model"]] = model
arguments[["fixed"]] = fixed
arguments[["random"]] = random
arguments[["idiag"]] = idiag
arguments[["nwg"]] = nwg
arguments[["randomY"]] = randomY
arguments[["link"]] = link
arguments[["intnodes"]] = intnodes
arguments[["epsY"]] = epsY
arguments[["cor"]] = substitute(cor)
arguments[["nsim"]] = nsim
arguments[["range"]] = range
arguments[["subject"]] = subject
arguments[["mixture"]] = mixture
arguments[["ng"]] = ng
arguments[["iEst"]] = iEst
##technical options
arguments[["maxiter"]] = maxiter
arguments[["verbose"]] = verbose
arguments[["argumentsIn"]] = argumentsIn
arguments[["funIn"]] = funIn
arguments[["nproc"]] = nproc
arguments[["convB"]] = convB
arguments[["convL"]] = convL
arguments[["convG"]] = convG
##on ajoute des valeurs de base pour nos nouveaux estimateurs
arguments[["B"]] = rep(0, nIn+nOut) ## ici B = valeurs pour les 2 modeles
##initial values
if(missing(B)){
arguments[["B"]][iEst] = strMod$best[(nMB+1):nStr]
} else {
if(length(B) != length(iEst)) stop("B should be of length ", length(iEst))
arguments[["B"]][iEst] = B
}
funOut = "conditional"
}
##Xextern
if(!missing(classmb)){
##we need it all in a function in order to be able to use parametric bootstrap later
##nEst : number of MB parameters in output model
nEst1G = ncol(model.matrix(classmb, data))
nEst = nEst1G*(ng-1)
iEst = 1:nEst
iKeepIn = 1:nIn
iKeepOut = 1:nIn+nEst
nOut = nEst
##Id of varcov estimates to keep in bootstrap
iVCKeep = iVCIn
##Id of varcov estimates (none)
iVCOut = c()
nVCIn = 0
##We need what is inside of longitudinal to still exist in the worker
argumentsIn[["longitudinal"]] = eval(argumentsIn[["longitudinal"]])
conditionalX = function(classmb,
data,
ng,
B,
iKeepIn,
iEst,
argumentsIn,
funIn,
nproc,
maxiter,
verbose){
argumentsInEdit = argumentsIn
argumentsInEdit[["B"]] = B[iKeepOut]
argumentsInEdit[["maxiter"]] = 0
argumentsInEdit[["verbose"]] = F
argumentsInEdit[["nproc"]] = nproc
argumentsInEdit[["data"]] = data
model = do.call(funIn, argumentsInEdit)
B = B[iEst]
predCl = predictClass(model, data)
##First : let's compute P(\tilde C|C) (\tilde C : A)
pAlY = sapply(1:ng, function(g){
return(as.numeric(predCl[,2] == g))
})
pClY = as.matrix(predCl[,3:(2+ng)])
if(any(is.nan(pClY))) stop("NaN in posterior classification probability")
betas = c(model$best[1:(ng-1)], 0)
pC = sapply(betas, function(b, betas){
exp(b)/sum(exp(betas))
}, betas=betas)
pAlC = t(pClY)%*%pAlY/(model$ns*pC)
if(det(pAlC) == 0 | is.na(det(pAlC))) stop("Computed error matrix is singular. One class might be empty")
##Then : let's add the classification to the dataset for each individual
indivProb = pAlY%*%t(pAlC)
indivProb = cbind(predCl[,1], indivProb)
colnames(indivProb) = c(subject, paste0("class", 1:ng))
data = merge(data, indivProb, by = subject)
for(id in unique(data[[subject]])){
if(sum(data[[subject]] == id) > 1) data = data[-which(data[[subject]] == id)[-1],]
}
##negative log likelihood function
nLL = function(beta, y, X) {
beta = matrix(ncol = ncol(y)-1, byrow = T, beta)
denom <- apply(exp(X%*%beta),1,sum) + 1
num_mat <- y*cbind(exp(X%*%beta), 1)
num <- apply(num_mat,1,sum)
vrais <- sum(log(num/denom))
return(-vrais)
}
##frame
argmf = list(
formula = classmb,
data = data
)
mf = do.call(model.frame, argmf)
ns = nrow(mf)
y = data[, paste0("class", 1:ng)]
y = as.matrix(sapply(y, as.numeric))
nBy = ncol(y)-1
X = model.matrix(classmb, data)
nEst = ncol(X)*nBy
opt = mla(b=B, fn=nLL, y=y, X=X, print.info = verbose, nproc = nproc, maxiter = maxiter)
namesX = c("intercept", colnames(X)[colnames(X) != "(Intercept)"])
names(opt$b) = c(sapply(namesX, FUN = function(i){
return(paste0(i, " ", colnames(y)[-ng]))
}))
Names = list(Xnsnames = namesX,
ID = subject)
N = c(nEst)
res = list(best = opt$b,
V = opt$v,
conv = opt$istop,
loglik = -opt$fn.value,
ns = length(unique(data[[subject]])),
ng = ng,
idprob = rep(1, ncol(X)),
nv2 = ncol(X),
gconv = c(opt$ca, opt$cb, opt$rdm),
pprob = NULL,
Names = Names,
N = N,
call = match.call())
return(res)
}
##Finally : we need to build the model arguments
arguments[["classmb"]] = classmb
arguments[["data"]] = data
arguments[["ng"]] = ng
arguments[["iKeepIn"]] = iKeepIn
arguments[["iEst"]] = iEst
arguments[["argumentsIn"]] = argumentsIn
arguments[["funIn"]] = funIn
arguments[["nproc"]] = nproc
arguments[["maxiter"]] = maxiter
arguments[["verbose"]] = verbose
##on ajoute des valeurs de base pour nos nouveaux estimateurs
arguments[["B"]] = rep(0, nIn+nOut)
##initial values
if(!missing(B)){
if(length(B) != length(iEst)) stop("B should be of length ", length(iEst))
arguments[["B"]][iEst] = B
}
funOut = "conditionalX"
}
}
if(varest != "paramBoot"){
arguments[["B"]][iKeepOut] = model$best[iKeepIn]
if(verbose){cat("Model estimation...\n\n")}
##Model Estimation
modOut = do.call(funOut, c(arguments))
}
if(varest == "Hessian"){
if(method != "twoStageJoint") stop("Hessian variance estimation method only avaliable for 'twoStageJoint' method")
if(verbose){cat("Variance estimation...\n\n")}
nb11 = length(model$best)
V11 = matrix(0, nb11, nb11)
V11[upper.tri(V11, diag=T)] = model$V
V11[lower.tri(V11, diag=F)] = t(V11)[lower.tri(V11, diag=F)]
n1 = model$ns
V11 = V11[iKeepIn, iKeepIn]
nb22 = length(modOut$best)
V22 = matrix(0, nb22, nb22)
V22[upper.tri(V22, diag=T)] = modOut$V
V22[lower.tri(V22, diag=F)] = t(V22)[lower.tri(V22, diag=F)]
saveV22 = V22
V22 = V22[iEst, iEst]
n2 = modOut$ns
modOut
if(nproc == 1){
I12 = -hessienne(modOut)
} else {
I12 = -hessienne(modOut, method = "deriva", nproc = nproc)
}
I12 = I12[iEst, iKeepOut]
V = V22*n2 + (V22*n2) %*% (I12/n2) %*% ((n2/n1)*(V11*n1)) %*% t(I12/n2) %*% (V22*n2)
V = V/n2
saveV22[iEst, iEst] = V
V = saveV22
modOut$V = V[upper.tri(V, diag = TRUE)]
}
##Get Bootstrap Models
if(varest == "paramBoot"){
if(verbose){cat("Bootstrap estimation...\n\n")}
est = estimates(model)
Vin = matrix(0, length(est), length(est))
Vin[upper.tri(Vin, diag = T)] = model$V
Vin[lower.tri(Vin, diag = F)] = t(Vin)[lower.tri(t(Vin), diag=F)]
est = est[iKeepIn]
Vin = Vin[iKeepIn,iKeepIn]
coefss = rmvnorm(M, est, Vin)
coefss = as.data.frame(coefss)
colnames(coefss) = names(model$best)[iKeepIn]
##we just need to build back varcov into the coefs instead of cholesky matrix
coefss = apply(coefss, 1, function(coefs, model, data, iVCKeep){
if(funIn == "mpjlcmm"){
varcovMods = longitudinal
} else {
varcovMods = list(model)
}
chols = coefs[iVCKeep]
model$cholesky = chols
varcov = c()
countChol = 0
for(varcovMod in varcovMods){
ncolRandMod = ncol(model.matrix(formula(varcovMod$call$random), data))
if(varcovMod$idiag){
nChol = ncolRandMod-as.integer(funIn == "multlcmm")
vc = chols[1:nChol+countChol]
varcov = c(varcov, (vc^2))
} else {
nChol = ncolRandMod*(ncolRandMod+1)/2-as.integer(funIn == "multlcmm")
##cholMatrix
cholMatrix = matrix(0, ncolRandMod, ncolRandMod)
cholsToMatrix = chols[1:nChol+countChol]
if(funIn == "multlcmm") cholsToMatrix = c(1, cholsToMatrix)
cholMatrix[upper.tri(cholMatrix, diag = T)] = cholsToMatrix
vc = t(cholMatrix)%*%cholMatrix
varcov = c(varcov, vc[upper.tri(vc, diag = T)])
if(funIn == "multlcmm") varcov = varcov[-1]
}
countChol = countChol + nChol
}
coefs[iVCKeep] = varcov
return(coefs)
}, model = model, data = data, iVCKeep = iVCKeep)
if(nproc > 1)
{
clust <- parallel::makeCluster(nproc)
##load all loaded packages
packages = loadedNamespaces()
for(pack in packages){
clusterExport(clust, "pack", environment())
clusterEvalQ(clust, require(pack, character.only = T))
}
survivalMissing = missing(survival)
fixedMissing = missing(fixed)
modOuts <- parApply(clust, coefss, 2, function(coefs, arguments, iKeepOut, funOut, iEst, survivalMissing, fixedMissing, logscale){
arguments[["B"]][iKeepOut] = coefs
arguments[["nproc"]] = 1
##Model Estimation
modOut = do.call(funOut, c(arguments))
######## chgmt Viviane ########
## if(!fixedMissing){
## ##cholesky not varcov as output in best
## modOut$best = estimates(modOut)
## ## abs value of link prm :
## n = length(modOut$best)
## if(inherits(modOut, "lcmm")){ #For lcmm
## nLink <- n - sum(modOut$N[c(1:4,6)])
## if(modOut$linktype == 2){ #with spline link
## ##nSpl = n-sum(modOut$N[1:4]) #count number of link function parameters (the only one not in N)
## ##modOut$best[n-2:nSpl+1+1] = abs(modOut$best[n-2:nSpl+1+1]) #all but the first
## modOut$best[sum(modOut$N[1:4])+2:nLink] = abs(modOut$best[sum(modOut$N[1:4])+2:nLink]) #all but the first
## } else { #rest of lcmm
## modOut$best[sum(modOut$N[1:4])+nLink] = abs(modOut$best[sum(modOut$N[1:4])+nLink]) #only the last one
## }
## } else if (inherits(modOut, "hlme")){ #hlme
## modOut$best[n] = abs(modOut$best[n]) #only the last one
## } else if(inherits(modOut, "multlcmm")){
## nPreLink = sum(modOut$N[3:8]) #number of parameters before the one for the link function
## numSPL = 0
## for (ny in 1:modOut$N[8]){
## if (modOut$linktype[ny]==0) nLink = 2
## if (modOut$linktype[ny]==1) nLink = 4
## if (modOut$linktype[ny]==2){
## numSPL <- numSPL+1
## nLink = modOut$nbnodes[numSPL]+2
## modOut$best[nPreLink+2:nLink] = abs(modOut$best[nPreLink+2:nLink])
## } else {
## modOut$best[nPreLink+nLink] = abs(modOut$best[nPreLink+nLink])
## }
## nPreLink = nPreLink+nLink
## }
## } else if(inherits(modOut, "mpjlcmm")){
## #Residual Error : need to be the same sign across bootstrap iterations
## nPre = sum(modOut$N[1:(2+modOut$nbevt)])
## sumny = 0
## for(k in 1:modOut$K){
## if(modOut$contrainte[k] == 2){ ## attention : length(nrisq) = nbevt!!
## nPre = nPre + sum(modOut$Nprm[2+modOut$nbevt+1:7*modOut$K-modOut$K+k]) #add number of parameters before the one for the link function
## } else {
## nPre = nPre + sum(modOut$Nprm[2+modOut$nbevt+1:8*modOut$K-modOut$K+k]) #add the number of parameters for this K
## }
## for(y in 1:modOut$ny[k]){
## sumny = sumny+1
## if(modOut$contrainte[k] == 1 & modOut$linktype[sumny] == 2){ #for lcmm with spl link
## nSpl = sum(modOut$Nprm[2+modOut$nbevt+7*modOut$K+sumny]) #count number of link function parameters (the only one not in N)
## modOut$best[nPre-2:nSpl+1+1] = abs(modOut$best[nPre-2:nSpl+1+1]) #all but the first
## } else if (modOut$contrainte[k] == 0 | modOut$contrainte[k] == 1){ #hlme ou le reste de lcmm
## modOut$best[nPre] = abs(modOut$best[nPre]) #only the last one
## } else if (modOut$contrainte[k] == 2) {
## numSPL = 0
## if (modOut$linktype[sumny]==0) nLink = 2
## if (modOut$linktype[sumny]==1) nLink = 4
## if (modOut$linktype[sumny]==2){
## numSPL <- numSPL+1
## nLink = modOut$call$longitudinal[[k]]$nbnodes[numSPL]+2
## modOut$best[nPre+2:nLink] = abs(modOut$best[nPre+2:nLink])
## } else {
## modOut$best[nPre+nLink] = abs(modOut$best[nPre+nLink])
## }
## nPre = nPre+nLink
## }
## }
## }
## }
## }
## ##nw & cor & alea
## idnw <- NULL
## idalea <- NULL
## idcor <- NULL
## if(inherits(modOut, "multlcmm")){#with multlcmm
## if(modOut$N[5] != 0) idnw = sum(modOut$N[3:4])+1:modOut$N[5]
## if(modOut$N[6] != 0) idalea = sum(modOut$N[3:5])+1:modOut$N[6]
## if(modOut$N[7] != 0) idcor = sum(modOut$N[3:6])+modOut$N[7]
## } else if(inherits(modOut, "mpjlcmm")){#with mpjlcmm
## nPre = sum(modOut$N[1:3])
## idnw = c()
## idcor = c()
## idalea = c()
## for(k in 1:modOut$K){
## nprenw = nPre+sum(modOut$Nprm[c(3+1:3*modOut$K-modOut$K+k)])
## if(modOut$Nprm[3+4*modOut$K-modOut$K+k] != 0){
## idnw = c(idnw, nprenw+1:modOut$Nprm[3+4*modOut$K-modOut$K+k])
## }
## nprecor = nPre+sum(modOut$Nprm[c(3+1:4*modOut$K-modOut$K+k)])
## if(modOut$Nprm[3+5*modOut$K-modOut$K+k] != 0){
## idcor = c(idcor, nprecor+modOut$Nprm[3+5*modOut$K-modOut$K+k])
## }
## nprealea = nPre+sum(modOut$Nprm[c(3+1:6*modOut$K-modOut$K+k)])
## if(modOut$Nprm[3+7*modOut$K-modOut$K+k] != 0){
## idalea = c(idalea, nprealea+1:modOut$Nprm[3+7*modOut$K-modOut$K+k])
## }
## nPre = nPre+sum(modOut$Nprm[c(3+1:8*modOut$K-modOut$K+k)])
## }
## } else {#with lcmm & hlme
## if(modOut$N[4] != 0) idnw = sum(modOut$N[1:3])+1:modOut$N[4]
## if(modOut$N[5] != 0) idcor = sum(modOut$N[1:4])+modOut$N[5]
## }
## modOut$best[idnw] = abs(modOut$best[idnw])
## modOut$best[idcor] = abs(modOut$best[idcor])
## modOut$best[idalea] = abs(modOut$best[idalea])
## #Survival Base Function : need to be the same sign across bootstrap iterations
## if(!survivalMissing & !logscale){
## iSurvConstraint = 1:nSurvConstraint+modOut$N[1]
## modOut$best[iSurvConstraint] = abs(modOut$best[iSurvConstraint])
## }
######## fin chgmt Viviane ########
if(inherits(modOut, c("hlme", "lcmm", "multlcmm", "Jointlcmm", "mpjlcmm"))) modOut <- absprm(modOut) ## remplace ce que j'ai mis en commentaire
## si model multinomial, pas besoin de valeur absolue dans les prm
return(modOut)
}, arguments = arguments, iKeepOut = iKeepOut, funOut = funOut, iEst = iEst, survivalMissing = survivalMissing, fixedMissing = fixedMissing, logscale = logscale)
parallel::stopCluster(clust)
##format output
bests = as.data.frame(matrix(NA, nrow = nEst, ncol = M))
Vs = list()
conv = c()
modOut = modOuts[[M]]
for (i in 1:M){
##output V and betas
bests[,i] = modOut$best[iEst]
V = matrix(NA, nOut, nOut)
V[upper.tri(V, diag = T)] = modOut$V
V[lower.tri(V, diag = F)] = t(V)[lower.tri(t(V), diag = F)]
Vs = c(Vs, list(V[iEst, iEst]))
conv = c(conv, modOut$conv)
}
} else {
##estimate final models, extract usefull information
bests = as.data.frame(matrix(NA, nrow = nEst, ncol = M))
Vs = list()
conv = c()
for(i in 1:M){
if(verbose) cat("==================== Bootstrap Iteration", i, "====================\n")
arguments[["B"]][iKeepOut] = coefss[,i]
##Model Estimation
if(verbose){
modOut = do.call(funOut, c(arguments))
} else {
captured_log = capture.output({modOut = do.call(funOut, c(arguments))})
}
######## chgmt Viviane ########
## if(!missing(fixed)){
## #cholesky not varcov as output in best
## modOut$best = estimates(modOut)
## n = length(modOut$best)
## if(inherits(modOut, "lcmm")){ #For lcmm
## if(modOut$linktype == 2){ #with spline link
## nSpl = n-sum(modOut$N[1:4]) #count number of link function parameters (the only one not in N)
## modOut$best[n-2:nSpl+1+1] = abs(modOut$best[n-2:nSpl+1+1]) #all but the first
## } else { #rest of lcmm
## modOut$best[n] = abs(modOut$best[n]) #only the last one
## }
## } else if (inherits(modOut, "hlme")){ #hlme
## modOut$best[n] = abs(modOut$best[n]) #only the last one
## } else if(inherits(modOut, "multlcmm")){
## nPreLink = modOut$N[1:5] #number of parameters before the one for the link function
## numSPL = 0
## for (ny in 1:modOut$N[8]){
## if (modOut$linktype[ny]==0) nLink = 2
## if (modOut$linktype[ny]==1) nLink = 4
## if (modOut$linktype[ny]==2){
## numSPL <- numSPL+1
## nLink = modOut$nbnodes[numSPL]+2
## modOut$best[nPreLink+2:nLink] = abs(modOut$best[nPreLink+2:nLink])
## } else {
## modOut$best[nPreLink+nLink] = abs(modOut$best[nPreLink+nLink])
## }
## nPreLink = nPreLink+nLink
## }
## } else if(inherits(modOut, "mpjlcmm")){
## #Residual Error : need to be the same sign across bootstrap iterations
## nPre = sum(modOut$N[1:3])
## sumny = 0
## for(k in 1:modOut$K){
## if(modOut$contrainte[k] == 2){
## nPre = nPre + sum(modOut$Nprm[3+1:7*modOut$K-modOut$K+k]) #add number of parameters before the one for the link function
## } else {
## nPre = nPre + sum(modOut$Nprm[3+1:8*modOut$K-modOut$K+k]) #add the number of parameters for this K
## }
## for(y in 1:modOut$ny[k]){
## sumny = sumny+1
## if(modOut$contrainte[k] == 1 & modOut$linktype[sumny] == 2){ #for lcmm with spl link
## nSpl = sum(modOut$Nprm[3+7*modOut$K+1*modOut$K-modOut$K+k]) #count number of link function parameters (the only one not in N)
## modOut$best[nPre-2:nSpl+1+1] = abs(modOut$best[nPre-2:nSpl+1+1]) #all but the first
## } else if (modOut$contrainte[k] == 0 | modOut$contrainte[k] == 1){ #hlme ou le reste de lcmm
## modOut$best[nPre] = abs(modOut$best[nPre]) #only the last one
## } else if (modOut$contrainte[k] == 2) {
## numSPL = 0
## if (modOut$linktype[sumny]==0) nLink = 2
## if (modOut$linktype[sumny]==1) nLink = 4
## if (modOut$linktype[sumny]==2){
## numSPL <- numSPL+1
## nLink = modOut$call$longitudinal[[k]]$nbnodes[numSPL]+2
## modOut$best[nPre+2:nLink] = abs(modOut$best[nPre+2:nLink])
## } else {
## modOut$best[nPre+nLink] = abs(modOut$best[nPre+nLink])
## }
## nPre = nPre+nLink
## }
## }
## }
## }
## }
## ##nw & cor & alea
## idnw <- NULL
## idalea <- NULL
## idcor <- NULL
## if(inherits(modOut, "multlcmm")){#with multlcmm
## if(modOut$N[5] != 0) idnw = sum(modOut$N[3:4])+1:modOut$N[5]
## if(modOut$N[6] != 0) idalea = sum(modOut$N[3:5])+1:modOut$N[6]
## if(modOut$N[7] != 0) idcor = sum(modOut$N[3:6])+modOut$N[7]
## } else if(inherits(modOut, "mpjlcmm")){#with mpjlcmm
## nPre = sum(modOut$N[1:3])
## idnw = c()
## idcor = c()
## idalea = c()
## for(k in 1:modOut$K){
## nprenw = nPre+sum(modOut$Nprm[c(3+1:3*modOut$K-modOut$K+k)])
## if(modOut$Nprm[3+4*modOut$K-modOut$K+k] != 0){
## idnw = c(idnw, nprenw+1:modOut$Nprm[3+4*modOut$K-modOut$K+k])
## }
## nprecor = nPre+sum(modOut$Nprm[c(3+1:4*modOut$K-modOut$K+k)])
## if(modOut$Nprm[3+5*modOut$K-modOut$K+k] != 0){
## idcor = c(idcor, nprecor+modOut$Nprm[3+5*modOut$K-modOut$K+k])
## }
## nprealea = nPre+sum(modOut$Nprm[c(3+1:6*modOut$K-modOut$K+k)])
## if(modOut$Nprm[3+7*modOut$K-modOut$K+k] != 0){
## idalea = c(idalea, nprealea+1:modOut$Nprm[3+7*modOut$K-modOut$K+k])
## }
## nPre = nPre+sum(modOut$Nprm[c(3+1:8*modOut$K-modOut$K+k)])
## }
## } else {#with lcmm
## if(inherits(modOut, "lcmm")){
## if(modOut$N[4] != 0) idnw = sum(modOut$N[1:3])+1:modOut$N[4]
## if(modOut$N[6] != 0) idcor = sum(modOut$N[1:4])+ntr+modOut$N[5]
## }
## else { #hlme
## if(inherits(modOut, "hlme")){
## if(modOut$N[4] != 0) idnw = sum(modOut$N[1:3])+1:modOut$N[4]
## if(modOut$N[5] != 0) idcor = sum(modOut$N[1:4])+modOut$N[5]
## }
## }
## }
## modOut$best[idnw] = abs(modOut$best[idnw])
## modOut$best[idcor] = abs(modOut$best[idcor])
## modOut$best[idalea] = abs(modOut$best[idalea])
## #Survival Base Function constraint : need to be the same sign across bootstrap iterations
## if(!missing(survival) & !logscale){
## iSurvConstraint = 1:nSurvConstraint+modOut$N[1]
## modOut$best[iSurvConstraint] = abs(modOut$best[iSurvConstraint])
## }
######## chgmt Viviane ########
if(inherits(modOut, c("hlme", "lcmm", "multlcmm", "Jointlcmm", "mpjlcmm"))) modOut <- absprm(modOut) ## remplace ce que j'ai mis en commentaire
## si model multinomial, pas besoin de valeur absolue dans les prm
##output V and betas
bests[,i] = modOut$best[iEst]
V = matrix(NA, nOut, nOut)
V[upper.tri(V, diag = T)] = modOut$V
V[lower.tri(V, diag = F)] = t(V)[lower.tri(t(V), diag = F)]
Vs = c(Vs, list(V[iEst, iEst]))
conv = c(conv, modOut$conv)
}
}
Mconv = sum(conv %in% c(1,3))
if(Mconv <= 1){
stop("No parametric boostrap iteration could converge.")
}
##compute variance
bests = bests[,conv %in% c(1,3)]
mb = apply(bests, 1, mean, na.rm = T)
Vs = Vs[conv %in% c(1,3)]
V2 = Reduce("+", Vs)/Mconv
V1 = (Mconv+1)/((Mconv-1)*Mconv) * Reduce("+", lapply(bests, function(best){
(best-mb)%*%t(best-mb)
}))
covar = V2 + V1
V = matrix(0, nOut, nOut)
V[iEst, iEst] = covar
##Create output model, based on last estimated bootstrap model structure
modOut$best[iKeepOut] = model$best[iKeepIn]
modOut$best[iEst] = mb
modOut$V = V[upper.tri(V, diag = T)]
modOut$Mconv = Mconv
##replace chol for varcov in best
if(!missing(fixed)){
modOut$cholesky[-(1:sum(nVCIn))*(nVCIn != 0)] = modOut$best[iVCOut]
if(idiag & inherits(modOut, "mpjlcmm")){ #idiag but not multlcmm because diff $chol structure
if(sum(nVCIn) == 0){
modOut$best[iVCOut] = modOut$cholesky^2
} else {
modOut$best[iVCOut] = modOut$cholesky[-(1:sum(nVCIn))]^2
}
} else if(inherits(modOut, "mpjlcmm")) { #only mpj (twoStage), because $contrainte is used to get multlcmm
isMult = as.integer(modOut$contrainte[modOut$K] == 2)
NVC = sqrt(2*(length(modOut$cholesky)+isMult-sum(nVCIn))+1/4)-1/2
cholMatrix = matrix(0, NVC, NVC)
chols = modOut$best[iVCOut]
if(modOut$contrainte[modOut$K] == 2) chols = c(1, chols) #mpjlcmm & mult : no 1 in $chol
cholMatrix[upper.tri(cholMatrix, diag = T)] = chols
vc = t(cholMatrix)%*%cholMatrix
vc = vc[upper.tri(vc, diag = T)]
if(modOut$contrainte[modOut$K] == 2) vc = vc[-1]
modOut$best[iVCOut] = vc
} else { #for conditional, but mult can be idiag. Because $chol structure is full for mult even with idiag
NVC = sqrt(2*(length(modOut$cholesky)-nVCIn)+1/4)-1/2
cholMatrix = matrix(0, NVC, NVC)
chols = modOut$best[iVCOut]
if(inherits(modOut, "multlcmm")) chols = c(1, chols)
if(idiag){
chols_idiag = c()
for(i in 1:NVC){
chols_idiag = c(chols_idiag, rep(0, i-1), chols[i])
}
chols = chols_idiag
}
cholMatrix[upper.tri(cholMatrix, diag = T)] = chols
vc = t(cholMatrix)%*%cholMatrix
if(idiag){
vc = diag(vc)
} else {
vc = vc[upper.tri(vc, diag = T)]
}
if(inherits(modOut, "multlcmm")) vc = vc[-1]
modOut$best[iVCOut] = vc
}
}
##Computations on final model
z = modOut$call
z$posfix = NULL
z$B = modOut$best
z$maxiter = 0
z$verbose = FALSE
modelEdit = eval(z)
if(method == "twoStageJoint") modOut$pprob = modelEdit$pprob
if(!missing(survival)) modOut$predSurv = modelEdit$predSurv
modOut$loglik = modelEdit$loglik
if(method == "conditional") modOut$best = modOut$best[iEst]
}
modOut$call$data = substitute(data)
##On remet longicall comme considere dans environnement present, pas dans celui de mpjlcmm
if(inherits(modOut, "mpjlcmm")){
longicall = vector("list", modOut$K)
for(k in 1:modOut$K){
longicall[[k]] = eval(modOut$call$longitudinal)[[k]]$call
}
modOut$longicall = longicall
}
cost = proc.time()-ptm
##Object Class Creation and transformation from mpjlcmm
if(!missing(survival)){
if(method == "twoStageJoint"){
best = modOut$best[iEst]
V = matrix(NA, nOut, nOut)
V[upper.tri(V, diag = T)] = modOut$V
V = V[iEst, iEst]
V = V[upper.tri(V, diag = T)]
##select output
N = modOut$N[c(2:3, 11+modOut$K+modOut$nbevt)]
Nprm = modOut$Nprm[2:(2+modOut$nbevt)]
Names = modOut$Names[c("Xnsnames", "ID", "Tnames", "TimeDepVar.name")]
Names[["Xnsnames"]] = Names[["Xnsnames"]][as.logical(modOut$idspecif)]
levels = modOut$levels[c("levelsdata", "levelssurv")]
modOut = list(nbevt = modOut$nbevt, ng = modOut$ng, ns = modOut$ns, idcom = modOut$idcom,
idspecif = modOut$idspecif, idtdv = modOut$idtdv, loglik = modOut$loglik,
best = best, V = V, gconv = modOut$gconv, conv = modOut$conv, call = cl,
niter = modOut$niter, N = N, Nprm = Nprm, pprob = pprob, Names = Names,
logspecif = modOut$logspecif, predSurv = modOut$predSurv, typrisq = modOut$typrisq,
hazardtype = modOut$hazardtype, hazardnodes = modOut$hazardnodes, nz = modOut$nz,
scoretest = modOut$scoretest, na.action = modOut$na.action, levels = levels,
AIC = 2*(length(best)-length(posfix)-modOut$loglik),
BIC = (length(best)-length(posfix))*log(modOut$ns)-2*modOut$loglik,
varest = varest, runtime = cost[3])
if(varest == "paramBoot"){
modOut$Mconv = round(Mconv/M*100, 1)
modOut$conv = as.integer(modOut$Mconv > 90)
}
}
if(method == "conditional"){
data$dummyY = 1
subLoglik = hlme(dummyY~1,
data = data,
subject = subject,
maxiter = 0,
B = c(1, 0.000001))$loglik
loglik = modOut$loglik - subLoglik
best = modOut$best[iEst]
V = matrix(NA, nOut, nOut)
V[upper.tri(V, diag = T)] = modOut$V
V = V[iEst, iEst]
V = V[upper.tri(V, diag = T)]
N = modOut$N[c(2:3, 10)]
Nprm = modOut$N[2:3]
Names = modOut$Names[c("Xnames2", "ID", "Tnames", "TimeDepVar.name")]
names(Names)[names(Names) == "Xnames2"] = "Xnsnames"
levels = modOut$levels[c("levelsdata", "levelssurv")]
modOut = list(nbevt = 1, ng = modOut$ng, ns = modOut$ns, idcom = modOut$idcom,
idspecif = modOut$idspecif, idtdv = modOut$idtdv, loglik = loglik,
best = best, V = V, gconv = modOut$gconv, conv = modOut$conv, call = cl,
niter = modOut$niter, N = N, Nprm = Nprm, pprob = pprob, Names = Names,
logspecif = modOut$logspecif, predSurv = modOut$predSurv, typrisq = modOut$hazard[[1]],
hazardtype = modOut$hazard[[2]], hazardnodes = modOut$hazard[[3]], nz = modOut$hazard[[4]],
scoretest = modOut$scoretest, na.action = modOut$na.action, levels = levels,
AIC = 2*(length(best)-length(posfix)-loglik),
BIC = (length(best)-length(posfix))*log(modOut$ns)-2*loglik,
varest = varest, runtime = cost[3])
if(varest == "paramBoot"){
modOut$Mconv = round(Mconv/M*100, 1)
modOut$conv = as.integer(modOut$Mconv > 90)
}
}
class(modOut) = c("externSurv", "externVar")
}
if(!missing(fixed)){
if(method == "twoStageJoint"){
##Get info
gconv = modOut$gconv
conv = modOut$conv
niter = modOut$niter
##With exclusively longitudinal external outcome
modUpdate = update(modOut)
modOut = modUpdate[[length(modUpdate)]]
modOut$gconv = gconv
modOut$conv = conv
modOut$niter = niter
modOut$pprob = pprob
if(inherits(modOut, "multlcmm")) modOut$N[3] = modOut$N[3]-modOut$N[1]
modOut$N[1] = 0
modOut$idprob0 = rep(0, length(modOut$idprob0))
modOut$call = cl
modOut$varest = varest
modOut$runtime = cost[3]
V = matrix(NA, length(modOut$best), length(modOut$best))
V[upper.tri(V, diag = T)] = modOut$V
V[lower.tri(V, diag = F)] = t(V)[lower.tri(t(V), diag = F)]
V = V[-c(1:(modOut$ng-1)), -c(1:(modOut$ng-1))]
V = V[upper.tri(V, diag = T)]
modOut$V = V
modOut$best = modOut$best[-c(1:nMB)]
if(varest == "paramBoot"){
modOut$Mconv = round(Mconv/M*100, 1)
modOut$conv = as.integer(modOut$Mconv > 90)
}
class(modOut) = c(class(modOut), "externVar")
}
if(method == "conditional"){
##Get info
modOut$pprob = pprob
if(inherits(modOut, "multlcmm")) modOut$N[3] = modOut$N[3]-modOut$N[1]
modOut$call = cl
modOut$varest = varest
modOut$runtime = cost[3]
if(varest == "paramBoot"){
modOut$Mconv = round(Mconv/M*100, 1)
modOut$conv = as.integer(modOut$Mconv > 90)
}
class(modOut) = c(class(modOut), "externVar")
}
}
if(!missing(classmb)){
V = matrix(NA, nOut, nOut)
V[upper.tri(V, diag = T)] = modOut$V
V[lower.tri(V, diag = F)] = t(V)[lower.tri(t(V), diag = F)]
V = V[iEst, iEst]
V = V[upper.tri(V, diag = T)]
best = modOut$best[iEst]
##Select output
N = modOut$N[1]
Names = modOut$Names[c("Xnsnames", "ID")]
levels = modOut$levels[c("levelsdata", "levelsclassmb")]
modOut = list(ng = modOut$ng, ns = modOut$ns, idprob = modOut$idprob, nv2 = modOut$nv2,
loglik = modOut$loglik, best = best, V = V, gconv = modOut$gconv,
conv = modOut$conv, call = cl, niter = modOut$niter, N = N,
pprob = modOut$pprob, Names = Names, na.action = modOut$na.action,
AIC = 2*(length(best)-length(posfix)-modOut$loglik),
BIC = (length(best)-length(posfix))*log(modOut$ns)-2*modOut$loglik,
varest = varest, method = method, runtime = cost[3])
if(varest == "paramBoot"){
modOut$Mconv = round(Mconv/M*100, 1)
modOut$conv = as.integer(modOut$Mconv > 90)
}
class(modOut) = c("externX", "externVar")
}
if(verbose){cat(paste("The externVar program took", round(cost[3],2), "seconds\n"))}
return(modOut)
}
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.