R/summary.lcmm.R

In lcmm: Extended Mixed Models Using Latent Classes and Latent Processes

#' Summary of a \code{hlme}, \code{lcmm}, \code{Jointlcmm}, \code{multlcmm},
#' \code{mpjlcmm}, \code{externSurv}, \code{externX},
#' \code{epoce} or \code{Diffepoce} objects
#' 
#' The function provides a summary of \code{hlme}, \code{lcmm}, \code{multlcmm}
#' and \code{Jointlcmm} estimations, or \code{epoce} and \code{Diffepoce}
#' computations.
#' 
#' 
#' @aliases summary.hlme summary.lcmm summary.Jointlcmm summary.multlcmm
#' summary.epoce summary.Diffepoce summary.mpjlcmm summary.externSurv summary.externX
#' @param object an object inheriting from classes \code{hlme}, \code{lcmm},
#' \code{multlcmm} for fitted latent class mixed-effects, or class
#' \code{Jointlcmm}, \code{mpjlcmm} for a Joint latent class mixed model or \code{epoce} or
#' \code{Diffepoce} for predictive accuracy computations or \code{externSurv}, \code{externX}
#' for secondary regression models.
#' @param \dots further arguments to be passed to or from other methods.  They
#' are ignored in this function.
#' @return For \code{epoce} or \code{Diffepoce} objects, returns NULL. For
#' \code{hlme}, \code{lcmm}, \code{Jointlcmm} or \code{multlcmm} returns also a
#' matrix containing the fixed effect estimates in the longitudinal model,
#' their standard errors, Wald statistics and p-values
#' @author Cecile Proust-Lima, Viviane Philipps, Amadou Diakite and Benoit
#' Liquet
#' @seealso \code{\link{hlme}}, \code{\link{lcmm}}, \code{\link{multlcmm}},
#' \code{\link{Jointlcmm}}, \code{epoce}, \code{Diffepoce}
#' @keywords print
#' 
#' @export
#'
summary.lcmm <- function(object,...)
{
    x <- object
    if (!inherits(x, "lcmm")) stop("use only with \"lcmm\" objects")
    
    if(inherits(x, "externVar")){
      cat("Secondary linear mixed model", "\n")
    } else {
      cat("General latent class mixed model", "\n")
    }
    cat("     fitted by maximum likelihood method", "\n")
    
    cl <- x$call
    cl$B <- NULL
    if(is.data.frame(cl$data))
        {
            cl$data <- NULL
            x$call$data <- NULL    
        }
    cat(" \n")
    dput(cl)
    cat(" \n")

    posfix <- eval(cl$posfix)

    cat("Statistical Model:", "\n")
    cat(paste("     Dataset:", as.character(as.expression(x$call$data))),"\n")
    cat(paste("     Number of subjects:", x$ns),"\n")

    cat(paste("     Number of observations:", x$N[5]),"\n")
    if(length(x$na.action))cat(paste("     Number of observations deleted:",length(x$na.action)),"\n")
    cat(paste("     Number of latent classes:", x$ng), "\n")
    cat(paste("     Number of parameters:", length(x$best))," \n")
    if(length(posfix)) cat(paste("     Number of estimated parameters:", length(x$best)-length(posfix))," \n")
    if (x$linktype==0) {
        ntrtot <- 2
        cat("     Link function: linear"," \n")
    }
    if (x$linktype==1)
        {
            ntrtot <- 4
            cat("     Link function: Standardised Beta CdF"," \n")
        }
    if (x$linktype==2) {
        ntrtot <- length(x$linknodes)+2
        cat("     Link function: Quadratic I-splines with nodes"," \n")
                                        #cat(paste("      ",x$linknodes)," \n")
        cat(     x$linknodes," \n")
    }


    if (x$linktype==3) {
        ntrtot <- sum(x$ide==1)
        cat("     Link function: thresholds"," \n")
    }

    cat(" \n")
    cat("Iteration process:", "\n")

    if(x$conv==1) cat("     Convergence criteria satisfied")
    if(x$conv==2) cat("     Maximum number of iteration reached without convergence")
    if(x$conv==3) cat("     Convergence with restrained Hessian matrix")
    if(x$conv==4|x$conv==12) {
        cat("     The program stopped abnormally. No results can be displayed.\n")
    }
    else{

        cat(" \n")
        if(inherits(x, "externVar")) {
          if(x$varest == "paramBoot"){
            cat("     Proportion of convergence on bootstrap iterations (%)=", x$Mconv, "\n")
          } else {
            cat("     Number of iterations: ", x$niter, "\n")
            cat("     Convergence criteria: parameters=", signif(x$gconv[1],2), "\n")
            cat("                         : likelihood=", signif(x$gconv[2],2), "\n")
            cat("                         : second derivatives=", signif(x$gconv[3],2), "\n")
          }
        } else {
          cat("     Number of iterations: ", x$niter, "\n")
          cat("     Convergence criteria: parameters=", signif(x$gconv[1],2), "\n")
          cat("                         : likelihood=", signif(x$gconv[2],2), "\n")
          cat("                         : second derivatives=", signif(x$gconv[3],2), "\n")
        }
        cat(" \n")
        cat("Goodness-of-fit statistics:", "\n")
        cat(paste("     maximum log-likelihood:", round(x$loglik,2))," \n")
        cat(paste("     AIC:", round(x$AIC,2))," \n")
        cat(paste("     BIC:", round(x$BIC,2))," \n")
        cat(" \n")
        
        ncor <- x$N[6]

        if (x$Ydiscrete==1 & ncor==0){
            cat(paste("     Discrete posterior log-likelihood:", round(x$discrete_loglik,2))," \n")
            cat(paste("     Discrete AIC:", round(-2*x$discrete_loglik+2*(length(x$best)-length(posfix)),2))," \n")
            cat(" \n")
            cat(paste("     Mean discrete AIC per subject:",round((-x$discrete_loglik+length(x$best)-length(posfix))/as.double(x$ns),4))," \n")
            cat(paste("     Mean UACV per subject:",round(x$UACV,4))," \n")
            cat(paste("     Mean discrete LL per subject:",round(x$discrete_loglik/as.double(x$ns),4))," \n")
        }
        cat(" \n")

        cat("Maximum Likelihood Estimates:", "\n")
        cat(" \n")

        NPROB <- x$N[1]
        NEF   <- x$N[2]
        NVC   <- x$N[3]
        NW    <- x$N[4]
        NPM   <- length(x$best)


        ## shorten names if > 20 characters
        names_best <- names(x$best)
        if(any(sapply(names_best, nchar)>20))
        {
            islong <- which(sapply(names_best, nchar)>20)
            split_names_best <- strsplit(names_best, split=":", fixed=TRUE)
            short_names_best <- lapply(split_names_best, gsub, pattern="\\(.*\\)", replacement="(...)")
            new_names <- lapply(short_names_best, paste, collapse=":")
            names_best[islong] <- unlist(new_names)[islong]
            names(x$best) <- names_best
            
            islong <- which(sapply(x$Xnames, nchar)>20)
            if(length(islong))
            {
                x$Xnames[islong] <- sapply(x$Xnames[islong], gsub, pattern="\\(.*\\)", replacement="(...)")
            }
        }
        
        
        se <- rep(NA,NPM)
        if (x$conv==1 | x$conv==3)
            {
                ##recuperation des indices de V
                id <- 1:NPM
                indice <- rep(id*(id+1)/2)
                se <-sqrt(x$V[indice])
                #if (NVC>0) se[(NPROB+NEF+1):(NPROB+NEF+NVC)]<-NA
                wald <- x$best/se
                pwald <- 1-pchisq(wald**2,1)
                coef <- x$best
            }
        else
            {
                se <- NA
                wald <- NA
                pwald <- NA
                coef <- x$best

                sech <- rep(NA,length(coef))
                waldch <- rep(NA,length(coef))
                pwaldch <- rep(NA,length(coef))
            }

        ##prendre abs pour les parametres mis au carre
        if(NW>0) coef[NPROB+NEF+NVC+1:NW] <- abs(coef[NPROB+NEF+NVC+1:NW])
        if(ncor>0) coef[NPROB+NEF+NVC+NW+ntrtot+ncor] <- abs(coef[NPROB+NEF+NVC+NW+ntrtot+ncor])

        
        ow <- options("warn")
        options(warn=-1) # to avoid warnings with conv=3
        ## convertir en character
        if(x$conv!=2)
            {
                coefch <- format(as.numeric(sprintf("%.5f",coef)),nsmall=5,scientific=FALSE)
                sech <- format(as.numeric(sprintf("%.5f",se)),nsmall=5,scientific=FALSE)
                waldch <- format(as.numeric(sprintf("%.3f",wald)),nsmall=3,scientific=FALSE)
                pwaldch <- format(as.numeric(sprintf("%.5f",pwald)),nsmall=5,scientific=FALSE)
            }
        else
            {
                coefch <- format(as.numeric(sprintf("%.5f",coef)),nsmall=5,scientific=FALSE)
            }
        options(ow)
        
        if(length(posfix))
            {
                coefch[posfix] <- paste(coefch[posfix],"*",sep="")
                sech[posfix] <- ""
                waldch[posfix] <- ""
                pwaldch[posfix] <- ""
            }


        ## fct pr determiner la longueur max d'une chaine de caracteres
        ## (avec gestion des NA)
        maxchar <- function(x)
            {
                xx <- na.omit(x)
                if(length(xx))
                    {
                        res <- max(nchar(xx))
                    }
                else
                    {
                        res <- 2
                    }
                return(res)
            }
        
        
        if(NPROB>0)
            {
                cat("Fixed effects in the class-membership model:\n" )
                cat("(the class of reference is the last class) \n")

                tmp <- cbind(coefch[1:NPROB],sech[1:NPROB],waldch[1:NPROB],pwaldch[1:NPROB])
                maxch <- apply(tmp,2,maxchar)
                if(any(c(1:NPROB) %in% posfix)) maxch[1] <- maxch[1]-1
                dimnames(tmp) <- list(names(coef)[1:NPROB],
                                      c(paste(paste(rep(" ",max(maxch[1]-4,0)),collapse=""),"coef",sep=""),
                                        paste(paste(rep(" ",max(maxch[2]-2,0)),collapse=""),"Se",sep=""),
                                        paste(paste(rep(" ",max(maxch[3]-4,0)),collapse=""),"Wald",sep=""),
                                        paste(paste(rep(" ",max(maxch[4]-7,0)),collapse=""),"p-value",sep="")))
                    
                    cat("\n")
                    print(tmp,quote=FALSE,na.print="")
                    cat("\n")
            }



        cat("Fixed effects in the longitudinal model:\n" )
        
        tmp <- matrix(c(paste(c(rep(" ",maxchar(coefch[NPROB+1:NEF])-ifelse(any(c(NPROB+1:NEF) %in% posfix),2,1)),0),collapse=""),"","",""),nrow=1,ncol=4)
        tTable <- matrix(c(0,NA,NA,NA),nrow=1,ncol=4)

        if (NEF>0)
            {
                tmp2 <- cbind(coefch[NPROB+1:NEF],
                              sech[NPROB+1:NEF],
                              waldch[NPROB+1:NEF],
                              pwaldch[NPROB+1:NEF])
                tmp <- rbind(tmp,tmp2)
                tTable <- rbind(tTable,cbind(round(coef[NPROB+1:NEF],5),
                                             round(se[NPROB+1:NEF],5),
                                             round(wald[NPROB+1:NEF],3),
                                             round(pwald[NPROB+1:NEF],5)))
            }
            

        if (x$ng>1)
            {
                interc <- "intercept class1 (not estimated)"
            }
        else
            {
                interc <- "intercept (not estimated)"
            }
            
        if(NEF>0)
            {
                maxch <- apply(tmp,2,maxchar)
                if(any(c(NPROB+1:NEF) %in% posfix)) maxch[1] <- maxch[1]-1

                dimnames(tmp) <- list(c(interc,names(coef)[NPROB+1:NEF]),
                                      c(paste(paste(rep(" ",max(maxch[1]-4,0)),collapse=""),"coef",sep=""),
                                        paste(paste(rep(" ",max(maxch[2]-2,0)),collapse=""),"Se",sep=""),
                                        paste(paste(rep(" ",max(maxch[3]-4,0)),collapse=""),"Wald",sep=""),
                                        paste(paste(rep(" ",max(maxch[4]-7,0)),collapse=""),"p-value",sep="")))
                }
            else
                {
                    dimnames(tmp) <- list(interc, c("coef", "Se", "Wald", "p-value"))
                }

            rownames(tTable) <- rownames(tmp)
            colnames(tTable) <-  c("coef", "Se", "Wald", "p-value")
            
            cat("\n")
            print(tmp,quote=FALSE,na.print="")
            cat("\n")


        if(NVC>0)
            {
                cat("\n")
                cat("Variance-covariance matrix of the random-effects:\n" )
                if(x$idiag==1)
                    {
                        if (NVC>1)
                            {
                                Mat.cov <- diag(coef[(NPROB+NEF+1):(NPROB+NEF+NVC)])
                            }
                        else
                            {
                                Mat.cov <- matrix(coef[(NPROB+NEF+1)],ncol=1)
                            }
                        colnames(Mat.cov) <- x$Xnames[x$idea0==1]
                        rownames(Mat.cov) <- x$Xnames[x$idea0==1]
                        Mat.cov[lower.tri(Mat.cov)] <- 0
                        Mat.cov[upper.tri(Mat.cov)] <- NA
                    }
                if(x$idiag==0)
                    {
                        Mat.cov <- matrix(0,ncol=sum(x$idea0),nrow=sum(x$idea0))
                        colnames(Mat.cov) <- x$Xnames[x$idea0==1]
                        rownames(Mat.cov) <- x$Xnames[x$idea0==1]
                        Mat.cov[upper.tri(Mat.cov,diag=TRUE)] <- coef[(NPROB+NEF+1):(NPROB+NEF+NVC)]
                        Mat.cov <- t(Mat.cov)
                        Mat.cov[upper.tri(Mat.cov)] <- NA
                    }
                    
                if(any(posfix %in% c(NPROB+NEF+1:NVC)))
                    {
                        Mat.cov <- apply(Mat.cov,2,format,digits=5,nsmall=5)
                        Mat.cov[upper.tri(Mat.cov)] <- ""
                        pf <- sort(intersect(c(NPROB+NEF+1:NVC),posfix))
                        p <- matrix(0,sum(x$idea0),sum(x$idea0))
                        if(x$idiag==FALSE) p[upper.tri(p,diag=TRUE)] <- c(NPROB+NEF+1:NVC)
                        if(x$idiag==TRUE & NVC>1) diag(p) <- c(NPROB+NEF+1:NVC)
                        if(x$idiag==TRUE & NVC==1) p <- matrix(c(NPROB+NEF+1),1,1)
                        Mat.cov[which(t(p) %in% pf)] <- paste(Mat.cov[which(t(p) %in% pf)],"*",sep="")
                        print(Mat.cov,quote=FALSE)
                    }
                else
                    {
                        prmatrix(round(Mat.cov,5),na.print="")
                    }
                cat("\n")
            }
        

        std <- NULL
        nom <- NULL
        if(NW>=1)
            {
                nom <- paste("Proportional coefficient class",c(1:(x$ng-1)),sep="")
                std <-cbind(coefch[NPROB+NEF+NVC+1:NW],sech[NPROB+NEF+NVC+1:NW])
            }
        if(ncor==2)
            {
                nom <- c(nom,"AR correlation parameter:","AR standard error:")
                    std <-rbind(std,c(coefch[(NPROB+NEF+NVC+NW+ntrtot+1)],
                                      sech[(NPROB+NEF+NVC+NW+ntrtot+1)]),
                                c(coefch[(NPROB+NEF+NVC+NW+ntrtot+2)],
                                  sech[(NPROB+NEF+NVC+NW+ntrtot+2)]))
            }
        if(ncor==1)
            {
                nom <- c(nom,"BM standard error:")
                std <-rbind(std,c(coefch[(NPROB+NEF+NVC+NW+ntrtot+1)],sech[(NPROB+NEF+NVC+NW+ntrtot+1)]))
            }
        
        if (!is.null(std))
            { 
                rownames(std) <- nom
                maxch <- apply(std,2,maxchar)
                if(NW>0 & any(c(NPROB+NEF+NVC+1:NW) %in% posfix))
                    {
                        maxch[1] <- maxch[1]-1
                    }
                else
                    {
                        if(ncor>0 & any(c(NPROB+NEF+NVC+NW+ntrtot+1:ncor) %in% posfix)) maxch[1] <- maxch[1]-1
                    }
            colnames(std) <- c(paste(paste(rep(" ",max(maxch[1]-4,0)),collapse=""),"coef",sep=""),
                               paste(paste(rep(" ",max(maxch[2]-2,0)),collapse=""),"Se",sep=""))
                print(std,quote=FALSE,na.print="")
                cat("\n")
            }

        
        cat("Residual standard error (not estimated) = 1\n")
        cat("\n")

        
        cat("Parameters of the link function:\n" )
        if (x$linktype==3 & ntrtot != (x$linknodes[2]-x$linknodes[1]))
            {
                temp <- (x$linknodes[1]:(x$linknodes[2]-1))*(1-x$ide)
                cat("(the following levels are not observed in the data: ",temp[temp!=0],"\n")
                cat("so that the number of parameters in the threshold transformation is reduced to",ntrtot,") \n")
            }

        tmp <- cbind(coefch[NPROB+NEF+NVC+NW+1:ntrtot],sech[NPROB+NEF+NVC+NW+1:ntrtot],waldch[NPROB+NEF+NVC+NW+1:ntrtot],pwaldch[NPROB+NEF+NVC+NW+1:ntrtot])
        rownames(tmp) <- names(coef)[NPROB+NEF+NVC+NW+1:ntrtot]
        maxch <- apply(tmp,2,maxchar)
        if(any(c(NPROB+NEF+NVC+NW+1:ntrtot) %in% posfix)) maxch[1] <- maxch[1]-1
        colnames(tmp) <- c(paste(paste(rep(" ",max(maxch[1]-4,0)),collapse=""),"coef",sep=""),
                           paste(paste(rep(" ",max(maxch[2]-2,0)),collapse=""),"Se",sep=""),
                           paste(paste(rep(" ",max(maxch[3]-4,0)),collapse=""),"Wald",sep=""),
                           paste(paste(rep(" ",max(maxch[4]-7,0)),collapse=""),"p-value",sep=""))
        if(inherits(x, "externVar")) colnames(tmp)[2] = paste(paste(rep(" ",max(maxch[2]-4,0)),collapse=""),"Se**",sep="")
        cat("\n")
        print(tmp,quote=FALSE,na.print="")
        cat("\n")

        if(length(posfix))
            {
                cat(" * coefficient fixed by the user \n \n")
            }
        if(inherits(x, "externVar")){
          if(x$varest == "none") cat(" ** total variance estimated witout correction for primary model uncertainty", "\n \n")
          if(x$varest == "Hessian") cat(" ** total variance estimated through the Hessian of the joint likelihood", "\n \n")
          if(x$varest == "paramBoot") cat(" ** total variance estimated through parametric bootstrap", "\n \n")
        }
        
        return(invisible(tTable))
    }
}