Nothing
#' Permutation of the latent classes
#'
#' This function allows a re-ordering of the latent classes of an estimated model.
#'
#' @param m an object inheriting from classes \code{hlme}, \code{lcmm}, \code{multlcmm} or \code{Jointlcmm}
#' @param order a vector (integer between 1 and ng) containing the new order of the latent classes
#' @param estim optional boolean specifying if the model should estimated with the reordered parameters as initial values. By default, the model is estimated. If FALSE, only the coefficients in \code{$best} are modified. All other outputs are not changed.
#' @return An object of the same class as m, with reordered classes, or the initial object with new coefficients if estim is FALSE.
#' @author Viviane Philipps and Cecile Proust-Lima
#' @examples
#'
#' ## Estimation of a hlme model with 2 classes
#' m2 <- 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))
#'
#' ## Exchange class 2 and class 1
#' m2b <- permut(m2,order=c(2,1))
#'
#' @export
permut <- function(m,order,estim=TRUE)
{
ng <- m$ng
if(ng==1) stop("Please use this function only with latent classes")
if(length(order)!=ng) stop(paste("Argument 'order' should be of length",ng))
if(!all(order %in% 1:ng)) stop(paste("Please specify classes between 1 and",ng,"in argument 'order'",collapse=" "))
if(!all(c(1:ng) %in% order)) stop("Please specify all the classes in argument 'order'")
if(inherits(m, "mpjlcmm")) {return(permutmpj(m,order=order,estim=estim))}
if(!inherits(m,c("hlme","lcmm","multlcmm","Jointlcmm"))) stop("Please use this function only with hlme, lcmm, multlcmm or Jointlcmm models")
bnew <- m$best
## coef de classmb reordonnes:
nv <- m$N[1]/(ng-1) # nb de variables dans classmb (avec l'intercept)
coefold <- rep(0,ng*nv) # coef 0 pour la classe G
coefold[sapply(c(0:(nv-1)*ng),"+",1:(ng-1))] <- m$best[1:m$N[1]] # coef de best pour les autres classes
coeford <- coefold[as.vector(sapply(c((0:(nv-1))*ng),"+",order))] #ordonner selon l'ordre donne dans l'argument order
coefref <- coeford-rep(coeford[sapply(c(0:(nv-1)*ng),"+",ng)],each=ng) # soustraire le coef de la nouvelle classe de reference
coefnew <- coefref[sapply(c(0:(nv-1)*ng),"+",1:(ng-1))] # enlever les coef de la nouvelle classe de reference
bnew[1:m$N[1]] <- coefnew
## coef survie pour Jointlcmm
if(inherits(m,"Jointlcmm"))
{
## risque de base
nbevt <- length(m$hazard[[1]])
nprisq <- rep(NA,nbevt)
nrisq <- rep(NA,nbevt)
typrisq <- m$hazard[[1]]
hazardtype <- m$hazard[[2]]
nz <- m$hazard[[4]]
for(ke in 1:nbevt)
{
if(typrisq[ke]==1) nprisq[ke] <- nz[ke]-1
if(typrisq[ke]==2) nprisq[ke] <- 2
if(typrisq[ke]==3) nprisq[ke] <- nz[ke]+2
if(hazardtype[ke]=="Common") nrisq[ke] <- nprisq[ke]
if(hazardtype[ke]=="PH") nrisq[ke] <- nprisq[ke]+m$ng-1
if(hazardtype[ke]=="Specific") nrisq[ke] <- nprisq[ke]*m$ng
if(hazardtype[ke]=="Specific")
{
## echanger les coef
bnew[m$N[1]+sum(nrisq[1:ke])-nrisq[ke]+1:nrisq[ke]] <- m$best[m$N[1]+sum(nrisq[1:ke])-nrisq[ke]+sapply(order,function(x,np) {(x-1)*np+1:np},np=nprisq[ke])]
}
if(hazardtype[ke]=="PH")
{
if(order[ng]==ng)
{
wref <- 0
}
else
{
wref <- m$best[m$N[1]+sum(nrisq[1:ke])-nrisq[ke]+nprisq[ke]+order[ng]]
}
if(typrisq[ke]==2) #Weibull
{
if(m$logspecif==1)
{
bnew[m$N[1]+sum(nrisq[1:ke])-nrisq[ke]+1] <- m$best[m$N[1]+sum(nrisq[1:ke])-nrisq[ke]+1] + wref # a devient a+wref
}
else
{
bnew[m$N[1]+sum(nrisq[1:ke])-nrisq[ke]+1] <- m$best[m$N[1]+sum(nrisq[1:ke])-nrisq[ke]+1] * exp(wref/(2*m$best[m$N[1]+sum(nrisq[1:ke])-nrisq[ke]+2]^2)) # a devient a*exp(wref/(2*b^2))
}
## b ne change pas
}
else #piecewise ou splines
{
if(m$logspecif==1)
{
bnew[m$N[1]+sum(nrisq[1:ke])-nrisq[ke]+1:nprisq[ke]] <- m$best[m$N[1]+sum(nrisq[1:ke])-nrisq[ke]+1:nprisq[ke]] + wref
}
else
{
bnew[m$N[1]+sum(nrisq[1:ke])-nrisq[ke]+1:nprisq[ke]] <- m$best[m$N[1]+sum(nrisq[1:ke])-nrisq[ke]+1:nprisq[ke]] * exp(wref/2)
}
}
## PH : wg devient wg-wref
coefold <- c(m$best[m$N[1]+sum(nrisq[1:ke])-nrisq[ke]+nprisq[ke]+1:(ng-1)],0) # car wG=0
coeford <- coefold[order] # ordonner selon order
coefref <- coeford-coeford[ng] # soustraire le coef de la nouvelle classe de reference
coefnew <- coefref[1:(ng-1)] # enlever le dernier
bnew[m$N[1]+sum(nrisq[1:ke])-nrisq[ke]+nprisq[ke]+1:(ng-1)] <- coefnew
}
}
## variables explicatives
if(m$N[3]>0)
{
idspecif <- matrix(m$idspecif,nrow=nbevt,byrow=TRUE)
avt <- 0
for(j in 1:ncol(idspecif))
{
if((m$idcom[j]==1) & (idspecif[1,j]==1)) # coef commun a tous les evts et pas en mixture
{
avt <- avt+1
}
if((m$idcom[j]==1) & (idspecif[1,j]==2)) # coef commum et en mixture
{
bnew[m$N[1]+m$N[2]+avt+1:ng] <- m$best[m$N[1]+m$N[2]+avt+order] # echanger les coef
avt <- avt+ng
}
if(m$idcom[j]==0) # coef specifique a chaque evt
{
for(ke in 1:nbevt)
{
if(idspecif[ke,j]==1) # pas en mixture
{
avt <- avt+1
}
if(idspecif[ke,j]==2) # coef en mixture
{
bnew[m$N[1]+m$N[2]+avt+1:ng] <- m$best[m$N[1]+m$N[2]+avt+order] # echanger les coef
avt <- avt+ng
}
}
}
}
}
#if(avt!=m$N[3]) stop("pb varxevt")
} # fin partie survie
## effets fixes modele mixte:
if(inherits(m,c("hlme","lcmm")))
{
avtnef <- m$N[1]
nef <- m$N[2]
ny <- 1
avttr <- sum(m$N[1:4]) # ne servira que pour lcmm
}
if(inherits(m,"multlcmm"))
{
avtnef <- m$N[1]
nef <- m$N[3]-m$N[1]-m$N[2]
ny <- m$N[8]
avttr <- sum(m$N[3:8])
}
if(inherits(m,"Jointlcmm"))
{
avtnef <- m$N[1]+m$N[2]+m$N[3]
nef <- m$N[4]
ny <- 1
avttr <- sum(m$N[1:7])
}
tmpnef <- 0
if(m$idg[1]==2) # intercept en mixture
{
if(inherits(m,"hlme") | (inherits(m,"Jointlcmm") & isTRUE(m$linktype==-1)))
{
bnew[avtnef+1:ng] <- m$best[avtnef+order] # echanger
tmpnef <- tmpnef+ng
}
else
{
coefold <- c(0,m$best[avtnef+1:(ng-1)]) # car intercept fixe a 0 dans la premiere classe
coeford <- coefold[order] # ordonner
newIref <- coeford[1] # servira pour transfo
coefref <- coeford-coeford[1] # soustraire le nouveau coef de ref
coefnew <- coefref[2:ng] # enlever le premier coef
bnew[avtnef+1:(ng-1)] <- coefnew
tmpnef <- tmpnef+ng-1
}
}
if(nef>tmpnef) # autres coef
{
for(k in 2:length(m$idg))
{
if(m$idg[k]==1)
{
tmpnef <- tmpnef+1
}
if(m$idg[k]==2)
{
bnew[avtnef+tmpnef+1:ng] <- m$best[avtnef+tmpnef+order] #echanger
tmpnef <- tmpnef+ng
}
}
}
## coef nw:
if(inherits(m,c("hlme","lcmm")))
{
avtnw <- m$N[1]+m$N[2]+m$N[3]
nw <- m$N[4]
nvc <- m$N[3]
}
if(inherits(m,"multlcmm"))
{
avtnw <- m$N[3]+m$N[4]
nw <- m$N[5]
nvc <- m$N[4]
}
if(inherits(m,"Jointlcmm"))
{
avtnw <- m$N[1]+m$N[2]+m$N[3]+m$N[4]+m$N[5]
nw <- m$N[6]
nvc <- m$N[5]
}
if(nw>0)
{
if(nvc>0)
{
reold <- m$best[(avtnw-nvc+1):avtnw]
renew <- reold*m$best[avtnw+order[ng]]^2
bnew[(avtnw-nvc+1):avtnw] <- renew
}
coefold <- c(m$best[avtnw+1:(ng-1)],1)
coeford <- coefold[order]
coefref <- coeford/coeford[ng]
coefnew <- coefref[1:(ng-1)]
bnew[avtnw+1:(ng-1)] <- coefnew
}
## coef transfo
if(any(m$linktype>-1))
{
if(inherits(m,c("lcmm","Jointlcmm"))) m$nbnodes <- length(m$linknodes)
sumntr <- 0
for(k in 1:ny)
{
if(m$linktype[k]==0) #linear (Y-a)/b, a devient a + b*newIref
{
bnew[avttr+sumntr+1] <- bnew[avttr+sumntr+1] + bnew[avttr+sumntr+2]*newIref
sumntr <- sumntr + 2
}
if(m$linktype[k]==1) #beta (b(Y)-a)/b, a devient a + b*newIref
{
bnew[avttr+sumntr+3] <- bnew[avttr+sumntr+3] + bnew[avttr+sumntr+4]*newIref
sumntr <- sumntr + 4
}
if(m$linktype[k]==2) #splines a+sum(bk*Ik(Y)), a devient a - newIref
{
bnew[avttr+sumntr+1] <- bnew[avttr+sumntr+1] - newIref
sumntr <- sumntr + m$nbnodes[k]+2
}
}
}
## modele avec les nouveaux coef
if(estim==TRUE)
{
z <- m$call
z$B <- bnew
mnew <- eval(z, envir=parent.frame())
}
else
{
mnew <- m
mnew$best <- bnew
}
return(mnew)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.