Nothing
predictMMSE <- function(model,VarTime,Timelim,nTime,Xprofile,methInteg="GH",nsim=20,draws=FALSE,ndraws=2000)
{
if(missing(model)) stop("The argument 'model' should be specified")
if(!inherits(model,"hlme")) stop("Use with hlme object only")
if(!(methInteg %in% c("GH","MC",0,1))) stop("Only Gauss-Hermite (0 or 'GH') and Monte-Carlo (1 or 'MC') integration methods are available")
if(methInteg %in% c("GH",0) & !(nsim %in% c(5,7,9,15,20,30,40,50))) stop("For Gauss-Hermite integration method, 'nsim' should be either 5,7,9,15,20,30,40 or 50")
if(isTRUE(draws) & model$conv!=1) stop("No confidence interval can be computed since the program did not converge properly ")
nodes <- c(0,10,20,23,26,28,30)
prm.transfo <- c(0,10.33921,8.06659,10.82235,17.69599,11.59447,14.93948,13.16887,13.37304)
nbclasses <- model$ng
time <- seq(Timelim[1],Timelim[2],length.out=nTime)
newdata <- as.data.frame(sapply(Xprofile,rep,length.out=nTime))
newdata <- cbind(time,newdata)
colnames(newdata) <- c(VarTime,names(Xprofile))
mu <- lcmm::predictY(model,newdata=newdata,var.time=VarTime,draws=FALSE)[[1]]
nobs <- length(mu[,1])
ch <- matrix(0,sum(model$idea0),sum(model$idea0))
ch[lower.tri(ch,diag=TRUE)] <- model$cholesky
if(is.null(model$call$random)) {call_random <- ~-1} else {call_random <- model$call$random }
z <- all.names(call_random)
ind_factor <- which(z=="factor")
olddata <- eval(model$call$data)
if(length(ind_factor))
{
nom.factor <- z[ind_factor+1]
for (v in nom.factor)
{
mod <- levels(as.factor(olddata[,v]))
if (!all(levels(as.factor(newdata1[,v])) %in% mod)) stop(paste("invalid level in factor", v))
newdata1[,v] <- factor(newdata1[,v], levels=mod)
}
}
call_random <- gsub("factor","",call_random)
Z <- model.matrix(formula(paste("~",call_random,sep="")),data=newdata)
VC <- Z %*% ch %*% t(ch) %*% t(Z) + model$best["stderr"]^2*diag(nobs)
if(length(model$N)>4)
{
if(model$N[5]>0)
{
corChar <- as.character(model$call$cor)
corTime <- newdata[,corChar[2]]
if(corChar[1]=="BM")
{
Corr <- model$best["cor1"]^2 * outer(corTime,corTime,Vectorize(function(x,y) min(abs(x),abs(y))))
}
if(corChar[1]=="AR")
{
Corr <- model$best["cor2"]^2 * (1/outer(corTime,corTime,Vectorize(function(x,y) exp(abs(x-y)))))^model$best["cor1"]
}
VC <- VC + Corr
}
}
nbzitr <- length(nodes)
methInteg <- ifelse(methInteg %in% c("GH",0),0,1)
Ymarg <- rep(0,nobs)
if(methInteg==0)
{
VC0 <- diag(VC)
VC1 <- rep(0,nobs*(nobs+1)/2)
}
else
{
VC0 <- rep(0,nobs)
VC1 <- VC[upper.tri(VC,diag=TRUE)]
}
if(nbclasses==1)
{
if(!isTRUE(draws))
{
res <- .Fortran(C_backtransformation,as.double(mu),as.double(VC0),as.double(VC1),as.integer(nobs),as.double(prm.transfo),
as.integer(nbzitr),as.double(nodes),as.integer(nsim),as.integer(methInteg),y=as.double(Ymarg))
res$y[which(res$y==9999)] <- NA
res.pred <-cbind(time, matrix(res$y,ncol=1))
colnames(res.pred) <- c("time","MMSEpred")
class(res.pred) <- "predMMSE"
return(res.pred)
}
else
{
ndraws <- as.integer(ndraws)
ydraws <- NULL
Mat <- matrix(0,ncol=length(model$best),nrow=length(model$best))
# que la partie sup utilis?e donc OK si rien en bas
Mat[upper.tri(Mat,diag=TRUE)]<- model$V
Chol <- chol(Mat)
Chol <- t(Chol)
doOneDraw <- function()
{
bdraw <- rnorm(length(model$best))
bdraw <- model$best + Chol %*% bdraw
modeldraw <- model
modeldraw$best <- bdraw
mudraw <- lcmm::predictY(modeldraw,newdata=newdata,var.time=VarTime,draws=FALSE)[[1]]
Ymarg <- rep(0,nobs)
out <- .Fortran(C_backtransformation,as.double(mudraw),as.double(VC0),as.double(VC1),as.integer(nobs),as.double(prm.transfo),
as.integer(nbzitr),as.double(nodes),as.integer(nsim),as.integer(methInteg),y=as.double(Ymarg))
out$y[out$y==9999] <- NA
return(out$y)
}
ydraws <- replicate(ndraws,doOneDraw())
ydistr <- apply(ydraws,1,FUN=function(x) quantile(x[!is.na(x)],probs=c(0.025,0.5,0.975)))
Ypred_50 <- matrix(ydistr[2,],ncol=1,byrow=F)
Ypred_2.5 <- matrix(ydistr[1,],ncol=1,byrow=F)
Ypred_97.5 <- matrix(ydistr[3,],ncol=1,byrow=F)
Ypred <- cbind(time,Ypred_50,Ypred_2.5,Ypred_97.5)
colnames(Ypred) <- c("time","MMSEdistr_50%","MMSEdistr_2.5%","MMSEdistr_97.5%")
class(Ypred) <- "predMMSE"
return(Ypred)
}
}
else # ie plusieurs classes
{
if(!isTRUE(draws))
{
res <- NULL
for(g in 1:nbclasses)
{
gmu <- mu[,g]
if(model$N[4]==0 | g==nbclasses) wg <- 1
else wg <- (model$best[paste("varprop class",g,sep=" ")])^2
gVC0 <- wg*VC0
gVC1 <- wg*VC1
gres <- .Fortran(C_backtransformation,as.double(gmu),as.double(gVC0),as.double(gVC1),as.integer(nobs),as.double(prm.transfo),
as.integer(nbzitr),as.double(nodes),as.integer(nsim),as.integer(methInteg),y=as.double(Ymarg))
gres$y[which(gres$y==9999)] <- NA
res <- cbind(res,gres$y)
}
res <- cbind(time,res)
colnames(res) <- c("time",paste("MMSEclass",1:nbclasses,sep=""))
class(res) <- "predMMSE"
return(res)
}
else
{
res <- NULL
ndraws <- as.integer(ndraws)
ydraws <- NULL
Mat <- matrix(0,ncol=length(model$best),nrow=length(model$best))
Mat[upper.tri(Mat,diag=TRUE)]<- model$V
Chol <- chol(Mat)
Chol <- t(Chol)
for (j in 1:ndraws)
{
bdraw <- rnorm(length(model$best))
bdraw <- model$best + Chol %*% bdraw
modeldraw <- model
modeldraw$best <- bdraw
mudraw <- lcmm::predictY(modeldraw,newdata=newdata,var.time=VarTime,draws=FALSE)[[1]]
for(g in 1:nbclasses)
{
gmudraw <- mudraw[,g]
if(model$N[4]==0 | g==nbclasses) wg <- 1
else wg <- (model$best[paste("varprop class",g,sep=" ")])^2
gVC0 <- wg*VC0
gVC1 <- wg*VC1
Ymarg <- rep(0,nobs)
gres <- .Fortran(C_backtransformation,as.double(gmudraw),as.double(gVC0),as.double(gVC1),as.integer(nobs),as.double(prm.transfo),
as.integer(nbzitr),as.double(nodes),as.integer(nsim),as.integer(methInteg),y=as.double(Ymarg))
gres$y[which(gres$y==9999)] <- NA
res <- cbind(res,gres$y)
}
}
classe <- factor(paste("class",rep(1:nbclasses,ndraws),sep=""))
ydistr <- apply(res,1,function(v) unlist(tapply(v,classe,FUN=function(x) quantile(x[!is.na(x)],probs=c(0.025,0.5,0.975)))))
ydistr <- t(ydistr)
ydistr <- ydistr[,order(rep(c(2,1,3),nbclasses))]
ydistr <- cbind(time,ydistr)
colnames(ydistr) <- c("time",paste("MMSEdistr_50_class",1:model$ng,sep=""),paste("MMSEdistr_2.5_class",1:model$ng,sep=""),paste("MMSEdistr_97.5_class",1:model$ng,sep=""))
class(ydistr) <- "predMMSE"
return(ydistr)
}
}
}
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.