#' @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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.