Nothing
R/hlme.R
In lcmm: Extended Mixed Models Using Latent Classes and Latent Processes
Defines functions
hlme
Documented in
hlme
#' Estimation of latent class linear mixed models
#'
#' This function fits linear mixed models and latent class linear mixed models
#' (LCLMM) also known as growth mixture models or heterogeneous linear mixed
#' models. The LCLMM consists in assuming that the population is divided in a
#' finite number of latent classes. Each latent class is characterised by a
#' specific trajectory modelled by a class-specific linear mixed model. Both
#' the latent class membership and the trajectory can be explained according to
#' covariates. This function is limited to a mixture of Gaussian outcomes. For
#' other types of outcomes, please see function \code{lcmm}. For multivariate
#' longitudinal outcomes, please see \code{multlcmm}.
#'
#'
#' A. THE VECTOR OF PARAMETERS B
#'
#' The parameters in the vector of initial values \code{B} or equivalently in
#' the vector of maximum likelihood estimates \code{best} are included in the
#' following order:
#'
#' (1) ng-1 parameters are required for intercepts in the latent class
#' membership model, and when covariates are included in \code{classmb}, ng-1
#' paramaters should be entered for each covariate;
#'
#' (2) for all covariates in \code{fixed}, one parameter is required if the
#' covariate is not in \code{mixture}, ng paramaters are required if the
#' covariate is also in \code{mixture};
#'
#' (3) the variance of each random-effect specified in \code{random} (including
#' the intercept) when \code{idiag=TRUE}, or the inferior triangular
#' variance-covariance matrix of all the random-effects when
#' \code{idiag=FALSE};
#'
#' (4) only when \code{nwg=TRUE}, ng-1 parameters are required for the ng-1
#' class-specific proportional coefficients in the variance covariance matrix
#' of the random-effects;
#'
#' (5) when \code{cor} is specified, 1 parameter corresponding to the variance
#' of the Brownian motion should be entered with \code{cor=BM} and 2 parameters
#' corresponding to the correlation and the variance parameters of the
#' autoregressive process should be entered
#'
#' (6) the standard error of the residual error.
#'
#' B. CAUTIONS
#'
#' Some caution should be made when using the program:
#'
#' (1) As the log-likelihood of a latent class model can have multiple maxima,
#' a careful choice of the initial values is crucial for ensuring convergence
#' toward the global maximum. The program can be run without entering the
#' vector of initial values (see point 2). However, we recommend to
#' systematically enter initial values in \code{B} and try different sets of
#' initial values.
#'
#' (2) The automatic choice of initial values we provide requires the
#' estimation of a preliminary linear mixed model. The user should be aware
#' that first, this preliminary analysis can take time for large datatsets and
#' second, that the generated initial values can be very not likely and even
#' may converge slowly to a local maximum. This is the reason why several
#' alternatives exist. The vector of initial values can be directly specified
#' in \code{B} the initial values can be generated (automatically or randomly)
#' from a model with \code{ng=}. Finally, function \code{gridsearch} performs
#' an automatic grid search.
#'
#' (3) Convergence criteria are very strict as they are based on the
#' derivatives of the log-likelihood in addition to the parameter stability and
#' log-likelihood stability. In some cases, the program may not converge and
#' reach the maximum number of iterations fixed at 100. In this case, the user
#' should check that parameter estimates at the last iteration are not on the
#' boundaries of the parameter space. If the parameters are on the boundaries
#' of the parameter space, the identifiability of the model is critical. This
#' may happen especially with splines parameters that may be too close to 0
#' (lower boundary) or classmb parameters that are too high or low (perfect
#' classification). When identifiability of some parameters is suspected, the
#' program can be run again from the former estimates by fixing the suspected
#' parameters to their value with option posfix. This usually solves the
#' problem. An alternative is to remove the parameters of the Beta of Splines
#' link function from the inverse of the Hessian with option partialH. If not,
#' the program should be run again with other initial values, with a higher
#' maximum number of iterations or less strict convergence tolerances.
#'
#' @param fixed two-sided linear formula object for the fixed-effects in the
#' linear mixed model. The response outcome is on the left of \code{~} and the
#' covariates are separated by \code{+} on the right of \code{~}. By default,
#' an intercept is included. If no intercept, \code{-1} should be the first
#' term included on the right of \code{~}.
#' @param mixture one-sided formula object for the class-specific fixed effects
#' in the linear mixed model (to specify only for a number of latent classes
#' greater than 1). Among the list of covariates included in \code{fixed}, the
#' covariates with class-specific regression parameters are entered in
#' \code{mixture} separated by \code{+}. By default, an intercept is included.
#' If no intercept, \code{-1} should be the first term included.
#' @param random optional one-sided formula for the random-effects in the
#' linear mixed model. Covariates with a random-effect are separated by
#' \code{+}. By default, an intercept is included. If no intercept, \code{-1}
#' should be the first term included.
#' @param subject name of the covariate representing the grouping structure
#' specified with ''.
#' @param classmb optional one-sided formula describing the covariates in the
#' class-membership multinomial logistic model. Covariates included are
#' separated by \code{+}. By default, classmb=~1 if ng>1.
#' @param ng optional number of latent classes considered. If \code{ng=1} (by
#' default) no \code{mixture} nor \code{classmb} should be specified. If
#' \code{ng>1}, \code{mixture} is required.
#' @param idiag optional logical for the structure of the variance-covariance
#' matrix of the random-effects. If \code{FALSE}, a non structured matrix of
#' variance-covariance is considered (by default). If \code{TRUE} a diagonal
#' matrix of variance-covariance is considered.
#' @param nwg optional logical indicating if the variance-covariance of the
#' random-effects is class-specific. If \code{FALSE} the variance-covariance
#' matrix is common over latent classes (by default). If \code{TRUE} a
#' class-specific proportional parameter multiplies the variance-covariance
#' matrix in each class (the proportional parameter in the last latent class
#' equals 1 to ensure identifiability).
#' @param cor optional brownian motion or autoregressive process modeling the
#' correlation between the observations. "BM" or "AR" should be specified,
#' followed by the time variable between brackets. By default, no correlation
#' is added.
#' @param data optional data frame containing the variables named in
#' \code{fixed}, \code{mixture}, \code{random}, \code{classmb} and
#' \code{subject}.
#' @param B optional specification for the initial values for the parameters.
#' Three options are allowed: (1) a vector of initial values is entered (the
#' order in which the parameters are included is detailed in \code{details}
#' section). (2) nothing is specified. A preliminary analysis involving the
#' estimation of a standard linear mixed model is performed to choose initial
#' values. (3) when ng>1, a hlme object is entered. It should correspond to
#' the exact same structure of model but with ng=1. The program will
#' automatically generate initial values from this model. This specification
#' avoids the preliminary analysis indicated in (2). Note that due to possible
#' local maxima, the \code{B} vector should be specified and several different
#' starting points should be tried.
#' @param convB optional threshold for the convergence criterion based on the
#' parameter stability. By default, convB=0.0001.
#' @param convL optional threshold for the convergence criterion based on the
#' log-likelihood stability. By default, convL=0.0001.
#' @param convG optional threshold for the convergence criterion based on the
#' derivatives. By default, convG=0.0001.
#' @param prior optional name of a covariate containing a prior information
#' about the latent class membership. The covariate should be an integer with
#' values in 0,1,...,ng. Value 0 indicates no prior for the subject while a
#' value in 1,...,ng indicates that the subject belongs to the corresponding
#' latent class.
#' @param pprior optional vector specifying the names of the covariates containing the
#' prior probabilities to belong to each latent class. These probabilities should be
#' between 0 and 1 and should sum up to 1 for each subject.
#' @param maxiter optional maximum number of iterations for the Marquardt
#' iterative algorithm. By default, maxiter=500.
#' @param subset a specification of the rows to be used: defaults to all rows.
#' This can be any valid indexing vector for the rows of data or if that is not
#' supplied, a data frame made up of the variable used in formula.
#' @param na.action Integer indicating how NAs are managed. The default is 1
#' for 'na.omit'. The alternative is 2 for 'na.fail'. Other options such as
#' 'na.pass' or 'na.exclude' are not implemented in the current version.
#' @param posfix Optional vector specifying the indices in vector B of the
#' parameters that should not be estimated. Default to NULL, all parameters are
#' estimated.
#' @param verbose logical indicating if information about computation should be
#' reported. Default to TRUE.
#' @param returndata logical indicating if data used for computation should be
#' returned. Default to FALSE, data are not returned.
#' @param var.time optional character indicating the name of the time variable.
#' @param partialH optional logical indicating if parameters can be dropped from the
#' Hessian matrix to define convergence criteria.
#' @param nproc the number cores for parallel computation.
#' Default to 1 (sequential mode).
#' @param clustertype optional character indicating the type of cluster for parallel computation.
#' @return The list returned is:
#' \item{ns}{number of grouping units in the dataset}
#' \item{ng}{number of latent classes}
#' \item{loglik}{log-likelihood of the model} \
#' item{best}{vector of parameter estimates in the same order as
#' specified in \code{B} and detailed in section \code{details}}
#' \item{V}{if the model converged (conv=1 or 3), vector containing the upper triangle
#' matrix of variance-covariance estimates of \code{Best} with exception for
#' variance-covariance parameters of the random-effects for which \code{V} contains
#' the variance-covariance estimates of the Cholesky transformed parameters displayed in
#' \code{cholesky}.
#' If conv=2, \code{V} contains the second derivatives of the log-likelihood.}
#' \item{gconv}{vector of convergence criteria: 1. on the parameters, 2. on the
#' likelihood, 3. on the derivatives}
#' \item{conv}{status of convergence: =1 if the convergence criteria were satisfied,
#' =2 if the maximum number of iterations was reached, =3 if the convergence criteria were
#' satisfied with a partial Hessian matrix, =4 or 5 if a problem occured
#' during optimisation}
#' \item{call}{the matched call}
#' \item{niter}{number of Marquardt iterations}
#' \item{N}{internal information used in related functions}
#' \item{idiag}{internal information used in related functions}
#' \item{pred}{table of individual predictions and residuals; it
#' includes marginal predictions (pred_m), marginal residuals (resid_m),
#' subject-specific predictions (pred_ss) and subject-specific residuals
#' (resid_ss) averaged over classes, the observation (obs) and finally the
#' class-specific marginal and subject-specific predictions (with the number of
#' the latent class: pred_m_1,pred_m_2,...,pred_ss_1,pred_ss_2,...). If \code{var.time}
#' is specified, the corresponding measurement time is also included.}
#' \item{pprob}{table of posterior classification and posterior individual
#' class-membership probabilities}
#' \item{Xnames}{list of covariates included in the model}
#' \item{predRE}{table containing individual predictions of the random-effects
#' : a column per random-effect, a line per subject}
#' \item{cholesky}{vector containing the estimates of the Cholesky transformed
#' parameters of the variance-covariance matrix of the random-effects}
#' \item{data}{the original data set (if returndata is TRUE)}
#' @author Cecile Proust-Lima, Benoit Liquet and Viviane Philipps
#'
#' \email{cecile.proust-lima@@inserm.fr}
#' @seealso \code{\link{postprob}}, \code{\link{plot.hlme}},
#' \code{\link{summary}}, \code{\link{predictY}}
#' @references
#'
#' Proust-Lima C, Philipps V, Liquet B (2017). Estimation of Extended Mixed
#' Models Using Latent Classes and Latent Processes: The R Package lcmm.
#' Journal of Statistical Software, 78(2), 1-56. doi:10.18637/jss.v078.i02
#'
#' Verbeke G and Lesaffre E (1996). A linear mixed-effects model with
#' heterogeneity in the random-effects population. Journal of the American
#' Statistical Association 91, 217-21
#'
#' Muthen B and Shedden K (1999). Finite mixture modeling with mixture outcomes
#' using the EM algorithm. Biometrics 55, 463-9
#'
#' Proust C and Jacqmin-Gadda H (2005). Estimation of linear mixed models with
#' a mixture of distribution for the random-effects. Computer Methods Programs
#' Biomedicine 78, 165-73
#' @examples
#'
#'
#' ##### Example of a latent class model estimated for a varying number
#' # of latent classes:
#' # The model includes a subject- (ID) and class-specific linear
#' # trend (intercept and Time in fixed, random and mixture components)
#' # and a common effect of X1 and its interaction with time over classes
#' # (in fixed).
#' # The variance of the random intercept and slope are assumed to be equal
#' # over classes (nwg=F).
#' # The covariate X3 predicts the class membership (in classmb).
#' #
#' # !CAUTION: initialization of mixed models with latent classes is
#' # of most importance because of the problem of multimodality of the likelihood.
#' # Calls m2a-m2d illustrate the different implementations for the
#' # initial values.
#'
#' ### homogeneous linear mixed model (standard linear mixed model)
#' ### with correlated random-effects
#' m1<-hlme(Y~Time*X1,random=~Time,subject='ID',ng=1,data=data_hlme)
#' summary(m1)
#'
#' ### latent class linear mixed model with 2 classes
#'
#' # a. automatic specification from G=1 model estimates:
#' m2a<-hlme(Y~Time*X1,mixture=~Time,random=~Time,classmb=~X2+X3,subject='ID',
#' ng=2,data=data_hlme,B=m1)
#'
#' # b. vector of initial values provided by the user:
#' m2b<-hlme(Y~Time*X1,mixture=~Time,random=~Time,classmb=~X2+X3,subject='ID',
#' ng=2,data=data_hlme,B=c(0.11,-0.74,-0.07,20.71,
#' 29.39,-1,0.13,2.45,-0.29,4.5,0.36,0.79,0.97))
#'
#' # c. random draws from G = 1 model estimates:
#' m2c<-hlme(Y~Time*X1,mixture=~Time,random=~Time,classmb=~X2+X3,subject='ID',
#' ng=2,data=data_hlme,B=random(m1))
#'
#' # d. gridsearch with 50 departures and 10 iterations of the algorithm
#' # (see function gridsearch for details)
#' \dontrun{
#' m2d <- gridsearch(rep = 50, maxiter = 10, minit = m1, hlme(Y ~ Time * X1,
#' mixture =~ Time, random =~ Time, classmb =~ X2 + X3, subject = 'ID', ng = 2,
#' data = data_hlme))
#'
#' }
#'
#'
#'
#' # summary of the estimation process
#' summarytable(m1, m2a, m2b, m2c)
#'
#' # summary of m2a
#' summary(m2a)
#'
#' # posterior classification
#' postprob(m2a)
#'
#' # plot of predicted trajectories using some newdata
#' newdata<-data.frame(Time=seq(0,5,length=100),
#' X1=rep(0,100),X2=rep(0,100),X3=rep(0,100))
#' plot(predictY(m2a,newdata,var.time="Time"),legend.loc="right",bty="l")
#'
#'
#'
#' @export
#'
#'
#'
hlme <-
function(fixed,mixture,random,subject,classmb,ng=1,idiag=FALSE,nwg=FALSE,cor=NULL,data,B,convB=0.0001,convL=0.0001,convG=0.0001,prior,pprior=NULL,maxiter=500,subset=NULL,na.action=1,posfix=NULL,verbose=FALSE,returndata=FALSE,var.time=NULL,partialH=FALSE,nproc=1,clustertype=NULL){
ptm<-proc.time()
cl <- match.call()
args <- as.list(match.call(hlme))[-1]
nom.subject <- as.character(subject)
#### INCLUSION PRIOR
nom.prior <- as.character(args$prior)
####
if(!missing(mixture) & ng==1) stop("No mixture can be specified with ng=1")
if(missing(mixture) & ng>1) stop("The argument mixture has to be specified for ng > 1")
if(!missing(classmb) & ng==1) stop("No classmb can be specified with ng=1")
if(missing(random)) random <- ~-1
if(missing(fixed)) stop("The argument Fixed must be specified in any model")
if(missing(classmb) & ng==1) classmb <- ~-1
if(missing(classmb) & ng>1) classmb <- ~1
if(missing(mixture)) mixture <- ~-1
if(ng==1&nwg==TRUE) stop ("The argument nwg should be FALSE for ng=1")
if(!inherits(fixed,"formula")) stop("The argument fixed must be a formula")
if(!inherits(mixture,"formula")) stop("The argument mixture must be a formula")
if(!inherits(random,"formula")) stop("The argument random must be a formula")
if(!inherits(classmb,"formula")) stop("The argument classmb must be a formula")
if(missing(data)){ stop("The argument data should be specified and defined as a data.frame")}
if(nrow(data)==0) stop("Data should not be empty")
if(missing(subject)){ stop("The argument subject must be specified in any model even without random-effects")}
if(!is.numeric(data[[subject]])) stop("The argument subject must be numeric")
if(!(na.action%in%c(1,2)))stop("only 1 for 'na.omit' or 2 for 'na.fail' are required in na.action argument")
if(na.action==1){
na.action=na.omit
}else{
na.action=na.fail
}
## garder data tel quel pour le renvoyer
if(returndata==TRUE)
{
datareturn <- data
}
else
{
datareturn <- NULL
}
### test de l'argument cor
ncor0 <- 0
cor.type <- cl$cor[1]
cor.time <- cl$cor[2]
cor <- paste(cor.type,cor.time,sep="-")
if (all.equal(cor,character(0))!=TRUE)
{
if (substr(cor,1,2)=="AR") { ncor0 <- 2 }
else if (substr(cor,1,2)=="BM") { ncor0 <- 1 }
else { stop("The argument cor must be of type AR or BM") }
if(!(strsplit(cor,"-")[[1]][2] %in% colnames(data))) stop("Unable to find time variable from argument cor in data")
else { cor.var.time <- strsplit(cor,"-")[[1]][2] }
}
### fin test argument cor
### ad 2/04/2012
X0.names2 <- c("intercept")
### ad
int.fixed <- 0
int.mixture <- 0
int.random <- 0
#7/05/2012
### Traitement des donnees manquantes
# fixed
m <- match.call()[c(1,match(c("data","subset","na.action"),names(match.call()),0))]
m$formula <- terms(fixed)
m$na.action=na.action
m[[1]] <- as.name("model.frame")
m <- eval(m, sys.parent())
na.fixed <- attr(m,"na.action")
# mixture
if(mixture[[2]] != "-1"){
m <- match.call()[c(1,match(c("data","subset","na.action"),names(match.call()),0))]
m$formula <- terms(mixture)
m$na.action <- na.action
m[[1]] <- as.name("model.frame")
m <- eval(m, sys.parent())
na.mixture <- attr(m,"na.action")
}else{
na.mixture <- NULL
}
# random
if(random[[2]] != "-1"){
m <- match.call()[c(1,match(c("data","subset","na.action"),names(match.call()),0))]
m$formula <- terms(random)
m$na.action <- na.action
m[[1]] <- as.name("model.frame")
m <- eval(m, sys.parent())
na.random <- attr(m,"na.action")
}else{
na.random <- NULL
}
# classmb
if(classmb[[2]] != "-1"){
m <- match.call()[c(1,match(c("data","subset","na.action"),names(match.call()),0))]
m$formula <- terms(classmb)
m$na.action <- na.action
m[[1]] <- as.name("model.frame")
m <- eval(m, sys.parent())
na.classmb <- attr(m,"na.action")
}else{
na.classmb <- NULL
}
#cor
if(ncor0!=0)
{
m <- match.call()[c(1,match(c("data","subset","na.action"),names(match.call()),0))]
m$formula <- as.formula(paste(cor.var.time,1,sep="~"))
m$na.action <- na.action
m[[1]] <- as.name("model.frame")
m <- eval(m,sys.parent())
na.cor <- attr(m,"na.action")
}
else { na.cor <- NULL }
### names of covariate in intial fit (sans les interactions)
X0.names2 <- unique(c(X0.names2,colnames(get_all_vars(formula(terms(fixed)),data=data))[-1]))
if(mixture[[2]] != "-1")X0.names2 <- unique(c(X0.names2,colnames(get_all_vars(formula(terms(mixture)),data=data))))
if(random[[2]] != "-1")X0.names2 <- unique(c(X0.names2,colnames(mtemp <- get_all_vars(formula(terms(random)),data=data))))
#7/05/2012
if(classmb[[2]] != "-1")X0.names2 <- unique(c(X0.names2,colnames(get_all_vars(formula(terms(classmb)),data=data))))
#7/05/2012
## Table sans donnees manquante: newdata
na.action <- unique(c(na.fixed,na.mixture,na.random,na.classmb,na.cor))
## dans na.action, on a les indices des NA dans le subset de data
#prendre le subset :
newdata <- data
if(!isTRUE(all.equal(as.character(cl$subset),character(0))))
{
cc <- cl
cc <- cc[c(1,which(names(cl)=="subset"))]
cc[[1]] <- as.name("model.frame")
cc$formula <- formula(paste("~",paste(colnames(data),collapse="+")))
cc$data <- data
cc$na.action <- na.pass
newdata <- eval(cc)
}
#enlever les NA
if(!is.null(na.action)){
newdata <- newdata[-na.action,]
}
attributes(newdata)$terms <- NULL
## Construction de nouvelles var explicatives sur la nouvelle table
## fixed
X_fixed <- model.matrix(fixed,data=newdata)
if(any(colnames(X_fixed)=="(Intercept)")){
ii <- which(colnames(X_fixed)=="(Intercept)")
colnames(X_fixed)[ii] <- "intercept"
int.fixed <- 1
}
nom.fixed <- inddepvar.fixed <- inddepvar.fixed.nom <- colnames(X_fixed)
if(int.fixed>0)inddepvar.fixed <- inddepvar.fixed[-ii]
## mixture
if(mixture[[2]] != "-1"){
X_mixture <- model.matrix(mixture,data=newdata)
if(any(colnames(X_mixture)=="(Intercept)")){
ii <- which(colnames(X_mixture)=="(Intercept)")
colnames(X_mixture)[ii] <- "intercept"
int.mixture <- 1
}
nom.mixture <- inddepvar.mixture <- inddepvar.mixture.nom <- colnames(X_mixture)
if(int.mixture>0)inddepvar.mixture <- inddepvar.mixture[-ii]
id.X_mixture <- 1
}else{
inddepvar.mixture <- nom.mixture <- inddepvar.mixture.nom <- NULL
id.X_mixture <- 0
}
## random
if(random[[2]] != "-1"){
X_random <- model.matrix(random,data=newdata)
if(any(colnames(X_random)=="(Intercept)")){
ii <- which(colnames(X_random)=="(Intercept)")
colnames(X_random)[ii] <- "intercept"
int.random <- 1
}
inddepvar.random <- inddepvar.random.nom <- colnames(X_random)
if(int.random>0) inddepvar.random <- inddepvar.random[-ii]
id.X_random <- 1
}else{
## ad: add inddepvar.random.nom2 <- NULL 10/04/2012
inddepvar.random <- inddepvar.random.nom <- NULL
id.X_random <- 0
}
## classmb
if(classmb[[2]] != "-1"){
if(attr(terms(classmb),"intercept")==0)
{
classmb <- paste("~",classmb[2],"+1")
}
X_classmb <- model.matrix(as.formula(classmb),data=newdata)
if(any(colnames(X_classmb)=="(Intercept)")){
ii <- which(colnames(X_classmb)=="(Intercept)")
colnames(X_classmb)[ii] <- "intercept"
}
id.X_classmb <- 1
}
else{
id.X_classmb <- 0
}
#7/05/2012
############## COVARIATES ##########################
# intercept is always in inddepvar.classmb
var.exp <- NULL
var.exp <- c(var.exp,colnames(X_fixed))
if(id.X_mixture == 1) var.exp <- c(var.exp,colnames(X_mixture))
if(id.X_random == 1)var.exp <- c(var.exp,colnames(X_random))
if(id.X_classmb == 1)var.exp <- c(var.exp,colnames(X_classmb))
var.exp <- unique(var.exp)
if(ncor0>0)
{ if(!(cor.var.time %in% var.exp))
{var.exp <- c(var.exp, cor.var.time)} #si la varaible de temps dans cor n'est dan sles variables expl, on l'ajoute
}
timeobs <- rep(0, nrow(newdata))
if(!is.null(var.time))
{
timeobs <- newdata[,var.time]
if(any(is.na(timeobs))) stop(paste("Cannot use",var.time,"as time variable because it contains missing data"))
}
#if(!(all(nom.mixture %in% nom.fixed))) stop("The covariates in mixture should be also included in the argument fixed")
# controler si les variables de mixture sont toutes dans fixed :
z.fixed <- strsplit(nom.fixed,split=":",fixed=TRUE)
z.fixed <- lapply(z.fixed,sort)
if(id.X_mixture==1)
{
z.mixture <- strsplit(nom.mixture,split=":",fixed=TRUE)
z.mixture <- lapply(z.mixture,sort)
}
else z.mixture <- list()
if(!all(z.mixture %in% z.fixed)) stop("The covariates in mixture should also be included in the argument fixed")
## var dependante
Y.name <- as.character(attributes(terms(fixed))$variables[2])
Y0 <- newdata[,Y.name]
## var expli
X0 <- X_fixed
oldnames <- colnames(X0)
z.X0 <- strsplit(colnames(X0),split=":",fixed=TRUE)
z.X0 <- lapply(z.X0,sort)
if(id.X_mixture == 1)
{
z.mixture <- strsplit(colnames(X_mixture),split=":",fixed=TRUE)
z.mixture <- lapply(z.mixture,sort)
for(i in 1:length(colnames(X_mixture)))
{
#if((colnames(X_mixture)[i] %in% colnames(X0))==F){
if(!isTRUE(z.mixture[i] %in% z.X0))
{
X0 <- cbind(X0,X_mixture[,i])
colnames(X0) <- c(oldnames, colnames(X_mixture)[i])
oldnames <- colnames(X0)
z.X0 <- strsplit(colnames(X0),split=":",fixed=TRUE)
z.X0 <- lapply(z.X0,sort)
}
}
}
else
{
z.mixture <- list()
}
if(id.X_random == 1)
{
z.random <- strsplit(colnames(X_random),split=":",fixed=TRUE)
z.random <- lapply(z.random,sort)
for(i in 1:length(colnames(X_random)))
{
#if((colnames(X_random)[i] %in% colnames(X0))==F){
if(!isTRUE(z.random[i] %in% z.X0))
{
X0 <- cbind(X0,X_random[,i])
colnames(X0) <- c(oldnames, colnames(X_random)[i])
oldnames <- colnames(X0)
z.X0 <- strsplit(colnames(X0),split=":",fixed=TRUE)
z.X0 <- lapply(z.X0,sort)
}
}
}
else
{
z.random <- list()
}
if(id.X_classmb == 1)
{
z.classmb <- strsplit(colnames(X_classmb),split=":",fixed=TRUE)
z.classmb <- lapply(z.classmb,sort)
for(i in 1:length(colnames(X_classmb)))
{
# if((colnames(X_classmb)[i] %in% colnames(X0))==F){
if(!isTRUE(z.classmb[i] %in% z.X0))
{
X0 <- cbind(X0,X_classmb[,i])
colnames(X0) <- c(oldnames, colnames(X_classmb)[i])
oldnames <- colnames(X0)
z.X0 <- strsplit(colnames(X0),split=":",fixed=TRUE)
z.X0 <- lapply(z.X0,sort)
}
}
}
else
{
z.classmb <- list()
}
if(ncor0>0)
{
if(!(cor.var.time %in% colnames(X0))) #cor.var.time jamais ne interaction donc ok (pas de z.cor)
{
X0 <- cbind(X0, newdata[,cor.var.time])
colnames(X0) <- c(oldnames, cor.var.time)
}
}
#colnames(X0) <- var.exp # a remettre si on enleve les z.fixed etc
#X0 <- X0[,-which(colnames(X0)=="intercept")]
### ad
if((any(is.na(X0))==TRUE)|(any(is.na(Y0))==TRUE))stop("The data should not contain any missing value")
#
n <- dim(data)[1]
#if ((int.fixed+int.random)>0) X0<- cbind(intercept=rep(1,n),X0)
#if (!((int.fixed+int.random)>0) & ("intercept" %in% colnames(X0))) X0 <- as.data.frame(X0[,-which(colnames(X0)=="intercept")])
nom.X0 <- colnames(X0)
nvar.exp <- length(nom.X0)
IND <- newdata[,nom.subject]
#IDnum <- as.numeric(IND)
#### INCLUSION PRIOR
#if(missing(prior)){ PRIOR <- seq(0,length=length(IDnum))}
if(missing(prior)){ PRIOR <- seq(0,length=length(IND))}
if(!missing(prior)){
PRIOR <- newdata[,nom.prior]
PRIOR[(is.na(PRIOR))] <- 0
}
####
##pprior : proba d'appartenance a chaque classe
if(is.null(pprior))
{
pprior0 <- matrix(1, length(IND), ng)
}
else
{
if(ng==1) stop("pprior is only useful if ng>1")
##if(classmb != ~-1) stop("classmb should be ~-1 if pprior is specified")
if(!is.character(pprior)) stop("pprior should be a character vector")
if(length(pprior) != ng) stop("pprior should be of length ng")
if(!all(pprior %in% colnames(newdata))) stop("pprior variables should be included in data")
pprior0 <- newdata[,pprior]
}
ng0 <- ng
idiag0 <- as.integer(idiag)
nwg0 <- as.integer(nwg)
idea0 <- rep(0,nvar.exp)
idprob0 <- rep(0,nvar.exp)
idg0 <- rep(0,nvar.exp)
z.X0 <- strsplit(nom.X0,split=":",fixed=TRUE)
z.X0 <- lapply(z.X0,sort)
for (i in 1:nvar.exp){
#idea0[i] <- nom.X0[i]%in%inddepvar.random.nom
#idprob0[i] <- nom.X0[i]%in%inddepvar.classmb.nom
#if(nom.X0[i]%in%nom.fixed & !(nom.X0[i]%in%nom.mixture)) idg0[i] <- 1
#if(nom.X0[i]%in%nom.fixed & nom.X0[i]%in%nom.mixture) idg0[i] <- 2
idea0[i] <- z.X0[i] %in% z.random
idprob0[i] <- z.X0[i] %in% z.classmb
if((z.X0[i] %in% z.fixed) & !(z.X0[i] %in% z.mixture)) idg0[i] <- 1
if((z.X0[i] %in% z.fixed) & (z.X0[i] %in% z.mixture)) idg0[i] <- 2
}
idcor0 <- rep(0,length(z.X0))
if (ncor0!=0) idcor0 <- as.numeric(nom.X0 %in% cor.var.time)
# on ordonne les donn es suivants la variable IND
#matYX <- cbind(IDnum,IND,PRIOR,Y0,X0)
#matYXord <- matYX[sort.list(matYX[,1]),]
#Y0 <- matYXord[,4]
#X0 <- matYXord[,-c(1,2,3,4)]
#IDnum <- matYXord[,1]
#IND <- matYXord[,2]
matYX <- cbind(IND,timeobs,PRIOR,pprior0,Y0,X0)
matYXord <- matYX[sort.list(matYX[,1]),]
Y0 <- as.numeric(matYXord[,4+ng])
X0 <- apply(matYXord[,-c(1:(4+ng)),drop=FALSE],2,as.numeric)
IND <- matYXord[,1]
timeobs <- matYXord[,2]
#### INCLUSION PRIOR
PRIOR <- as.numeric(matYXord[,3])
PRIOR <-as.integer(as.vector(PRIOR))
####
X0 <- as.numeric(as.matrix(X0))
Y0 <- as.numeric(as.matrix(Y0))
#nmes0<-as.vector(table(IDnum))
nmes0 <- as.vector(table(IND))
ns0 <- length(nmes0)
##### INCLUSION PRIOR
# definition de prior a 0 pour l'analyse G=1
prior2 <- as.integer(rep(0,ns0))
prior0 <- prior2
# si prior pas missing alors mettre dedans la classe a priori. Attention tester q les valeurs sont dans 0, G
if(!missing(prior)){
prior0 <- PRIOR[cumsum(nmes0)]
}
INDuniq <- IND[cumsum(nmes0)]
seqnG <- 0:ng0
if (!(all(prior0 %in% seqnG))) stop ("The argument prior should contain integers between 0 and ng")
#####
## pprior de dim ns
if(!is.null(pprior))
{
pprior0 <- matYX[cumsum(nmes0),3+1:ng, drop=FALSE]
if(any(is.na(pprior0))) stop("pprior contains missing values")
#if(any(pprior0 < 0)) stop("pprior should be non negative")
#if(any(pprior0 > 1)) stop("pprior should be < 1")
#if(!all.equal(as.numeric(apply(pprior0,1,sum)), rep(1,ns0))) stop("pprior should sum up to 1 for all subjects")
}
else
{
pprior0 <- matrix(1, ns0, ng0)
}
pprior0 <- t(pprior0)
loglik <- as.double(0)
ni <- 0
istop <- 0
gconv <-rep(0,3)
ppi0 <- rep(0,ns0*ng0)
nv0<-nvar.exp
nobs0<-length(Y0)
resid_m <- rep(0,nobs0)
resid_ss <- rep(0,nobs0)
pred_m_g <- rep(0,nobs0*ng0)
pred_ss_g <- rep(0,nobs0*ng0)
nea0 <- sum(idea0==1)
predRE <- rep(0,nea0*ns0)
varRE <- rep(0,nea0*(nea0+1)/2*ns0)
##---------------------------------------------------------------------------
##definition du vecteur de parametres + initialisation
##---------------------------------------------------------------------------
## gestion de B=random(mod)
Brandom <- FALSE
tryB <- try(as.numeric(B), silent=TRUE)
if(inherits(tryB, "try-error"))
{
if(length(cl$B)==2)
{
if(!inherits(eval(cl$B[[2]], parent.env(environment())),"hlme")) stop("The model specified in B should be of class hlme")
if(as.character(cl$B[1])!="random") stop("Please use random() to specify random initial values")
Brandom <- TRUE
B <- eval(cl$B[[2]], parent.env(environment()))
if(B$conv != 1) stop("Model in argument B did not converge properly")
#if(length(posfix)) stop("Argument posfix is not compatible with random intial values")
}
}
#####cas 1 : ng=1
b <- NULL
b1 <- NULL
NPROB <- 0
if(ng0==1| missing(B)){
NEF<-sum(idg0!=0)
b1[1:NEF] <- 0
if(int.fixed > 0) b1[1] <- mean(Y0)
if(idiag0==1){
NVC<-sum(idea0==1)
b1[(NEF+1):(NEF+NVC)] <- 1}
if(idiag0==0){
kk<-sum(idea0==1)
NVC<-(kk*(kk+1))/2
indice<-cumsum(1:kk)
bidiag<-rep(0,NVC)
bidiag[indice] <- 1
b1[(NEF+1):(NEF+NVC)] <- bidiag
}
if(ncor0==1)
{b1[NEF+NVC+1] <- 1 }
if(ncor0==2)
{b1[(NEF+NVC+1):(NEF+NVC+ncor0)] <- c(0,1) }
b1[NEF+NVC+ncor0+1] <- 1
NPM <- length(b1)
NW <- 0
V <- rep(0,NPM*(NPM+1)/2)
}
#####cas 2 : ng>=2
if(ng0>1){
NPROB <- (sum(idprob0==1))*(ng0-1)
if(NPROB>0) b[1:NPROB] <- 0
NEF <- sum(idg0==1)+(sum(idg0==2))*ng0
if(idiag0==1) NVC <- sum(idea0==1)
if(idiag0==0){
kk <- sum(idea0==1)
NVC <- (kk*(kk+1))/2}
NW <- nwg0*(ng0-1)
if(NW>0) b[(NPROB+NEF+NVC+1):(NPROB+NEF+NVC+NW)] <- 1
if(ncor0==1)
{b[NEF+NVC+NW+1] <- 1 }
if(ncor0==2)
{b[(NPROB+NEF+NVC+NW+1):(NPROB+NEF+NVC+NW+ncor0)] <- c(0,1) }
NPM <- NPROB+NEF+NVC+NW+ncor0+1
V <- rep(0,NPM*(NPM+1)/2)
}
## prm fixes
fix0 <- rep(0,NPM)
if(length(posfix))
{
if(any(!(posfix %in% 1:NPM))) stop("Indexes in posfix are not correct")
fix0[posfix] <- 1
}
if(length(posfix)==NPM) stop("No parameter to estimate")
## pour H restreint
Hr0 <- as.numeric(partialH)
pbH0 <- rep(0,NPM)
if(is.logical(partialH))
{
if(partialH) pbH0 <- rep(1,NPM)
#pbH0[posfix] <- 0
if(sum(pbH0)==0 & Hr0==1) stop("No partial Hessian matrix can be defined")
}
else
{
if(!all(Hr0 %in% 1:NPM)) stop("Indexes in partialH are not correct")
pbH0[Hr0] <- 1
#pbH0[posfix] <- 0
}
indexHr <- NULL
if(sum(pbH0)>0)
{
if(length(posfix)) pbH1 <- pbH0[-posfix] else pbH1 <- pbH0
indexHr <- which(pbH1==1)
}
if(missing(B)){
if(ng0==1 ){
b <- b1
} else {
stop("Please specify initial values with argument 'B'")
}
}
else
{
if(is.vector(B))
{
if(length(B)!=NPM)
{
stop(paste("Vector B should be of length",NPM))
}
else
{
b <-B
if(NVC>0)
{
## remplacer varcov des EA par les prm a estimer
if(idiag0==1)
{
b[NPROB+NEF+1:NVC] <- sqrt(b[NPROB+NEF+1:NVC])
}
else
{
varcov <- matrix(0,nrow=nea0,ncol=nea0)
varcov[upper.tri(varcov,diag=TRUE)] <- b[NPROB+NEF+1:NVC]
varcov <- t(varcov)
varcov[upper.tri(varcov,diag=TRUE)] <- b[NPROB+NEF+1:NVC]
ch <- chol(varcov)
b[NPROB+NEF+1:NVC] <- ch[upper.tri(ch,diag=TRUE)]
}
}
}
}
else
{
if(!inherits(B,"hlme")) stop("B should be either a vector or an object of class hlme")
## B est le meme modele mais pr ng=1 :
if(ng>1 & B$ng==1)
{
NEF2 <- sum(idg0!=0)
NPM2 <- NEF2+NVC+ncor0+1
if(length(B$best)!=NPM2) stop("B is not correct")
if(Brandom==FALSE)
{
## B deterministe
l <- 0
t <- 0
for (i in 1:nvar.exp)
{
if(idg0[i]==1)
{
l <- l+1
t <- t+1
b[NPROB+t] <- B$best[l]
}
if(idg0[i]==2)
{
l <- l+1
for (g in 1:ng)
{
t <- t+1
if(B$conv==1) b[NPROB+t] <- B$best[l]+(g-(ng+1)/2)*sqrt(B$V[l*(l+1)/2])
else b[NPROB+t] <- B$best[l]+(g-(ng+1)/2)*B$best[l]
}
}
}
if(NVC>0)
{
if(idiag==TRUE)
{
b[(NPROB+NEF+1):(NPROB+NEF+NVC)] <- B$cholesky[(1:nea0)*(2:(nea0+1))/2]
}
else
{
b[(NPROB+NEF+1):(NPROB+NEF+NVC)] <- B$cholesky
}
}
if (ncor0>0) {b[(NPROB+NEF+NVC+NW+1):(NPROB+NEF+NVC+NW+ncor0)] <- B$best[(NPM2-ncor0):(NPM2-1)]}
b[NPROB+NEF+NVC+NW+ncor0+1] <- B$best[NPM2]
}
else
{
## B random
bb <- rep(0,NPM-NPROB-NW)
vbb <- matrix(0,NPM-NPROB-NW,NPM-NPROB-NW)
VB <- matrix(0,NPM2,NPM2)
VB[upper.tri(VB,diag=TRUE)] <- B$V
VB <- t(VB)
VB[upper.tri(VB,diag=TRUE)] <- B$V
nbg <- idg0[which(idg0!=0)]
nbg[which(nbg==2)] <- ng
nbgnef <- unlist(sapply(nbg,function(k) if(k>1) rep(2,k) else k))
vbb[which(nbgnef==1),setdiff(1:ncol(vbb),which(nbgnef!=1))] <- VB[which(nbg==1),setdiff(1:ncol(VB),which(nbg!=1))]
vbb[(NEF+1):nrow(vbb),(NEF+1):ncol(vbb)] <- VB[(NEF2+1):nrow(VB),(NEF2+1):ncol(VB)]
l <- 0
t <- 0
for (i in 1:nvar.exp)
{
if(idg0[i]==1)
{
l <- l+1
t <- t+1
bb[t] <- B$best[l]
}
if(idg0[i]==2)
{
l <- l+1
for (g in 1:ng)
{
t <- t+1
bb[t] <- B$best[l]
vbb[t,t] <- VB[l,l]
}
}
}
if(NVC>0)
{
if(idiag==TRUE)
{
bb[NEF+1:NVC] <- B$cholesky[(1:nea0)*(2:(nea0+1))/2]
}
else
{
bb[NEF+1:NVC] <- B$cholesky
}
}
if (ncor0>0)
{
bb[NEF+NVC+1:ncor0] <- B$best[(NPM2-ncor0):(NPM2-1)]
}
bb[NEF+NVC+ncor0+1] <- B$best[NPM2]
up <- vbb[upper.tri(vbb,diag=TRUE)]
vbb <- t(vbb)
vbb[upper.tri(vbb,diag=TRUE)] <- up
##Chol <- chol(vbb)
##Chol <- t(Chol)
b[c((NPROB+1):(NPROB+NEF+NVC),(NPROB+NEF+NVC+NW+1):NPM)] <- rmvnorm(n=1,mean=bb,sigma = vbb) #bb + Chol %*% rnorm(NPM-NPROB-NW)
if(NPROB>0) b[1:NPROB] <- 0
if(NW>0) b[NPROB+NEF+NVC+1:NW] <- 1
}
}
}
}
## faire wRandom et b0Random
NEF2 <- sum(idg0!=0)
NPM2 <- NEF2+NVC+ncor0+1
wRandom <- rep(0,NPM)
b0Random <- NULL
if(NPROB>0) b0Random <- rep(0,ng-1) # nprob
l <- 0
t <- 0
for (i in 1:nvar.exp)
{
if(idg0[i]==1)
{
l <- l+1
t <- t+1
wRandom[NPROB+t] <- l
}
if(idg0[i]==2)
{
l <- l+1
for (g in 1:ng)
{
t <- t+1
wRandom[NPROB+t] <- l
}
}
}
if(NVC>0)
{
wRandom[NPROB+NEF+1:NVC] <- NEF2+1:NVC
}
if(NW>0)
{
b0Random <- c(b0Random,rep(1,ng-1))
}
if (ncor0>0) {wRandom[NPROB+NEF+NVC+NW+1:ncor0] <- NEF2+NVC+1:ncor0}
wRandom[NPM] <- NPM2
## wRandom et b0Random ok.
##------------------------------------------
##------nom au vecteur best
##--------------------------------------------
if(ng0>=2 & NPROB>0){
nom <-rep(nom.X0[idprob0==1],each=ng0-1)
nom1 <- paste(nom," class",c(1:(ng0-1)),sep="")
names(b)[1:NPROB]<-nom1
}
if(ng0==1) names(b)[1:NEF] <- nom.X0[idg0!=0]
if(ng0>1){
nom1<- NULL
for (i in 1:nvar.exp) {
if(idg0[i]==2){ nom <- paste(nom.X0[i]," class",c(1:ng0),sep="")
nom1 <- cbind(nom1,t(nom))}
if(idg0[i]==1) nom1 <- cbind(nom1,nom.X0[i])
}
names(b)[(NPROB+1):(NPROB+NEF)]<- nom1
}
if(NVC!=0)names(b)[(NPROB+NEF+1):(NPROB+NEF+NVC)] <- paste("varcov",c(1:(NVC)))
if(NW!=0)names(b)[(NPROB+NEF+NVC+1):(NPROB+NEF+NVC+NW)] <- paste("varprop class",c(1:(ng0-1)))
names(b)[NPM] <- "stderr"
if(ncor0!=0) {names(b)[(NPROB+NEF+NVC+NW+1):(NPROB+NEF+NVC+NW+ncor0)] <- paste("cor",1:ncor0,sep="") }
names_best <- names(b)
## pas de nprob si pprior
## if(!is.null(pprior))
## {
## if(NPROB>0)
## {
## b <- b[-c(1:NPROB)]
## names_best <- names_best[-c(1:NPROB)]
## fix0 <- fix0[-c(1:NPROB)]
## pbH0 <- pbH0[-c(1:NPROB)]
## NPM <- NPM - NPROB
## indexHr <- NULL
## if(sum(pbH0)>0) indexHr <- which(pbH0==1)
## }
## NPROB <- 0
## }
N <- NULL
N[1] <- NPROB
N[2] <- NEF
N[3] <- NVC
N[4] <- NW
N[5] <- ncor0
idiag <- as.integer(idiag0)
idea <- as.integer(idea0)
nv <- as.integer(nv0)
################ Sortie ###########################
nfix <- sum(fix0)
bfix <- 0
if(nfix>0)
{
bfix <- b[which(fix0==1)]
b <- b[which(fix0==0)]
NPM <- NPM-nfix
}
if(maxiter==0)
{
vrais <- loglikhlme(b,Y0,X0,prior0,pprior0,idprob0,idea0,idg0,idcor0,
ns0,ng0,nv0,nobs0,nea0,nmes0,idiag0,nwg0,ncor0,
NPM,fix0,nfix,bfix)
out <- list(conv=2, V=rep(NA, NPM*(NPM+1)/2), best=b,
ppi2=rep(NA,ns0*ng0), predRE=rep(NA,ns0*nea0),
varRE=rep(NA,ns0*nea0*(nea0+1)/2),
resid_m=rep(NA,nobs0), resid_ss=rep(NA,nobs0),
pred_m_g=rep(NA,nobs0*ng0), pred_ss_g=rep(NA,nobs0*ng0),
gconv=rep(NA,3), niter=0, loglik=vrais)
}
else
{
res <- mla(b=b, m=length(b), fn=loglikhlme,
clustertype=clustertype,.packages=NULL,
epsa=convB,epsb=convL,epsd=convG,partialH=indexHr,
digits=8,print.info=verbose,blinding=FALSE,
multipleTry=25,file="",
nproc=nproc,maxiter=maxiter,minimize=FALSE,
Y0=Y0,X0=X0,prior0=prior0,pprior0=pprior0,idprob0=idprob0,idea0=idea0,
idg0=idg0,idcor0=idcor0,ns0=ns0,ng0=ng0,nv0=nv0,nobs=nobs0,
nea0=nea0,nmes0=nmes0,idiag=idiag0,nwg0=nwg0,ncor0=ncor0,
npm0=NPM,fix0=fix0,nfix0=nfix,bfix0=bfix)
out <- list(conv=res$istop, V=res$v, best=res$b,
ppi2=rep(NA,ns0*ng0), predRE=rep(NA,ns0*nea0),
varRE=rep(NA,ns0*nea0*(nea0+1)/2),
resid_m=rep(NA,nobs0), resid_ss=rep(NA,nobs0),
pred_m_g=rep(NA,nobs0*ng0), pred_ss_g=rep(NA,nobs0*ng0),
gconv=c(res$ca, res$cb, res$rdm), niter=res$ni,
loglik=res$fn.value)
}
if(out$conv %in% c(1,2,3))
{
estim0 <- 0
ll <- 0
ppi0 <- rep(0,ns0*ng0)
resid_m <- rep(0,nobs0)
resid_ss <- rep(0,nobs0)
pred_m_g <- rep(0,nobs0*ng0)
pred_ss_g <- rep(0,nobs0*ng0)
predRE <- rep(0,ns0*nea0)
varRE <- rep(0,ns0*nea0*(nea0+1)/2)
post <- .Fortran(C_loglikhlme,
as.double(Y0),
as.double(X0),
as.integer(prior0),
as.double(pprior0),
as.integer(idprob0),
as.integer(idea0),
as.integer(idg0),
as.integer(idcor0),
as.integer(ns0),
as.integer(ng0),
as.integer(nv0),
as.integer(nobs0),
as.integer(nea0),
as.integer(nmes0),
as.integer(idiag0),
as.integer(nwg0),
as.integer(ncor0),
as.integer(NPM),
as.double(out$best),
ppi=as.double(ppi0),
resid_m=as.double(resid_m),
resid_ss=as.double(resid_ss),
pred_m_g=as.double(pred_m_g),
pred_ss_g=as.double(pred_ss_g),
predRE=as.double(predRE),
varRE=as.double(varRE),
as.integer(fix0),
as.integer(nfix),
as.double(bfix),
as.integer(estim0),
loglik=as.double(ll))
out$ppi2 <- post$ppi
out$predRE <- post$predRE
out$varRE <- post$varRE
out$resid_m <- post$resid_m
out$resid_ss <- post$resid_ss
out$pred_m_g <- post$pred_m_g
out$pred_ss_g <- post$pred_ss_g
}
## creer best a partir de b et bfix
best <- rep(NA,length(fix0))
best[which(fix0==0)] <- out$best
best[which(fix0==1)] <- bfix
names(best) <- names_best
out$best <- best
NPM <- NPM+nfix
### mettre 0 pr les prm fixes
if(length(posfix))
{
if(out$conv==3)
{
mr <- NPM-sum(pbH0)-length(posfix)
Vr <- matrix(0,mr,mr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
Vr <- t(Vr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
V <- matrix(NA,NPM,NPM)
V[setdiff(1:NPM,c(which(pbH0==1),posfix)),setdiff(1:NPM,c(which(pbH0==1),posfix))] <- Vr
V[,posfix] <- 0
V[posfix,] <- 0
V <- V[upper.tri(V,diag=TRUE)]
}
else
{
mr <- NPM-length(posfix)
Vr <- matrix(0,mr,mr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
Vr <- t(Vr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
V <- matrix(0,NPM,NPM)
V[setdiff(1:NPM,posfix),setdiff(1:NPM,posfix)] <- Vr
V <- V[upper.tri(V,diag=TRUE)]
}
}
else
{
if(out$conv==3)
{
mr <- NPM-sum(pbH0)
Vr <- matrix(0,mr,mr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
Vr <- t(Vr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
V <- matrix(NA,NPM,NPM)
V[setdiff(1:NPM,which(pbH0==1)),setdiff(1:NPM,which(pbH0==1))] <- Vr
V <- V[upper.tri(V,diag=TRUE)]
}
else
{
V <- out$V
}
}
### Creation du vecteur cholesky
Cholesky <- rep(0,(nea0*(nea0+1)/2))
if(idiag0==0 & NVC>0){
Cholesky[1:NVC] <- out$best[(NPROB+NEF+1):(NPROB+NEF+NVC)]
### Construction de la matrice U
U <- matrix(0,nrow=nea0,ncol=nea0)
U[upper.tri(U,diag=TRUE)] <- Cholesky[1:NVC]
z <- t(U) %*% U
out$best[(NPROB+NEF+1):(NPROB+NEF+NVC)] <- z[upper.tri(z,diag=TRUE)]
}
if(idiag0==1 & NVC>0){
id <- 1:nea0
indice <- rep(id+id*(id-1)/2)
Cholesky[indice] <- out$best[(NPROB+NEF+1):(NPROB+NEF+nea0)]
out$best[(NPROB+NEF+1):(NPROB+NEF+NVC)] <- out$best[(NPROB+NEF+1):(NPROB+NEF+NVC)]**2
}
####################################################
if (nea0>0) {
predRE <- matrix(out$predRE,ncol=nea0,byrow=TRUE)
predRE <- data.frame(INDuniq,predRE)
colnames(predRE) <- c(nom.subject,nom.X0[idea0!=0])
}
if(ng0>1) {
ppi<- matrix(out$ppi2,ncol=ng0,byrow=TRUE)
}
else {
ppi <- matrix(rep(1,ns0),ncol=ng0)
}
chooseClass <- function(ppi)
{
res <- which.max(ppi)
if(!length(res)) res <- 0
return(res)
}
classif <- apply(ppi,1,chooseClass)
if(any(!is.finite(ppi)))
{
classif <- rep(NA,ns0)
iok <- which(is.finite(ppi[,1]))
classif[iok] <- apply(ppi[iok,,drop=FALSE],1,which.max)
}
ppi <- data.frame(INDuniq,classif,ppi)
temp <- paste("prob",1:ng0,sep="")
colnames(ppi) <- c(nom.subject,"class",temp)
rownames(ppi) <- 1:ns0
pred_m_g <- matrix(out$pred_m_g,nrow=nobs0)
pred_ss_g <- matrix(out$pred_ss_g,nrow=nobs0)
pred_m <- Y0-out$resid_m
pred_ss <- Y0-out$resid_ss
pred <- data.frame(IND,pred_m,out$resid_m,pred_ss,out$resid_ss,Y0,pred_m_g,pred_ss_g)
temp<-paste("pred_m",1:ng0,sep="")
temp1<-paste("pred_ss",1:ng0,sep="")
colnames(pred)<-c(nom.subject,"pred_m","resid_m","pred_ss","resid_ss","obs",temp,temp1)
if(!is.null(var.time))
{
pred <- data.frame(IND,pred_m,out$resid_m,pred_ss,out$resid_ss,Y0,pred_m_g,pred_ss_g,timeobs)
colnames(pred)<-c(nom.subject,"pred_m","resid_m","pred_ss","resid_ss","obs",temp,temp1,var.time)
}
### ad 2/04/2012
if (!("intercept" %in% nom.X0)) X0.names2 <- X0.names2[-1]
### ad
## levels = modalites des variables dans X0 (si facteurs)
levelsdata <- vector("list", length(X0.names2))
levelsfixed <- vector("list", length(X0.names2))
levelsrandom <- vector("list", length(X0.names2))
levelsmixture <- vector("list", length(X0.names2))
levelsclassmb <- vector("list", length(X0.names2))
names(levelsdata) <- X0.names2
names(levelsfixed) <- X0.names2
names(levelsrandom) <- X0.names2
names(levelsmixture) <- X0.names2
names(levelsclassmb) <- X0.names2
for(v in X0.names2)
{
if(v == "intercept") next
if(is.factor(data[,v]))
{
levelsdata[[v]] <- levels(data[,v])
}
#else
#{
if(length(grep(paste("factor\\(",v,"\\)",sep=""), fixed)))
{
levelsfixed[[v]] <- levels(as.factor(data[,v]))
}
# else
# {
if(length(grep(paste("factor\\(",v,"\\)",sep=""), random)))
{
levelsrandom[[v]] <- levels(as.factor(data[,v]))
}
# else
# {
if(length(grep(paste("factor\\(",v,"\\)",sep=""), mixture)))
{
levelsmixture[[v]] <- levels(as.factor(data[,v]))
}
# else
# {
if(length(grep(paste("factor\\(",v,"\\)",sep=""), classmb)))
{
levelsclassmb[[v]] <- levels(as.factor(data[,v]))
}
# }
# }
# }
#}
}
levels <- list(levelsdata=levelsdata, levelsfixed=levelsfixed, levelsrandom=levelsrandom,
levelsmixture=levelsmixture, levelsclassmb=levelsclassmb)
cost<-proc.time()-ptm
res <-list(ns=ns0,ng=ng0,idea0=idea0,idprob0=idprob0,idg0=idg0,idcor0=idcor0,loglik=out$loglik,best=out$best,V=V,gconv=out$gconv,conv=out$conv,call=cl,niter=out$niter,N=N,idiag=idiag0,pred=pred,pprob=ppi,predRE=predRE,varRE=out$varRE,Xnames=nom.X0,Xnames2=X0.names2,cholesky=Cholesky,na.action=na.action,AIC=2*(length(out$best)-length(posfix)-out$loglik),BIC=(length(out$best)-length(posfix))*log(ns0)-2*out$loglik,data=datareturn,wRandom=wRandom,b0Random=b0Random, levels=levels, var.time=var.time, runtime=cost[3])
class(res) <-c("hlme")
if(verbose==TRUE) cat("The program took", round(cost[3],2), "seconds \n")
res
}
Try the lcmm package in your browser
Any scripts or data that you put into this service are public.
lcmm documentation built on Oct. 7, 2023, 1:08 a.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 hlme
Documented in hlme
#' Estimation of latent class linear mixed models
#'
#' This function fits linear mixed models and latent class linear mixed models
#' (LCLMM) also known as growth mixture models or heterogeneous linear mixed
#' models. The LCLMM consists in assuming that the population is divided in a
#' finite number of latent classes. Each latent class is characterised by a
#' specific trajectory modelled by a class-specific linear mixed model. Both
#' the latent class membership and the trajectory can be explained according to
#' covariates. This function is limited to a mixture of Gaussian outcomes. For
#' other types of outcomes, please see function \code{lcmm}. For multivariate
#' longitudinal outcomes, please see \code{multlcmm}.
#'
#'
#' A. THE VECTOR OF PARAMETERS B
#'
#' The parameters in the vector of initial values \code{B} or equivalently in
#' the vector of maximum likelihood estimates \code{best} are included in the
#' following order:
#'
#' (1) ng-1 parameters are required for intercepts in the latent class
#' membership model, and when covariates are included in \code{classmb}, ng-1
#' paramaters should be entered for each covariate;
#'
#' (2) for all covariates in \code{fixed}, one parameter is required if the
#' covariate is not in \code{mixture}, ng paramaters are required if the
#' covariate is also in \code{mixture};
#'
#' (3) the variance of each random-effect specified in \code{random} (including
#' the intercept) when \code{idiag=TRUE}, or the inferior triangular
#' variance-covariance matrix of all the random-effects when
#' \code{idiag=FALSE};
#'
#' (4) only when \code{nwg=TRUE}, ng-1 parameters are required for the ng-1
#' class-specific proportional coefficients in the variance covariance matrix
#' of the random-effects;
#'
#' (5) when \code{cor} is specified, 1 parameter corresponding to the variance
#' of the Brownian motion should be entered with \code{cor=BM} and 2 parameters
#' corresponding to the correlation and the variance parameters of the
#' autoregressive process should be entered
#'
#' (6) the standard error of the residual error.
#'
#' B. CAUTIONS
#'
#' Some caution should be made when using the program:
#'
#' (1) As the log-likelihood of a latent class model can have multiple maxima,
#' a careful choice of the initial values is crucial for ensuring convergence
#' toward the global maximum. The program can be run without entering the
#' vector of initial values (see point 2). However, we recommend to
#' systematically enter initial values in \code{B} and try different sets of
#' initial values.
#'
#' (2) The automatic choice of initial values we provide requires the
#' estimation of a preliminary linear mixed model. The user should be aware
#' that first, this preliminary analysis can take time for large datatsets and
#' second, that the generated initial values can be very not likely and even
#' may converge slowly to a local maximum. This is the reason why several
#' alternatives exist. The vector of initial values can be directly specified
#' in \code{B} the initial values can be generated (automatically or randomly)
#' from a model with \code{ng=}. Finally, function \code{gridsearch} performs
#' an automatic grid search.
#'
#' (3) Convergence criteria are very strict as they are based on the
#' derivatives of the log-likelihood in addition to the parameter stability and
#' log-likelihood stability. In some cases, the program may not converge and
#' reach the maximum number of iterations fixed at 100. In this case, the user
#' should check that parameter estimates at the last iteration are not on the
#' boundaries of the parameter space. If the parameters are on the boundaries
#' of the parameter space, the identifiability of the model is critical. This
#' may happen especially with splines parameters that may be too close to 0
#' (lower boundary) or classmb parameters that are too high or low (perfect
#' classification). When identifiability of some parameters is suspected, the
#' program can be run again from the former estimates by fixing the suspected
#' parameters to their value with option posfix. This usually solves the
#' problem. An alternative is to remove the parameters of the Beta of Splines
#' link function from the inverse of the Hessian with option partialH. If not,
#' the program should be run again with other initial values, with a higher
#' maximum number of iterations or less strict convergence tolerances.
#'
#' @param fixed two-sided linear formula object for the fixed-effects in the
#' linear mixed model. The response outcome is on the left of \code{~} and the
#' covariates are separated by \code{+} on the right of \code{~}. By default,
#' an intercept is included. If no intercept, \code{-1} should be the first
#' term included on the right of \code{~}.
#' @param mixture one-sided formula object for the class-specific fixed effects
#' in the linear mixed model (to specify only for a number of latent classes
#' greater than 1). Among the list of covariates included in \code{fixed}, the
#' covariates with class-specific regression parameters are entered in
#' \code{mixture} separated by \code{+}. By default, an intercept is included.
#' If no intercept, \code{-1} should be the first term included.
#' @param random optional one-sided formula for the random-effects in the
#' linear mixed model. Covariates with a random-effect are separated by
#' \code{+}. By default, an intercept is included. If no intercept, \code{-1}
#' should be the first term included.
#' @param subject name of the covariate representing the grouping structure
#' specified with ''.
#' @param classmb optional one-sided formula describing the covariates in the
#' class-membership multinomial logistic model. Covariates included are
#' separated by \code{+}. By default, classmb=~1 if ng>1.
#' @param ng optional number of latent classes considered. If \code{ng=1} (by
#' default) no \code{mixture} nor \code{classmb} should be specified. If
#' \code{ng>1}, \code{mixture} is required.
#' @param idiag optional logical for the structure of the variance-covariance
#' matrix of the random-effects. If \code{FALSE}, a non structured matrix of
#' variance-covariance is considered (by default). If \code{TRUE} a diagonal
#' matrix of variance-covariance is considered.
#' @param nwg optional logical indicating if the variance-covariance of the
#' random-effects is class-specific. If \code{FALSE} the variance-covariance
#' matrix is common over latent classes (by default). If \code{TRUE} a
#' class-specific proportional parameter multiplies the variance-covariance
#' matrix in each class (the proportional parameter in the last latent class
#' equals 1 to ensure identifiability).
#' @param cor optional brownian motion or autoregressive process modeling the
#' correlation between the observations. "BM" or "AR" should be specified,
#' followed by the time variable between brackets. By default, no correlation
#' is added.
#' @param data optional data frame containing the variables named in
#' \code{fixed}, \code{mixture}, \code{random}, \code{classmb} and
#' \code{subject}.
#' @param B optional specification for the initial values for the parameters.
#' Three options are allowed: (1) a vector of initial values is entered (the
#' order in which the parameters are included is detailed in \code{details}
#' section). (2) nothing is specified. A preliminary analysis involving the
#' estimation of a standard linear mixed model is performed to choose initial
#' values. (3) when ng>1, a hlme object is entered. It should correspond to
#' the exact same structure of model but with ng=1. The program will
#' automatically generate initial values from this model. This specification
#' avoids the preliminary analysis indicated in (2). Note that due to possible
#' local maxima, the \code{B} vector should be specified and several different
#' starting points should be tried.
#' @param convB optional threshold for the convergence criterion based on the
#' parameter stability. By default, convB=0.0001.
#' @param convL optional threshold for the convergence criterion based on the
#' log-likelihood stability. By default, convL=0.0001.
#' @param convG optional threshold for the convergence criterion based on the
#' derivatives. By default, convG=0.0001.
#' @param prior optional name of a covariate containing a prior information
#' about the latent class membership. The covariate should be an integer with
#' values in 0,1,...,ng. Value 0 indicates no prior for the subject while a
#' value in 1,...,ng indicates that the subject belongs to the corresponding
#' latent class.
#' @param pprior optional vector specifying the names of the covariates containing the
#' prior probabilities to belong to each latent class. These probabilities should be
#' between 0 and 1 and should sum up to 1 for each subject.
#' @param maxiter optional maximum number of iterations for the Marquardt
#' iterative algorithm. By default, maxiter=500.
#' @param subset a specification of the rows to be used: defaults to all rows.
#' This can be any valid indexing vector for the rows of data or if that is not
#' supplied, a data frame made up of the variable used in formula.
#' @param na.action Integer indicating how NAs are managed. The default is 1
#' for 'na.omit'. The alternative is 2 for 'na.fail'. Other options such as
#' 'na.pass' or 'na.exclude' are not implemented in the current version.
#' @param posfix Optional vector specifying the indices in vector B of the
#' parameters that should not be estimated. Default to NULL, all parameters are
#' estimated.
#' @param verbose logical indicating if information about computation should be
#' reported. Default to TRUE.
#' @param returndata logical indicating if data used for computation should be
#' returned. Default to FALSE, data are not returned.
#' @param var.time optional character indicating the name of the time variable.
#' @param partialH optional logical indicating if parameters can be dropped from the
#' Hessian matrix to define convergence criteria.
#' @param nproc the number cores for parallel computation.
#' Default to 1 (sequential mode).
#' @param clustertype optional character indicating the type of cluster for parallel computation.
#' @return The list returned is:
#' \item{ns}{number of grouping units in the dataset}
#' \item{ng}{number of latent classes}
#' \item{loglik}{log-likelihood of the model} \
#' item{best}{vector of parameter estimates in the same order as
#' specified in \code{B} and detailed in section \code{details}}
#' \item{V}{if the model converged (conv=1 or 3), vector containing the upper triangle
#' matrix of variance-covariance estimates of \code{Best} with exception for
#' variance-covariance parameters of the random-effects for which \code{V} contains
#' the variance-covariance estimates of the Cholesky transformed parameters displayed in
#' \code{cholesky}.
#' If conv=2, \code{V} contains the second derivatives of the log-likelihood.}
#' \item{gconv}{vector of convergence criteria: 1. on the parameters, 2. on the
#' likelihood, 3. on the derivatives}
#' \item{conv}{status of convergence: =1 if the convergence criteria were satisfied,
#' =2 if the maximum number of iterations was reached, =3 if the convergence criteria were
#' satisfied with a partial Hessian matrix, =4 or 5 if a problem occured
#' during optimisation}
#' \item{call}{the matched call}
#' \item{niter}{number of Marquardt iterations}
#' \item{N}{internal information used in related functions}
#' \item{idiag}{internal information used in related functions}
#' \item{pred}{table of individual predictions and residuals; it
#' includes marginal predictions (pred_m), marginal residuals (resid_m),
#' subject-specific predictions (pred_ss) and subject-specific residuals
#' (resid_ss) averaged over classes, the observation (obs) and finally the
#' class-specific marginal and subject-specific predictions (with the number of
#' the latent class: pred_m_1,pred_m_2,...,pred_ss_1,pred_ss_2,...). If \code{var.time}
#' is specified, the corresponding measurement time is also included.}
#' \item{pprob}{table of posterior classification and posterior individual
#' class-membership probabilities}
#' \item{Xnames}{list of covariates included in the model}
#' \item{predRE}{table containing individual predictions of the random-effects
#' : a column per random-effect, a line per subject}
#' \item{cholesky}{vector containing the estimates of the Cholesky transformed
#' parameters of the variance-covariance matrix of the random-effects}
#' \item{data}{the original data set (if returndata is TRUE)}
#' @author Cecile Proust-Lima, Benoit Liquet and Viviane Philipps
#'
#' \email{cecile.proust-lima@@inserm.fr}
#' @seealso \code{\link{postprob}}, \code{\link{plot.hlme}},
#' \code{\link{summary}}, \code{\link{predictY}}
#' @references
#'
#' Proust-Lima C, Philipps V, Liquet B (2017). Estimation of Extended Mixed
#' Models Using Latent Classes and Latent Processes: The R Package lcmm.
#' Journal of Statistical Software, 78(2), 1-56. doi:10.18637/jss.v078.i02
#'
#' Verbeke G and Lesaffre E (1996). A linear mixed-effects model with
#' heterogeneity in the random-effects population. Journal of the American
#' Statistical Association 91, 217-21
#'
#' Muthen B and Shedden K (1999). Finite mixture modeling with mixture outcomes
#' using the EM algorithm. Biometrics 55, 463-9
#'
#' Proust C and Jacqmin-Gadda H (2005). Estimation of linear mixed models with
#' a mixture of distribution for the random-effects. Computer Methods Programs
#' Biomedicine 78, 165-73
#' @examples
#'
#'
#' ##### Example of a latent class model estimated for a varying number
#' # of latent classes:
#' # The model includes a subject- (ID) and class-specific linear
#' # trend (intercept and Time in fixed, random and mixture components)
#' # and a common effect of X1 and its interaction with time over classes
#' # (in fixed).
#' # The variance of the random intercept and slope are assumed to be equal
#' # over classes (nwg=F).
#' # The covariate X3 predicts the class membership (in classmb).
#' #
#' # !CAUTION: initialization of mixed models with latent classes is
#' # of most importance because of the problem of multimodality of the likelihood.
#' # Calls m2a-m2d illustrate the different implementations for the
#' # initial values.
#'
#' ### homogeneous linear mixed model (standard linear mixed model)
#' ### with correlated random-effects
#' m1<-hlme(Y~Time*X1,random=~Time,subject='ID',ng=1,data=data_hlme)
#' summary(m1)
#'
#' ### latent class linear mixed model with 2 classes
#'
#' # a. automatic specification from G=1 model estimates:
#' m2a<-hlme(Y~Time*X1,mixture=~Time,random=~Time,classmb=~X2+X3,subject='ID',
#' ng=2,data=data_hlme,B=m1)
#'
#' # b. vector of initial values provided by the user:
#' m2b<-hlme(Y~Time*X1,mixture=~Time,random=~Time,classmb=~X2+X3,subject='ID',
#' ng=2,data=data_hlme,B=c(0.11,-0.74,-0.07,20.71,
#' 29.39,-1,0.13,2.45,-0.29,4.5,0.36,0.79,0.97))
#'
#' # c. random draws from G = 1 model estimates:
#' m2c<-hlme(Y~Time*X1,mixture=~Time,random=~Time,classmb=~X2+X3,subject='ID',
#' ng=2,data=data_hlme,B=random(m1))
#'
#' # d. gridsearch with 50 departures and 10 iterations of the algorithm
#' # (see function gridsearch for details)
#' \dontrun{
#' m2d <- gridsearch(rep = 50, maxiter = 10, minit = m1, hlme(Y ~ Time * X1,
#' mixture =~ Time, random =~ Time, classmb =~ X2 + X3, subject = 'ID', ng = 2,
#' data = data_hlme))
#'
#' }
#'
#'
#'
#' # summary of the estimation process
#' summarytable(m1, m2a, m2b, m2c)
#'
#' # summary of m2a
#' summary(m2a)
#'
#' # posterior classification
#' postprob(m2a)
#'
#' # plot of predicted trajectories using some newdata
#' newdata<-data.frame(Time=seq(0,5,length=100),
#' X1=rep(0,100),X2=rep(0,100),X3=rep(0,100))
#' plot(predictY(m2a,newdata,var.time="Time"),legend.loc="right",bty="l")
#'
#'
#'
#' @export
#'
#'
#'
hlme <-
function(fixed,mixture,random,subject,classmb,ng=1,idiag=FALSE,nwg=FALSE,cor=NULL,data,B,convB=0.0001,convL=0.0001,convG=0.0001,prior,pprior=NULL,maxiter=500,subset=NULL,na.action=1,posfix=NULL,verbose=FALSE,returndata=FALSE,var.time=NULL,partialH=FALSE,nproc=1,clustertype=NULL){
ptm<-proc.time()
cl <- match.call()
args <- as.list(match.call(hlme))[-1]
nom.subject <- as.character(subject)
#### INCLUSION PRIOR
nom.prior <- as.character(args$prior)
####
if(!missing(mixture) & ng==1) stop("No mixture can be specified with ng=1")
if(missing(mixture) & ng>1) stop("The argument mixture has to be specified for ng > 1")
if(!missing(classmb) & ng==1) stop("No classmb can be specified with ng=1")
if(missing(random)) random <- ~-1
if(missing(fixed)) stop("The argument Fixed must be specified in any model")
if(missing(classmb) & ng==1) classmb <- ~-1
if(missing(classmb) & ng>1) classmb <- ~1
if(missing(mixture)) mixture <- ~-1
if(ng==1&nwg==TRUE) stop ("The argument nwg should be FALSE for ng=1")
if(!inherits(fixed,"formula")) stop("The argument fixed must be a formula")
if(!inherits(mixture,"formula")) stop("The argument mixture must be a formula")
if(!inherits(random,"formula")) stop("The argument random must be a formula")
if(!inherits(classmb,"formula")) stop("The argument classmb must be a formula")
if(missing(data)){ stop("The argument data should be specified and defined as a data.frame")}
if(nrow(data)==0) stop("Data should not be empty")
if(missing(subject)){ stop("The argument subject must be specified in any model even without random-effects")}
if(!is.numeric(data[[subject]])) stop("The argument subject must be numeric")
if(!(na.action%in%c(1,2)))stop("only 1 for 'na.omit' or 2 for 'na.fail' are required in na.action argument")
if(na.action==1){
na.action=na.omit
}else{
na.action=na.fail
}
## garder data tel quel pour le renvoyer
if(returndata==TRUE)
{
datareturn <- data
}
else
{
datareturn <- NULL
}
### test de l'argument cor
ncor0 <- 0
cor.type <- cl$cor[1]
cor.time <- cl$cor[2]
cor <- paste(cor.type,cor.time,sep="-")
if (all.equal(cor,character(0))!=TRUE)
{
if (substr(cor,1,2)=="AR") { ncor0 <- 2 }
else if (substr(cor,1,2)=="BM") { ncor0 <- 1 }
else { stop("The argument cor must be of type AR or BM") }
if(!(strsplit(cor,"-")[[1]][2] %in% colnames(data))) stop("Unable to find time variable from argument cor in data")
else { cor.var.time <- strsplit(cor,"-")[[1]][2] }
}
### fin test argument cor
### ad 2/04/2012
X0.names2 <- c("intercept")
### ad
int.fixed <- 0
int.mixture <- 0
int.random <- 0
#7/05/2012
### Traitement des donnees manquantes
# fixed
m <- match.call()[c(1,match(c("data","subset","na.action"),names(match.call()),0))]
m$formula <- terms(fixed)
m$na.action=na.action
m[[1]] <- as.name("model.frame")
m <- eval(m, sys.parent())
na.fixed <- attr(m,"na.action")
# mixture
if(mixture[[2]] != "-1"){
m <- match.call()[c(1,match(c("data","subset","na.action"),names(match.call()),0))]
m$formula <- terms(mixture)
m$na.action <- na.action
m[[1]] <- as.name("model.frame")
m <- eval(m, sys.parent())
na.mixture <- attr(m,"na.action")
}else{
na.mixture <- NULL
}
# random
if(random[[2]] != "-1"){
m <- match.call()[c(1,match(c("data","subset","na.action"),names(match.call()),0))]
m$formula <- terms(random)
m$na.action <- na.action
m[[1]] <- as.name("model.frame")
m <- eval(m, sys.parent())
na.random <- attr(m,"na.action")
}else{
na.random <- NULL
}
# classmb
if(classmb[[2]] != "-1"){
m <- match.call()[c(1,match(c("data","subset","na.action"),names(match.call()),0))]
m$formula <- terms(classmb)
m$na.action <- na.action
m[[1]] <- as.name("model.frame")
m <- eval(m, sys.parent())
na.classmb <- attr(m,"na.action")
}else{
na.classmb <- NULL
}
#cor
if(ncor0!=0)
{
m <- match.call()[c(1,match(c("data","subset","na.action"),names(match.call()),0))]
m$formula <- as.formula(paste(cor.var.time,1,sep="~"))
m$na.action <- na.action
m[[1]] <- as.name("model.frame")
m <- eval(m,sys.parent())
na.cor <- attr(m,"na.action")
}
else { na.cor <- NULL }
### names of covariate in intial fit (sans les interactions)
X0.names2 <- unique(c(X0.names2,colnames(get_all_vars(formula(terms(fixed)),data=data))[-1]))
if(mixture[[2]] != "-1")X0.names2 <- unique(c(X0.names2,colnames(get_all_vars(formula(terms(mixture)),data=data))))
if(random[[2]] != "-1")X0.names2 <- unique(c(X0.names2,colnames(mtemp <- get_all_vars(formula(terms(random)),data=data))))
#7/05/2012
if(classmb[[2]] != "-1")X0.names2 <- unique(c(X0.names2,colnames(get_all_vars(formula(terms(classmb)),data=data))))
#7/05/2012
## Table sans donnees manquante: newdata
na.action <- unique(c(na.fixed,na.mixture,na.random,na.classmb,na.cor))
## dans na.action, on a les indices des NA dans le subset de data
#prendre le subset :
newdata <- data
if(!isTRUE(all.equal(as.character(cl$subset),character(0))))
{
cc <- cl
cc <- cc[c(1,which(names(cl)=="subset"))]
cc[[1]] <- as.name("model.frame")
cc$formula <- formula(paste("~",paste(colnames(data),collapse="+")))
cc$data <- data
cc$na.action <- na.pass
newdata <- eval(cc)
}
#enlever les NA
if(!is.null(na.action)){
newdata <- newdata[-na.action,]
}
attributes(newdata)$terms <- NULL
## Construction de nouvelles var explicatives sur la nouvelle table
## fixed
X_fixed <- model.matrix(fixed,data=newdata)
if(any(colnames(X_fixed)=="(Intercept)")){
ii <- which(colnames(X_fixed)=="(Intercept)")
colnames(X_fixed)[ii] <- "intercept"
int.fixed <- 1
}
nom.fixed <- inddepvar.fixed <- inddepvar.fixed.nom <- colnames(X_fixed)
if(int.fixed>0)inddepvar.fixed <- inddepvar.fixed[-ii]
## mixture
if(mixture[[2]] != "-1"){
X_mixture <- model.matrix(mixture,data=newdata)
if(any(colnames(X_mixture)=="(Intercept)")){
ii <- which(colnames(X_mixture)=="(Intercept)")
colnames(X_mixture)[ii] <- "intercept"
int.mixture <- 1
}
nom.mixture <- inddepvar.mixture <- inddepvar.mixture.nom <- colnames(X_mixture)
if(int.mixture>0)inddepvar.mixture <- inddepvar.mixture[-ii]
id.X_mixture <- 1
}else{
inddepvar.mixture <- nom.mixture <- inddepvar.mixture.nom <- NULL
id.X_mixture <- 0
}
## random
if(random[[2]] != "-1"){
X_random <- model.matrix(random,data=newdata)
if(any(colnames(X_random)=="(Intercept)")){
ii <- which(colnames(X_random)=="(Intercept)")
colnames(X_random)[ii] <- "intercept"
int.random <- 1
}
inddepvar.random <- inddepvar.random.nom <- colnames(X_random)
if(int.random>0) inddepvar.random <- inddepvar.random[-ii]
id.X_random <- 1
}else{
## ad: add inddepvar.random.nom2 <- NULL 10/04/2012
inddepvar.random <- inddepvar.random.nom <- NULL
id.X_random <- 0
}
## classmb
if(classmb[[2]] != "-1"){
if(attr(terms(classmb),"intercept")==0)
{
classmb <- paste("~",classmb[2],"+1")
}
X_classmb <- model.matrix(as.formula(classmb),data=newdata)
if(any(colnames(X_classmb)=="(Intercept)")){
ii <- which(colnames(X_classmb)=="(Intercept)")
colnames(X_classmb)[ii] <- "intercept"
}
id.X_classmb <- 1
}
else{
id.X_classmb <- 0
}
#7/05/2012
############## COVARIATES ##########################
# intercept is always in inddepvar.classmb
var.exp <- NULL
var.exp <- c(var.exp,colnames(X_fixed))
if(id.X_mixture == 1) var.exp <- c(var.exp,colnames(X_mixture))
if(id.X_random == 1)var.exp <- c(var.exp,colnames(X_random))
if(id.X_classmb == 1)var.exp <- c(var.exp,colnames(X_classmb))
var.exp <- unique(var.exp)
if(ncor0>0)
{ if(!(cor.var.time %in% var.exp))
{var.exp <- c(var.exp, cor.var.time)} #si la varaible de temps dans cor n'est dan sles variables expl, on l'ajoute
}
timeobs <- rep(0, nrow(newdata))
if(!is.null(var.time))
{
timeobs <- newdata[,var.time]
if(any(is.na(timeobs))) stop(paste("Cannot use",var.time,"as time variable because it contains missing data"))
}
#if(!(all(nom.mixture %in% nom.fixed))) stop("The covariates in mixture should be also included in the argument fixed")
# controler si les variables de mixture sont toutes dans fixed :
z.fixed <- strsplit(nom.fixed,split=":",fixed=TRUE)
z.fixed <- lapply(z.fixed,sort)
if(id.X_mixture==1)
{
z.mixture <- strsplit(nom.mixture,split=":",fixed=TRUE)
z.mixture <- lapply(z.mixture,sort)
}
else z.mixture <- list()
if(!all(z.mixture %in% z.fixed)) stop("The covariates in mixture should also be included in the argument fixed")
## var dependante
Y.name <- as.character(attributes(terms(fixed))$variables[2])
Y0 <- newdata[,Y.name]
## var expli
X0 <- X_fixed
oldnames <- colnames(X0)
z.X0 <- strsplit(colnames(X0),split=":",fixed=TRUE)
z.X0 <- lapply(z.X0,sort)
if(id.X_mixture == 1)
{
z.mixture <- strsplit(colnames(X_mixture),split=":",fixed=TRUE)
z.mixture <- lapply(z.mixture,sort)
for(i in 1:length(colnames(X_mixture)))
{
#if((colnames(X_mixture)[i] %in% colnames(X0))==F){
if(!isTRUE(z.mixture[i] %in% z.X0))
{
X0 <- cbind(X0,X_mixture[,i])
colnames(X0) <- c(oldnames, colnames(X_mixture)[i])
oldnames <- colnames(X0)
z.X0 <- strsplit(colnames(X0),split=":",fixed=TRUE)
z.X0 <- lapply(z.X0,sort)
}
}
}
else
{
z.mixture <- list()
}
if(id.X_random == 1)
{
z.random <- strsplit(colnames(X_random),split=":",fixed=TRUE)
z.random <- lapply(z.random,sort)
for(i in 1:length(colnames(X_random)))
{
#if((colnames(X_random)[i] %in% colnames(X0))==F){
if(!isTRUE(z.random[i] %in% z.X0))
{
X0 <- cbind(X0,X_random[,i])
colnames(X0) <- c(oldnames, colnames(X_random)[i])
oldnames <- colnames(X0)
z.X0 <- strsplit(colnames(X0),split=":",fixed=TRUE)
z.X0 <- lapply(z.X0,sort)
}
}
}
else
{
z.random <- list()
}
if(id.X_classmb == 1)
{
z.classmb <- strsplit(colnames(X_classmb),split=":",fixed=TRUE)
z.classmb <- lapply(z.classmb,sort)
for(i in 1:length(colnames(X_classmb)))
{
# if((colnames(X_classmb)[i] %in% colnames(X0))==F){
if(!isTRUE(z.classmb[i] %in% z.X0))
{
X0 <- cbind(X0,X_classmb[,i])
colnames(X0) <- c(oldnames, colnames(X_classmb)[i])
oldnames <- colnames(X0)
z.X0 <- strsplit(colnames(X0),split=":",fixed=TRUE)
z.X0 <- lapply(z.X0,sort)
}
}
}
else
{
z.classmb <- list()
}
if(ncor0>0)
{
if(!(cor.var.time %in% colnames(X0))) #cor.var.time jamais ne interaction donc ok (pas de z.cor)
{
X0 <- cbind(X0, newdata[,cor.var.time])
colnames(X0) <- c(oldnames, cor.var.time)
}
}
#colnames(X0) <- var.exp # a remettre si on enleve les z.fixed etc
#X0 <- X0[,-which(colnames(X0)=="intercept")]
### ad
if((any(is.na(X0))==TRUE)|(any(is.na(Y0))==TRUE))stop("The data should not contain any missing value")
#
n <- dim(data)[1]
#if ((int.fixed+int.random)>0) X0<- cbind(intercept=rep(1,n),X0)
#if (!((int.fixed+int.random)>0) & ("intercept" %in% colnames(X0))) X0 <- as.data.frame(X0[,-which(colnames(X0)=="intercept")])
nom.X0 <- colnames(X0)
nvar.exp <- length(nom.X0)
IND <- newdata[,nom.subject]
#IDnum <- as.numeric(IND)
#### INCLUSION PRIOR
#if(missing(prior)){ PRIOR <- seq(0,length=length(IDnum))}
if(missing(prior)){ PRIOR <- seq(0,length=length(IND))}
if(!missing(prior)){
PRIOR <- newdata[,nom.prior]
PRIOR[(is.na(PRIOR))] <- 0
}
####
##pprior : proba d'appartenance a chaque classe
if(is.null(pprior))
{
pprior0 <- matrix(1, length(IND), ng)
}
else
{
if(ng==1) stop("pprior is only useful if ng>1")
##if(classmb != ~-1) stop("classmb should be ~-1 if pprior is specified")
if(!is.character(pprior)) stop("pprior should be a character vector")
if(length(pprior) != ng) stop("pprior should be of length ng")
if(!all(pprior %in% colnames(newdata))) stop("pprior variables should be included in data")
pprior0 <- newdata[,pprior]
}
ng0 <- ng
idiag0 <- as.integer(idiag)
nwg0 <- as.integer(nwg)
idea0 <- rep(0,nvar.exp)
idprob0 <- rep(0,nvar.exp)
idg0 <- rep(0,nvar.exp)
z.X0 <- strsplit(nom.X0,split=":",fixed=TRUE)
z.X0 <- lapply(z.X0,sort)
for (i in 1:nvar.exp){
#idea0[i] <- nom.X0[i]%in%inddepvar.random.nom
#idprob0[i] <- nom.X0[i]%in%inddepvar.classmb.nom
#if(nom.X0[i]%in%nom.fixed & !(nom.X0[i]%in%nom.mixture)) idg0[i] <- 1
#if(nom.X0[i]%in%nom.fixed & nom.X0[i]%in%nom.mixture) idg0[i] <- 2
idea0[i] <- z.X0[i] %in% z.random
idprob0[i] <- z.X0[i] %in% z.classmb
if((z.X0[i] %in% z.fixed) & !(z.X0[i] %in% z.mixture)) idg0[i] <- 1
if((z.X0[i] %in% z.fixed) & (z.X0[i] %in% z.mixture)) idg0[i] <- 2
}
idcor0 <- rep(0,length(z.X0))
if (ncor0!=0) idcor0 <- as.numeric(nom.X0 %in% cor.var.time)
# on ordonne les donn es suivants la variable IND
#matYX <- cbind(IDnum,IND,PRIOR,Y0,X0)
#matYXord <- matYX[sort.list(matYX[,1]),]
#Y0 <- matYXord[,4]
#X0 <- matYXord[,-c(1,2,3,4)]
#IDnum <- matYXord[,1]
#IND <- matYXord[,2]
matYX <- cbind(IND,timeobs,PRIOR,pprior0,Y0,X0)
matYXord <- matYX[sort.list(matYX[,1]),]
Y0 <- as.numeric(matYXord[,4+ng])
X0 <- apply(matYXord[,-c(1:(4+ng)),drop=FALSE],2,as.numeric)
IND <- matYXord[,1]
timeobs <- matYXord[,2]
#### INCLUSION PRIOR
PRIOR <- as.numeric(matYXord[,3])
PRIOR <-as.integer(as.vector(PRIOR))
####
X0 <- as.numeric(as.matrix(X0))
Y0 <- as.numeric(as.matrix(Y0))
#nmes0<-as.vector(table(IDnum))
nmes0 <- as.vector(table(IND))
ns0 <- length(nmes0)
##### INCLUSION PRIOR
# definition de prior a 0 pour l'analyse G=1
prior2 <- as.integer(rep(0,ns0))
prior0 <- prior2
# si prior pas missing alors mettre dedans la classe a priori. Attention tester q les valeurs sont dans 0, G
if(!missing(prior)){
prior0 <- PRIOR[cumsum(nmes0)]
}
INDuniq <- IND[cumsum(nmes0)]
seqnG <- 0:ng0
if (!(all(prior0 %in% seqnG))) stop ("The argument prior should contain integers between 0 and ng")
#####
## pprior de dim ns
if(!is.null(pprior))
{
pprior0 <- matYX[cumsum(nmes0),3+1:ng, drop=FALSE]
if(any(is.na(pprior0))) stop("pprior contains missing values")
#if(any(pprior0 < 0)) stop("pprior should be non negative")
#if(any(pprior0 > 1)) stop("pprior should be < 1")
#if(!all.equal(as.numeric(apply(pprior0,1,sum)), rep(1,ns0))) stop("pprior should sum up to 1 for all subjects")
}
else
{
pprior0 <- matrix(1, ns0, ng0)
}
pprior0 <- t(pprior0)
loglik <- as.double(0)
ni <- 0
istop <- 0
gconv <-rep(0,3)
ppi0 <- rep(0,ns0*ng0)
nv0<-nvar.exp
nobs0<-length(Y0)
resid_m <- rep(0,nobs0)
resid_ss <- rep(0,nobs0)
pred_m_g <- rep(0,nobs0*ng0)
pred_ss_g <- rep(0,nobs0*ng0)
nea0 <- sum(idea0==1)
predRE <- rep(0,nea0*ns0)
varRE <- rep(0,nea0*(nea0+1)/2*ns0)
##---------------------------------------------------------------------------
##definition du vecteur de parametres + initialisation
##---------------------------------------------------------------------------
## gestion de B=random(mod)
Brandom <- FALSE
tryB <- try(as.numeric(B), silent=TRUE)
if(inherits(tryB, "try-error"))
{
if(length(cl$B)==2)
{
if(!inherits(eval(cl$B[[2]], parent.env(environment())),"hlme")) stop("The model specified in B should be of class hlme")
if(as.character(cl$B[1])!="random") stop("Please use random() to specify random initial values")
Brandom <- TRUE
B <- eval(cl$B[[2]], parent.env(environment()))
if(B$conv != 1) stop("Model in argument B did not converge properly")
#if(length(posfix)) stop("Argument posfix is not compatible with random intial values")
}
}
#####cas 1 : ng=1
b <- NULL
b1 <- NULL
NPROB <- 0
if(ng0==1| missing(B)){
NEF<-sum(idg0!=0)
b1[1:NEF] <- 0
if(int.fixed > 0) b1[1] <- mean(Y0)
if(idiag0==1){
NVC<-sum(idea0==1)
b1[(NEF+1):(NEF+NVC)] <- 1}
if(idiag0==0){
kk<-sum(idea0==1)
NVC<-(kk*(kk+1))/2
indice<-cumsum(1:kk)
bidiag<-rep(0,NVC)
bidiag[indice] <- 1
b1[(NEF+1):(NEF+NVC)] <- bidiag
}
if(ncor0==1)
{b1[NEF+NVC+1] <- 1 }
if(ncor0==2)
{b1[(NEF+NVC+1):(NEF+NVC+ncor0)] <- c(0,1) }
b1[NEF+NVC+ncor0+1] <- 1
NPM <- length(b1)
NW <- 0
V <- rep(0,NPM*(NPM+1)/2)
}
#####cas 2 : ng>=2
if(ng0>1){
NPROB <- (sum(idprob0==1))*(ng0-1)
if(NPROB>0) b[1:NPROB] <- 0
NEF <- sum(idg0==1)+(sum(idg0==2))*ng0
if(idiag0==1) NVC <- sum(idea0==1)
if(idiag0==0){
kk <- sum(idea0==1)
NVC <- (kk*(kk+1))/2}
NW <- nwg0*(ng0-1)
if(NW>0) b[(NPROB+NEF+NVC+1):(NPROB+NEF+NVC+NW)] <- 1
if(ncor0==1)
{b[NEF+NVC+NW+1] <- 1 }
if(ncor0==2)
{b[(NPROB+NEF+NVC+NW+1):(NPROB+NEF+NVC+NW+ncor0)] <- c(0,1) }
NPM <- NPROB+NEF+NVC+NW+ncor0+1
V <- rep(0,NPM*(NPM+1)/2)
}
## prm fixes
fix0 <- rep(0,NPM)
if(length(posfix))
{
if(any(!(posfix %in% 1:NPM))) stop("Indexes in posfix are not correct")
fix0[posfix] <- 1
}
if(length(posfix)==NPM) stop("No parameter to estimate")
## pour H restreint
Hr0 <- as.numeric(partialH)
pbH0 <- rep(0,NPM)
if(is.logical(partialH))
{
if(partialH) pbH0 <- rep(1,NPM)
#pbH0[posfix] <- 0
if(sum(pbH0)==0 & Hr0==1) stop("No partial Hessian matrix can be defined")
}
else
{
if(!all(Hr0 %in% 1:NPM)) stop("Indexes in partialH are not correct")
pbH0[Hr0] <- 1
#pbH0[posfix] <- 0
}
indexHr <- NULL
if(sum(pbH0)>0)
{
if(length(posfix)) pbH1 <- pbH0[-posfix] else pbH1 <- pbH0
indexHr <- which(pbH1==1)
}
if(missing(B)){
if(ng0==1 ){
b <- b1
} else {
stop("Please specify initial values with argument 'B'")
}
}
else
{
if(is.vector(B))
{
if(length(B)!=NPM)
{
stop(paste("Vector B should be of length",NPM))
}
else
{
b <-B
if(NVC>0)
{
## remplacer varcov des EA par les prm a estimer
if(idiag0==1)
{
b[NPROB+NEF+1:NVC] <- sqrt(b[NPROB+NEF+1:NVC])
}
else
{
varcov <- matrix(0,nrow=nea0,ncol=nea0)
varcov[upper.tri(varcov,diag=TRUE)] <- b[NPROB+NEF+1:NVC]
varcov <- t(varcov)
varcov[upper.tri(varcov,diag=TRUE)] <- b[NPROB+NEF+1:NVC]
ch <- chol(varcov)
b[NPROB+NEF+1:NVC] <- ch[upper.tri(ch,diag=TRUE)]
}
}
}
}
else
{
if(!inherits(B,"hlme")) stop("B should be either a vector or an object of class hlme")
## B est le meme modele mais pr ng=1 :
if(ng>1 & B$ng==1)
{
NEF2 <- sum(idg0!=0)
NPM2 <- NEF2+NVC+ncor0+1
if(length(B$best)!=NPM2) stop("B is not correct")
if(Brandom==FALSE)
{
## B deterministe
l <- 0
t <- 0
for (i in 1:nvar.exp)
{
if(idg0[i]==1)
{
l <- l+1
t <- t+1
b[NPROB+t] <- B$best[l]
}
if(idg0[i]==2)
{
l <- l+1
for (g in 1:ng)
{
t <- t+1
if(B$conv==1) b[NPROB+t] <- B$best[l]+(g-(ng+1)/2)*sqrt(B$V[l*(l+1)/2])
else b[NPROB+t] <- B$best[l]+(g-(ng+1)/2)*B$best[l]
}
}
}
if(NVC>0)
{
if(idiag==TRUE)
{
b[(NPROB+NEF+1):(NPROB+NEF+NVC)] <- B$cholesky[(1:nea0)*(2:(nea0+1))/2]
}
else
{
b[(NPROB+NEF+1):(NPROB+NEF+NVC)] <- B$cholesky
}
}
if (ncor0>0) {b[(NPROB+NEF+NVC+NW+1):(NPROB+NEF+NVC+NW+ncor0)] <- B$best[(NPM2-ncor0):(NPM2-1)]}
b[NPROB+NEF+NVC+NW+ncor0+1] <- B$best[NPM2]
}
else
{
## B random
bb <- rep(0,NPM-NPROB-NW)
vbb <- matrix(0,NPM-NPROB-NW,NPM-NPROB-NW)
VB <- matrix(0,NPM2,NPM2)
VB[upper.tri(VB,diag=TRUE)] <- B$V
VB <- t(VB)
VB[upper.tri(VB,diag=TRUE)] <- B$V
nbg <- idg0[which(idg0!=0)]
nbg[which(nbg==2)] <- ng
nbgnef <- unlist(sapply(nbg,function(k) if(k>1) rep(2,k) else k))
vbb[which(nbgnef==1),setdiff(1:ncol(vbb),which(nbgnef!=1))] <- VB[which(nbg==1),setdiff(1:ncol(VB),which(nbg!=1))]
vbb[(NEF+1):nrow(vbb),(NEF+1):ncol(vbb)] <- VB[(NEF2+1):nrow(VB),(NEF2+1):ncol(VB)]
l <- 0
t <- 0
for (i in 1:nvar.exp)
{
if(idg0[i]==1)
{
l <- l+1
t <- t+1
bb[t] <- B$best[l]
}
if(idg0[i]==2)
{
l <- l+1
for (g in 1:ng)
{
t <- t+1
bb[t] <- B$best[l]
vbb[t,t] <- VB[l,l]
}
}
}
if(NVC>0)
{
if(idiag==TRUE)
{
bb[NEF+1:NVC] <- B$cholesky[(1:nea0)*(2:(nea0+1))/2]
}
else
{
bb[NEF+1:NVC] <- B$cholesky
}
}
if (ncor0>0)
{
bb[NEF+NVC+1:ncor0] <- B$best[(NPM2-ncor0):(NPM2-1)]
}
bb[NEF+NVC+ncor0+1] <- B$best[NPM2]
up <- vbb[upper.tri(vbb,diag=TRUE)]
vbb <- t(vbb)
vbb[upper.tri(vbb,diag=TRUE)] <- up
##Chol <- chol(vbb)
##Chol <- t(Chol)
b[c((NPROB+1):(NPROB+NEF+NVC),(NPROB+NEF+NVC+NW+1):NPM)] <- rmvnorm(n=1,mean=bb,sigma = vbb) #bb + Chol %*% rnorm(NPM-NPROB-NW)
if(NPROB>0) b[1:NPROB] <- 0
if(NW>0) b[NPROB+NEF+NVC+1:NW] <- 1
}
}
}
}
## faire wRandom et b0Random
NEF2 <- sum(idg0!=0)
NPM2 <- NEF2+NVC+ncor0+1
wRandom <- rep(0,NPM)
b0Random <- NULL
if(NPROB>0) b0Random <- rep(0,ng-1) # nprob
l <- 0
t <- 0
for (i in 1:nvar.exp)
{
if(idg0[i]==1)
{
l <- l+1
t <- t+1
wRandom[NPROB+t] <- l
}
if(idg0[i]==2)
{
l <- l+1
for (g in 1:ng)
{
t <- t+1
wRandom[NPROB+t] <- l
}
}
}
if(NVC>0)
{
wRandom[NPROB+NEF+1:NVC] <- NEF2+1:NVC
}
if(NW>0)
{
b0Random <- c(b0Random,rep(1,ng-1))
}
if (ncor0>0) {wRandom[NPROB+NEF+NVC+NW+1:ncor0] <- NEF2+NVC+1:ncor0}
wRandom[NPM] <- NPM2
## wRandom et b0Random ok.
##------------------------------------------
##------nom au vecteur best
##--------------------------------------------
if(ng0>=2 & NPROB>0){
nom <-rep(nom.X0[idprob0==1],each=ng0-1)
nom1 <- paste(nom," class",c(1:(ng0-1)),sep="")
names(b)[1:NPROB]<-nom1
}
if(ng0==1) names(b)[1:NEF] <- nom.X0[idg0!=0]
if(ng0>1){
nom1<- NULL
for (i in 1:nvar.exp) {
if(idg0[i]==2){ nom <- paste(nom.X0[i]," class",c(1:ng0),sep="")
nom1 <- cbind(nom1,t(nom))}
if(idg0[i]==1) nom1 <- cbind(nom1,nom.X0[i])
}
names(b)[(NPROB+1):(NPROB+NEF)]<- nom1
}
if(NVC!=0)names(b)[(NPROB+NEF+1):(NPROB+NEF+NVC)] <- paste("varcov",c(1:(NVC)))
if(NW!=0)names(b)[(NPROB+NEF+NVC+1):(NPROB+NEF+NVC+NW)] <- paste("varprop class",c(1:(ng0-1)))
names(b)[NPM] <- "stderr"
if(ncor0!=0) {names(b)[(NPROB+NEF+NVC+NW+1):(NPROB+NEF+NVC+NW+ncor0)] <- paste("cor",1:ncor0,sep="") }
names_best <- names(b)
## pas de nprob si pprior
## if(!is.null(pprior))
## {
## if(NPROB>0)
## {
## b <- b[-c(1:NPROB)]
## names_best <- names_best[-c(1:NPROB)]
## fix0 <- fix0[-c(1:NPROB)]
## pbH0 <- pbH0[-c(1:NPROB)]
## NPM <- NPM - NPROB
## indexHr <- NULL
## if(sum(pbH0)>0) indexHr <- which(pbH0==1)
## }
## NPROB <- 0
## }
N <- NULL
N[1] <- NPROB
N[2] <- NEF
N[3] <- NVC
N[4] <- NW
N[5] <- ncor0
idiag <- as.integer(idiag0)
idea <- as.integer(idea0)
nv <- as.integer(nv0)
################ Sortie ###########################
nfix <- sum(fix0)
bfix <- 0
if(nfix>0)
{
bfix <- b[which(fix0==1)]
b <- b[which(fix0==0)]
NPM <- NPM-nfix
}
if(maxiter==0)
{
vrais <- loglikhlme(b,Y0,X0,prior0,pprior0,idprob0,idea0,idg0,idcor0,
ns0,ng0,nv0,nobs0,nea0,nmes0,idiag0,nwg0,ncor0,
NPM,fix0,nfix,bfix)
out <- list(conv=2, V=rep(NA, NPM*(NPM+1)/2), best=b,
ppi2=rep(NA,ns0*ng0), predRE=rep(NA,ns0*nea0),
varRE=rep(NA,ns0*nea0*(nea0+1)/2),
resid_m=rep(NA,nobs0), resid_ss=rep(NA,nobs0),
pred_m_g=rep(NA,nobs0*ng0), pred_ss_g=rep(NA,nobs0*ng0),
gconv=rep(NA,3), niter=0, loglik=vrais)
}
else
{
res <- mla(b=b, m=length(b), fn=loglikhlme,
clustertype=clustertype,.packages=NULL,
epsa=convB,epsb=convL,epsd=convG,partialH=indexHr,
digits=8,print.info=verbose,blinding=FALSE,
multipleTry=25,file="",
nproc=nproc,maxiter=maxiter,minimize=FALSE,
Y0=Y0,X0=X0,prior0=prior0,pprior0=pprior0,idprob0=idprob0,idea0=idea0,
idg0=idg0,idcor0=idcor0,ns0=ns0,ng0=ng0,nv0=nv0,nobs=nobs0,
nea0=nea0,nmes0=nmes0,idiag=idiag0,nwg0=nwg0,ncor0=ncor0,
npm0=NPM,fix0=fix0,nfix0=nfix,bfix0=bfix)
out <- list(conv=res$istop, V=res$v, best=res$b,
ppi2=rep(NA,ns0*ng0), predRE=rep(NA,ns0*nea0),
varRE=rep(NA,ns0*nea0*(nea0+1)/2),
resid_m=rep(NA,nobs0), resid_ss=rep(NA,nobs0),
pred_m_g=rep(NA,nobs0*ng0), pred_ss_g=rep(NA,nobs0*ng0),
gconv=c(res$ca, res$cb, res$rdm), niter=res$ni,
loglik=res$fn.value)
}
if(out$conv %in% c(1,2,3))
{
estim0 <- 0
ll <- 0
ppi0 <- rep(0,ns0*ng0)
resid_m <- rep(0,nobs0)
resid_ss <- rep(0,nobs0)
pred_m_g <- rep(0,nobs0*ng0)
pred_ss_g <- rep(0,nobs0*ng0)
predRE <- rep(0,ns0*nea0)
varRE <- rep(0,ns0*nea0*(nea0+1)/2)
post <- .Fortran(C_loglikhlme,
as.double(Y0),
as.double(X0),
as.integer(prior0),
as.double(pprior0),
as.integer(idprob0),
as.integer(idea0),
as.integer(idg0),
as.integer(idcor0),
as.integer(ns0),
as.integer(ng0),
as.integer(nv0),
as.integer(nobs0),
as.integer(nea0),
as.integer(nmes0),
as.integer(idiag0),
as.integer(nwg0),
as.integer(ncor0),
as.integer(NPM),
as.double(out$best),
ppi=as.double(ppi0),
resid_m=as.double(resid_m),
resid_ss=as.double(resid_ss),
pred_m_g=as.double(pred_m_g),
pred_ss_g=as.double(pred_ss_g),
predRE=as.double(predRE),
varRE=as.double(varRE),
as.integer(fix0),
as.integer(nfix),
as.double(bfix),
as.integer(estim0),
loglik=as.double(ll))
out$ppi2 <- post$ppi
out$predRE <- post$predRE
out$varRE <- post$varRE
out$resid_m <- post$resid_m
out$resid_ss <- post$resid_ss
out$pred_m_g <- post$pred_m_g
out$pred_ss_g <- post$pred_ss_g
}
## creer best a partir de b et bfix
best <- rep(NA,length(fix0))
best[which(fix0==0)] <- out$best
best[which(fix0==1)] <- bfix
names(best) <- names_best
out$best <- best
NPM <- NPM+nfix
### mettre 0 pr les prm fixes
if(length(posfix))
{
if(out$conv==3)
{
mr <- NPM-sum(pbH0)-length(posfix)
Vr <- matrix(0,mr,mr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
Vr <- t(Vr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
V <- matrix(NA,NPM,NPM)
V[setdiff(1:NPM,c(which(pbH0==1),posfix)),setdiff(1:NPM,c(which(pbH0==1),posfix))] <- Vr
V[,posfix] <- 0
V[posfix,] <- 0
V <- V[upper.tri(V,diag=TRUE)]
}
else
{
mr <- NPM-length(posfix)
Vr <- matrix(0,mr,mr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
Vr <- t(Vr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
V <- matrix(0,NPM,NPM)
V[setdiff(1:NPM,posfix),setdiff(1:NPM,posfix)] <- Vr
V <- V[upper.tri(V,diag=TRUE)]
}
}
else
{
if(out$conv==3)
{
mr <- NPM-sum(pbH0)
Vr <- matrix(0,mr,mr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
Vr <- t(Vr)
Vr[upper.tri(Vr,diag=TRUE)] <- out$V[1:(mr*(mr+1)/2)]
V <- matrix(NA,NPM,NPM)
V[setdiff(1:NPM,which(pbH0==1)),setdiff(1:NPM,which(pbH0==1))] <- Vr
V <- V[upper.tri(V,diag=TRUE)]
}
else
{
V <- out$V
}
}
### Creation du vecteur cholesky
Cholesky <- rep(0,(nea0*(nea0+1)/2))
if(idiag0==0 & NVC>0){
Cholesky[1:NVC] <- out$best[(NPROB+NEF+1):(NPROB+NEF+NVC)]
### Construction de la matrice U
U <- matrix(0,nrow=nea0,ncol=nea0)
U[upper.tri(U,diag=TRUE)] <- Cholesky[1:NVC]
z <- t(U) %*% U
out$best[(NPROB+NEF+1):(NPROB+NEF+NVC)] <- z[upper.tri(z,diag=TRUE)]
}
if(idiag0==1 & NVC>0){
id <- 1:nea0
indice <- rep(id+id*(id-1)/2)
Cholesky[indice] <- out$best[(NPROB+NEF+1):(NPROB+NEF+nea0)]
out$best[(NPROB+NEF+1):(NPROB+NEF+NVC)] <- out$best[(NPROB+NEF+1):(NPROB+NEF+NVC)]**2
}
####################################################
if (nea0>0) {
predRE <- matrix(out$predRE,ncol=nea0,byrow=TRUE)
predRE <- data.frame(INDuniq,predRE)
colnames(predRE) <- c(nom.subject,nom.X0[idea0!=0])
}
if(ng0>1) {
ppi<- matrix(out$ppi2,ncol=ng0,byrow=TRUE)
}
else {
ppi <- matrix(rep(1,ns0),ncol=ng0)
}
chooseClass <- function(ppi)
{
res <- which.max(ppi)
if(!length(res)) res <- 0
return(res)
}
classif <- apply(ppi,1,chooseClass)
if(any(!is.finite(ppi)))
{
classif <- rep(NA,ns0)
iok <- which(is.finite(ppi[,1]))
classif[iok] <- apply(ppi[iok,,drop=FALSE],1,which.max)
}
ppi <- data.frame(INDuniq,classif,ppi)
temp <- paste("prob",1:ng0,sep="")
colnames(ppi) <- c(nom.subject,"class",temp)
rownames(ppi) <- 1:ns0
pred_m_g <- matrix(out$pred_m_g,nrow=nobs0)
pred_ss_g <- matrix(out$pred_ss_g,nrow=nobs0)
pred_m <- Y0-out$resid_m
pred_ss <- Y0-out$resid_ss
pred <- data.frame(IND,pred_m,out$resid_m,pred_ss,out$resid_ss,Y0,pred_m_g,pred_ss_g)
temp<-paste("pred_m",1:ng0,sep="")
temp1<-paste("pred_ss",1:ng0,sep="")
colnames(pred)<-c(nom.subject,"pred_m","resid_m","pred_ss","resid_ss","obs",temp,temp1)
if(!is.null(var.time))
{
pred <- data.frame(IND,pred_m,out$resid_m,pred_ss,out$resid_ss,Y0,pred_m_g,pred_ss_g,timeobs)
colnames(pred)<-c(nom.subject,"pred_m","resid_m","pred_ss","resid_ss","obs",temp,temp1,var.time)
}
### ad 2/04/2012
if (!("intercept" %in% nom.X0)) X0.names2 <- X0.names2[-1]
### ad
## levels = modalites des variables dans X0 (si facteurs)
levelsdata <- vector("list", length(X0.names2))
levelsfixed <- vector("list", length(X0.names2))
levelsrandom <- vector("list", length(X0.names2))
levelsmixture <- vector("list", length(X0.names2))
levelsclassmb <- vector("list", length(X0.names2))
names(levelsdata) <- X0.names2
names(levelsfixed) <- X0.names2
names(levelsrandom) <- X0.names2
names(levelsmixture) <- X0.names2
names(levelsclassmb) <- X0.names2
for(v in X0.names2)
{
if(v == "intercept") next
if(is.factor(data[,v]))
{
levelsdata[[v]] <- levels(data[,v])
}
#else
#{
if(length(grep(paste("factor\\(",v,"\\)",sep=""), fixed)))
{
levelsfixed[[v]] <- levels(as.factor(data[,v]))
}
# else
# {
if(length(grep(paste("factor\\(",v,"\\)",sep=""), random)))
{
levelsrandom[[v]] <- levels(as.factor(data[,v]))
}
# else
# {
if(length(grep(paste("factor\\(",v,"\\)",sep=""), mixture)))
{
levelsmixture[[v]] <- levels(as.factor(data[,v]))
}
# else
# {
if(length(grep(paste("factor\\(",v,"\\)",sep=""), classmb)))
{
levelsclassmb[[v]] <- levels(as.factor(data[,v]))
}
# }
# }
# }
#}
}
levels <- list(levelsdata=levelsdata, levelsfixed=levelsfixed, levelsrandom=levelsrandom,
levelsmixture=levelsmixture, levelsclassmb=levelsclassmb)
cost<-proc.time()-ptm
res <-list(ns=ns0,ng=ng0,idea0=idea0,idprob0=idprob0,idg0=idg0,idcor0=idcor0,loglik=out$loglik,best=out$best,V=V,gconv=out$gconv,conv=out$conv,call=cl,niter=out$niter,N=N,idiag=idiag0,pred=pred,pprob=ppi,predRE=predRE,varRE=out$varRE,Xnames=nom.X0,Xnames2=X0.names2,cholesky=Cholesky,na.action=na.action,AIC=2*(length(out$best)-length(posfix)-out$loglik),BIC=(length(out$best)-length(posfix))*log(ns0)-2*out$loglik,data=datareturn,wRandom=wRandom,b0Random=b0Random, levels=levels, var.time=var.time, runtime=cost[3])
class(res) <-c("hlme")
if(verbose==TRUE) cat("The program took", round(cost[3],2), "seconds \n")
res
}
Try the lcmm package in your browser
Any scripts or data that you put into this service are public.
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.