R/predictlink.lcmm.R

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

#' @export
predictlink.lcmm <- function(x,ndraws=2000,Yvalues,...)
{
    ## verification des arguments
    if(missing(x)) stop("The model should be specified.")
    if(!(inherits(x,"lcmm"))) stop("To use only with \"lcmm\" objects")
    if(!missing(Yvalues) & x$linktype==3) warning("With thresholds links, no \"Yvalues\" should be specified. Default values will be used. \n")
    if(x$conv!=1 & ndraws!=0) stop("No confidence intervals can be produced since the program did not converge properly")
    
    if(x$conv %in% c(1,2,3)) 
        {
            if(missing(Yvalues) | x$linktype==3)
                {
                    new.transf <- FALSE
                    Yvalues <- x$estimlink[,1]
                }
            else
                {
                    new.transf <- TRUE
                    Yvalues <- na.omit(Yvalues)
                    if(any(Yvalues<x$estimlink[1,1]) | any(Yvalues>x$estimlink[nrow(x$estimlink),1])) stop("The values specified in \"Yvalues\" are not in the range of the outcome")
                    Yvalues <- sort(Yvalues)
                }
            
            
            ##preparation des arguments pour Fortran
            npm <- length(x$best)
            best <- x$best
            if(x$idiag==0 & x$N[3]>0) best[(x$N[1]+x$N[2]+1):(x$N[1]+x$N[2]+x$N[3])] <- x$cholesky
            if(x$idiag==1 & x$N[3]>0) best[(x$N[1]+x$N[2]+1):(x$N[1]+x$N[2]+x$N[3])] <- sqrt(best[(x$N[1]+x$N[2]+1):(x$N[1]+x$N[2]+x$N[3])])
            
            if(x$linktype==0) ntrtot <- 2
            if(x$linktype==1) ntrtot <- 4
            if(x$linktype==2) ntrtot <- length(x$linknodes)+2
            if(x$linktype==3) ntrtot <- sum(x$ide)
            
            imoins <- x$N[1]+x$N[2]+x$N[3]+x$N[4]
            zitr <- x$linknodes
            maxnbzitr <- ifelse(x$linktype==2,length(x$linknodes),2)
            epsY <- x$epsY
            minY <- x$estimlink[1,1]
            maxY <- x$estimlink[nrow(x$estimlink),1]
            ny <- 1
            nsim <- length(Yvalues)
            
            if(x$linktype==3)
                {
                    ide <- x$ide
                    dimide <- length(ide)
                }
            else
                {
                    ide <- rep(0,1) #pas utilise si pas ordinal
                    dimide <- 1
                }
            
            ndraws <- as.integer(ndraws)

            posfix <- eval(x$call$posfix)

            if(ndraws>0)
                {
                    Mat <- matrix(0,ncol=npm,nrow=npm)
                    Mat[upper.tri(Mat,diag=TRUE)]<- x$V
                    if(length(posfix))
                        {
                            Mat2 <- Mat[-posfix,-posfix]
                            Chol2 <- chol(Mat2)
                            Chol <- matrix(0,npm,npm)
                            Chol[setdiff(1:npm,posfix),setdiff(1:npm,posfix)] <- Chol2
                            Chol <- t(Chol)
                        }
                    else
                        {
                            Chol <- chol(Mat)
                            Chol <- t(Chol)
                        }
                }
            
            
            
            ## calcul des valeurs trasnformees si necessaire
            if(isTRUE(new.transf) & ndraws==0)
                {
                    resFortran <- rep(0,nsim)
                    
                    out0 <- .Fortran(C_calculustransfo,
                                     as.double(best),
                                     as.integer(npm),
                                     as.integer(ny),
                                     as.integer(x$linktype),
                                     as.integer(ntrtot),
                                     as.integer(imoins),
                                     as.double(zitr),
                                     as.integer(maxnbzitr),
                                     as.double(Yvalues),
                                     as.integer(nsim),
                                     as.double(minY),
                                     as.double(maxY),
                                     as.double(epsY),
                                     as.integer(ide),
                                     as.integer(dimide),
                                     transfo=as.double(resFortran))
                    
                    transfY <- out0$transfo
                }
            else
                {
                    transfY <- x$estimlink[,2]
                }

            if(ndraws>0)
                {
                    if(x$conv==1)
                        {
                            ## boucle pour b=1,...B :
                            Hydraws <- NULL
                            for (j in 1:ndraws)
                                {
                                    bdraw <- rnorm(npm)
                                    bdraw <- best + Chol %*% bdraw
                                    
                                    resFortran <- rep(0,nsim)
                                    
                                    out <- .Fortran(C_calculustransfo,
                                                    as.double(bdraw),
                                                    as.integer(npm),
                                                    as.integer(ny),
                                                    as.integer(x$linktype),
                                                    as.integer(ntrtot),
                                                    as.integer(imoins),
                                                    as.double(zitr),
                                                    as.integer(maxnbzitr),
                                                    as.double(Yvalues),
                                                    as.integer(nsim),
                                                    as.double(minY),
                                                    as.double(maxY),
                                                    as.double(epsY),
                                                    as.integer(ide),
                                                    as.integer(dimide),
                                                    transfo=as.double(resFortran))
                                    
                                    Hydraws <- cbind(Hydraws,out$transfo)
                                }
                            
                            ## calcul des bornes IC
                            f <- function(x)
                                {
                                    quantile(x[!is.na(x)],probs=c(0.5,0.025,0.975))
                                }
                            
                            Hydistr <- apply(Hydraws,1,FUN=f)
                            borne_inf <- as.vector(Hydistr[2,])
                            borne_sup <- as.vector(Hydistr[3,])
                            mediane <- as.vector(Hydistr[1,])
                            
                            
                            ## resultat a renvoyer
                            res <- data.frame(Yvalues=Yvalues,transfY_50=mediane,transfY_2.5=borne_inf,transfY_97.5=borne_sup)
                        }
                    
                    if(x$conv==2|x$conv==3)
                        {
                            borne_inf <- rep(NA,length(Yvalues))
                            borne_sup <- rep(NA,length(Yvalues))
                            mediane <- rep(NA,length(Yvalues))
                            res <- data.frame(Yvalues=Yvalues,transfY_50=mediane,transfY_2.5=borne_inf,transfY_97.5=borne_sup)
                        }
                }
            else
                {
                    res <- data.frame(Yvalues=Yvalues,transY=transfY)
                }
            
        }
    else
        {
            cat("Output can not be produced since the program stopped abnormally.")
            res <- NA 
        }

    res.list <- NULL
    res.list$pred <- res
    res.list$object <- x
    class(res.list) <- "predictlink"
    return(res.list)
}