Nothing
####################### Conditional means - estimates of PSI_i ########################
#' Estimate conditional mean and variance of individual parameters using the
#' MCMC algorithm
#'
#' When the parameters of the model have been estimated, we can estimate the
#' individual parameters (psi_i).
#'
#' Let hattheta be the estimated value of theta computed with the SAEM
#' algorithm and let p(phi_i |y_i; hattheta) be the conditional distribution of
#' phi_i for 1<=i<=N . We use the MCMC procedure used in the SAEM algorithm to
#' estimate these conditional distributions. We empirically estimate the
#' conditional mean E(phi_i |y_i; hattheta) and the conditional standard
#' deviation sd(phi_i |y_i; hattheta).
#'
#' See PDF documentation for details of the computation. Briefly, the MCMC
#' algorithm is used to obtain samples from the individual conditional
#' distributions of the parameters. The algorithm is initialised for each
#' subject to the conditional estimate of the individual parameters obtained at
#' the end of the SAEMIX fit. A convergence criterion is used to ensure
#' convergence of the mean and variance of the conditional distributions. When
#' nsamp>1, several chains of the MCMC algorithm are run in parallel to obtain
#' samples from the conditional distributions, and the convergence criterion
#' must be achieved for all chains. When nsamp>1, the estimate of the
#' conditional mean is obtained by averaging over the different samples, and the
#' samples from the conditional distribution are output as an array of dimension
#' N x nb of parameters x nsamp in arrays phi.samp for the sampled phi_i and psi.samp
#' for the corresponding psi_i. The variance of the conditional phi_i for each
#' sample is given in a corresponding array phi.samp.var (the variance of the
#' conditional psi_i is not given but may be computed via the delta-method or by
#' transforming the confidence interval for phi_i).
#'
#' The shrinkage for any given parameter for the conditional estimate is
#' obtained as
#'
#' Sh=1-var(eta_i)/omega(eta)
#'
#' where var(eta_i) is the empirical variance of the estimates of the
#' individual random effects, and omega(eta) is the estimated variance.
#'
#' When the plot argument is set to TRUE, convergence graphs are produced
#' They can be used assess whether the mean of the individual estimates (on the scale of the parameters)
#' and the mean of the SD of the random effects have stabilised over the ipar.lmcmc (defaults to 50) iterations.
#' The evolution of these variables for each parameter is shown as a continuous line
#' while the shaded area represents the acceptable variation.
#'
#' The function adds or modifies the following elements in the results:
#' \describe{
#' \item{cond.mean.phi}{Conditional mean of the individual
#' distribution of the parameters (obtained as the mean of the samples)}
#' \item{cond.var.phi}{Conditional variance of the individual distribution of
#' the parameters (obtained as the mean of the estimated variance of the samples)}
#' \item{cond.shrinkage}{Estimate of the shrinkage for the conditional estimates}
#' \item{cond.mean.eta}{Conditional mean of the individual distribution of the
#' parameters (obtained as the mean of the samples)}
#' \item{phi.samp}{An array with 3 dimensions, giving nsamp samples
#' from the conditional distributions of the individual parameters}
#' \item{phi.samp.var}{The estimated individual variances for the sampled
#' parameters phi.samp}
#' }
#'
#' @details A warning is output if the maximum number of iterations is reached
#' without convergence (the maximum number of iterations is the sum of the elements in
#' saemix.options$nbiter.saemix).
#'
#' @param saemixObject an object returned by the \code{\link{saemix}} function
#' @param nsamp Number of samples to be drawn in the conditional distribution
#' for each subject. Defaults to 1
#' @param max.iter Maximum number of iterations for the computation of the
#' conditional estimates. Defaults to twice the total number of iterations
#' (see above)
#' @param plot a boolean indicating whether to display convergence plots (defaults to FALSE)
#' @param \dots optional arguments passed to the plots. Plots will appear if
#' the argument plot is set to TRUE
#'
#' @seealso \code{\link{SaemixData}},\code{\link{SaemixModel}},
#' \code{\link{SaemixObject}},\code{\link{saemixControl}},\code{\link{saemix}}
#'
#' @references E Comets, A Lavenu, M Lavielle M (2017). Parameter estimation in nonlinear mixed effect models using saemix,
#' an R implementation of the SAEM algorithm. Journal of Statistical Software, 80(3):1-41.
#'
#' E Kuhn, M Lavielle (2005). Maximum likelihood estimation in nonlinear mixed effects models.
#' Computational Statistics and Data Analysis, 49(4):1020-1038.
#'
#' E Comets, A Lavenu, M Lavielle (2011). SAEMIX, an R version of the SAEM algorithm. 20th meeting of the
#' Population Approach Group in Europe, Athens, Greece, Abstr 2173.
#'
#' @author Emmanuelle Comets \email{emmanuelle.comets@@inserm.fr}
#' @author Audrey Lavenu
#' @author Marc Lavielle
#'
#' @keywords model
#' @examples
#'
#' \donttest{
#' # Not run (strict time constraints for CRAN)
#' data(theo.saemix)
#'
#' saemix.data<-saemixData(name.data=theo.saemix,header=TRUE,sep=" ",na=NA,
#' name.group=c("Id"),name.predictors=c("Dose","Time"),
#' name.response=c("Concentration"),name.covariates=c("Weight","Sex"),
#' units=list(x="hr",y="mg/L",covariates=c("kg","-")), name.X="Time")
#'
#' model1cpt<-function(psi,id,xidep) {
#' dose<-xidep[,1]
#' tim<-xidep[,2]
#' ka<-psi[id,1]
#' V<-psi[id,2]
#' CL<-psi[id,3]
#' k<-CL/V
#' ypred<-dose*ka/(V*(ka-k))*(exp(-k*tim)-exp(-ka*tim))
#' return(ypred)
#' }
#'
#' saemix.model<-saemixModel(model=model1cpt,
#' description="One-compartment model with first-order absorption",
#' psi0=matrix(c(1.,20,0.5,0.1,0,-0.01),ncol=3, byrow=TRUE,dimnames=list(NULL,
#' c("ka","V","CL"))),transform.par=c(1,1,1),
#' covariate.model=matrix(c(0,1,0,0,0,0),ncol=3,byrow=TRUE),fixed.estim=c(1,1,1),
#' covariance.model=matrix(c(1,0,0,0,1,0,0,0,1),ncol=3,byrow=TRUE),
#' omega.init=matrix(c(1,0,0,0,1,0,0,0,1),ncol=3,byrow=TRUE), error.model="constant")
#'
#' saemix.options<-list(seed=632545,save=FALSE,save.graphs=FALSE)
#'
#' saemix.fit<-saemix(saemix.model,saemix.data,saemix.options)
#' saemix.fit<-conddist.saemix(saemix.fit,nsamp=3, plot=TRUE)
#'
#' # First sample from the conditional distribution
#' # (a N (nb of subject) by nb.etas (nb of parameters) matrix)
#' print(head(saemix.fit["results"]["phi.samp"][,,1]))
#'
#' # Second sample
#' print(head(saemix.fit["results"]["phi.samp"][,,2]))
#' }
#'
#' @export conddist.saemix
conddist.saemix<-function(saemixObject, nsamp=1, max.iter=NULL, plot=FALSE, ...) {
# Estimate conditional means and estimates for the individual parameters PSI_i using the MCMC algorithm
# nsamp= number of MCMC samples
# max.iter= max nb of iterations
# returns an array
displayPlot<-plot
saemix.data<-saemixObject["data"]
saemix.model<-saemixObject["model"]
N<-saemix.data["N"]
nb.parameters<-saemixObject["model"]["nb.parameters"]
if(is.null(max.iter)) kmax<-sum(saemixObject["options"]$nbiter.saemix)*2 else kmax<-max.iter
# using several Markov chains
chdat<-new(Class="SaemixRepData",data=saemixObject["data"], nb.chains=nsamp)
NM<-chdat["NM"]
IdM<-chdat["dataM"]$IdM
yM<-chdat["dataM"]$yM
XM<-chdat["dataM"][,c(saemix.data["name.predictors"],saemix.data["name.cens"],saemix.data["name.mdv"],saemix.data["name.ytype"]),drop=FALSE]
io<-matrix(data=0,nrow=N,ncol=max(saemixObject["data"]["nind.obs"]))
for(i in 1:N)
io[i,1:saemixObject["data"]["nind.obs"][i]]<-1
ioM<-do.call(rbind,rep(list(io),nsamp))
ind.ioM <- which(t(ioM)!=0)
DYF<-matrix(data=0,nrow=dim(ioM)[2],ncol=dim(ioM)[1])
ind.eta<-saemixObject["model"]["indx.omega"]
nb.etas<-length(ind.eta)
omega.eta<-saemixObject["results"]["omega"][ind.eta,ind.eta]
diag.omega<-mydiag(saemixObject["results"]["omega"])
domega2<-do.call(cbind,rep(list((sqrt(mydiag(omega.eta)))* saemixObject["options"]$rw.ini), nb.etas))
VK<-rep(c(1:nb.etas),2)
phi<-array(data=0,dim=c(N,nb.parameters, nsamp))
domega<-cutoff(mydiag(saemixObject["results"]["omega"][ind.eta,ind.eta]), .Machine$double.eps)
omega.eta<-saemixObject["results"]["omega"][ind.eta,ind.eta]
omega.eta<-omega.eta-mydiag(mydiag(saemixObject["results"]["omega"][ind.eta, ind.eta]))+mydiag(domega)
chol.omega<-chol(omega.eta)
pres<-saemixObject["results"]["respar"]
args<-list(ind.ioM=ind.ioM)
# Preparing plots
if(displayPlot) {
plot.opt<-saemixObject["prefs"]
plot.opt$xlab<-"Iteration"
plot.opt<-replace.plot.options(plot.opt,...)
change.ylab<-FALSE
if(plot.opt$ylab!=saemixObject["prefs"]$ylab & length(plot.opt$which.par)==1) change.ylab<-TRUE
change.main<-FALSE
if(plot.opt$main!=saemixObject["prefs"]$main & length(plot.opt$which.par)==1) change.main<-TRUE
np<-nb.etas
if(length(plot.opt$mfrow)==0) {
n1<-round(sqrt(np))
n2<-ceiling(np/n1)
if(n1>5 | n2>5) {
n1<-3
n2<-4
cat("Changing the plot layout\n")
}
plot.opt$mfrow<-c(n1,n2)
}
}
# Simulation MCMC
# initialisation a phiM=estimation des parametres individuels
mean.phiM<-do.call(rbind,rep(list(saemixObject["results"]["mean.phi"]),nsamp))
phiM<-do.call(rbind,rep(list(saemixObject["results"]["cond.mean.phi"]),nsamp))
etaM<-phiM[,ind.eta,drop=FALSE]-mean.phiM[,ind.eta,drop=FALSE]
Dargs<-list(transform.par=saemixObject["model"]["transform.par"], structural.model=saemixObject["model"]["model"],IdM=IdM,XM=XM,yM=yM,etype.exp=which(saemixObject["model"]["error.model"] == "exponential"), modeltype=saemixObject["model"]["modeltype"])
Uargs<-list(ind.ioM=ind.ioM)
somega<-solve(omega.eta)
# chol.omega<-chol.omega*3
# domega2<-domega2*3
# somega<-somega*3
U.y<-compute.LLy(phiM,Uargs,Dargs,DYF,pres)
phiMc<-phiM
econd<-sdcond<-array(0,dim=c(nb.parameters,N*kmax,nsamp)) # parametres individuels
ebar<-sdbar<-array(0,dim=c(nb.parameters,kmax, nsamp)) # moyenne des parametres individuels
if(saemixObject["options"]$warnings) cat("Estimating the conditional mean and variance of the distribution of individual parameters\n")
k<-0
while(k<kmax) { # Set a maximum nb of iterations
k<-k+1
for(u in 1:saemixObject["options"]$nbiter.mcmc[1]) { # 1er noyau
etaMc<-matrix(rnorm(NM*nb.etas),ncol=nb.etas)%*%chol.omega
phiMc[,ind.eta]<-mean.phiM[,ind.eta]+etaMc
Uc.y<-compute.LLy(phiMc,Uargs,Dargs,DYF,pres)
deltau<-Uc.y-U.y
ind<-which(deltau<(-1)*log(runif(NM)))
etaM[ind,]<-etaMc[ind,]
U.y[ind]<-Uc.y[ind]
}
U.eta<-0.5*rowSums(etaM*(etaM%*%somega))
# Second stage
if(saemixObject["options"]$nbiter.mcmc[2]>0) {
nt2<-nbc2<-matrix(data=0,nrow=nb.etas,ncol=1)
nrs2<-1
for (u in 1:saemixObject["options"]$nbiter.mcmc[2]) {
for(vk2 in 1:nb.etas) {
etaMc<-etaM
etaMc[,vk2]<-etaM[,vk2]+matrix(rnorm(NM*nrs2), ncol=nrs2)%*%mydiag(domega2[vk2,nrs2],nrow=1) # 2e noyau ? ou 1er noyau+permutation?
phiMc[,ind.eta]<-mean.phiM[,ind.eta]+etaMc
Uc.y<-compute.LLy(phiMc,Uargs,Dargs,DYF,pres)
Uc.eta<-0.5*rowSums(etaMc*(etaMc%*%somega))
deltu<-Uc.y-U.y+Uc.eta-U.eta
ind<-which(deltu<(-1)*log(runif(NM)))
etaM[ind,]<-etaMc[ind,]
U.y[ind]<-Uc.y[ind] # Warning: Uc.y, Uc.eta = vecteurs
U.eta[ind]<-Uc.eta[ind]
nbc2[vk2]<-nbc2[vk2]+length(ind)
nt2[vk2]<-nt2[vk2]+NM
}
}
domega2[,nrs2]<-domega2[,nrs2]*(1+saemixObject["options"]$stepsize.rw* (nbc2/nt2-saemixObject["options"]$proba.mcmc))
}
if(saemixObject["options"]$nbiter.mcmc[3]>0) {
nt2<-nbc2<-matrix(data=0,nrow=nb.etas,ncol=1)
nrs2<-k%%(nb.etas-1)+2
if(is.nan(nrs2)) nrs2<-1 # to deal with case nb.etas=1 ## Maybe need to skip this step if only one IIV ?
for (u in 1:saemixObject["options"]$nbiter.mcmc[3]) {
if(nrs2<nb.etas) {
vk<-c(0,sample(c(1:(nb.etas-1)),nrs2-1))
nb.iter2<-nb.etas
} else {
vk<-0:(nb.etas-1)
nb.iter2<-1
}
for(k2 in 1:nb.iter2) {
vk2<-VK[k2+vk]
etaMc<-etaM
etaMc[,vk2]<-etaM[,vk2]+matrix(rnorm(NM*nrs2), ncol=nrs2)%*%mydiag(domega2[vk2,nrs2])
phiMc[,ind.eta]<-mean.phiM[,ind.eta]+etaMc
Uc.y<-compute.LLy(phiMc,Uargs,Dargs,DYF,pres)
Uc.eta<-0.5*rowSums(etaMc*(etaMc%*%somega))
deltu<-Uc.y-U.y+Uc.eta-U.eta
ind<-which(deltu<(-log(runif(NM))))
etaM[ind,]<-etaMc[ind,]
U.y[ind]<-Uc.y[ind] # Warning: Uc.y, Uc.eta = vecteurs
U.eta[ind]<-Uc.eta[ind]
nbc2[vk2]<-nbc2[vk2]+length(ind)
nt2[vk2]<-nt2[vk2]+NM
}
}
domega2[,nrs2]<-domega2[,nrs2]*(1+saemixObject["options"]$stepsize.rw* (nbc2/nt2-saemixObject["options"]$proba.mcmc))
}
phiM[,ind.eta]<-mean.phiM[,ind.eta]+etaM
phik<-array(t(phiM),dim=c(nb.parameters,N,nsamp))
if(k==1) {
eik<-phik
varik<-0*eik
ebar[,k,]<-apply(eik,c(1,3),mean)
ij<-0
econd[,(ij+1):(ij+N),]<-eik
sdcond[,(ij+1):(ij+N),]<-varik
} else {
eik1<-eik
eik<-eik*(k-1)/k+phik/k
varik<-varik*(k-1)/k+(phik**2)/k + (eik1**2)*(k-1)/k-(eik**2)
sdik<-sqrt(varik)
sdik[sdik<0]<-0
ij<-(k-1)*N
econd[,(ij+1):(ij+N),]<-eik
sdcond[,(ij+1):(ij+N),]<-sdik
ebar[,k,]<-apply(eik,c(1,3),mean)
sdbar[,k,]<-apply(sdik,c(1,3),mean)
if(k>=saemixObject["options"]$ipar.lmcmc) {
ibeg<-max(1,k-saemixObject["options"]$ipar.lmcmc)
ekmax<-(ebar[,k,]+abs(ebar[,k,])* saemixObject["options"]$ipar.rmcmc)
ekmin<-(ebar[,k,]-abs(ebar[,k,])*saemixObject["options"]$ipar.rmcmc)
sdkmax<-(sdbar[,k,]+sdbar[,k,]* saemixObject["options"]$ipar.rmcmc)
sdkmin<-(sdbar[,k,]-sdbar[,k,]*saemixObject["options"]$ipar.rmcmc)
if(nsamp==1) {
maxek<-apply(ebar[,ibeg:k,],1,max)
minek<-apply(ebar[,ibeg:k,],1,min)
iek<-sum(maxek>ekmax)+sum(minek<ekmin)
maxsdk<-apply(sdbar[,ibeg:k,],1,max)
minsdk<-apply(sdbar[,ibeg:k,],1,min)
sdek<-sum(maxsdk>sdkmax)+ sum(minsdk<sdkmin)
} else {
maxek<-apply(ebar[,ibeg:k,],c(1,3),max)
minek<-apply(ebar[,ibeg:k,],c(1,3),min)
iek<-sum(maxek>ekmax)+sum(minek<ekmin)
ibeg2<-ifelse(ibeg==1,2,ibeg)
maxsdk<-apply(sdbar[,ibeg2:k,],c(1,3),max)
minsdk<-apply(sdbar[,ibeg2:k,],c(1,3),min)
sdek<-sum(maxsdk>sdkmax)+ sum(minsdk<sdkmin)
}
if(is.na(sdek)) sdek<-1
if(iek==0 & sdek==0) break
}
}
if(k>50 & k%%50==1) {
if(saemixObject["options"]$warnings) cat(".")
# if(displayPlot) {
# # par(mfrow=plot.opt$mfrow,ask=plot.opt$ask)
# par(mfrow=plot.opt$mfrow)
# limy<-cbind(apply(cbind(apply(ebar[,1:k,],1,min),ekmin),1,min), apply(cbind(apply(ebar[,1:k,],1,max),ekmax),1,max))
# for(ipar in ind.eta) {
# laby<-saemixObject["model"]["name.modpar"][ipar]
# plot(1:k,ebar[ipar,1:k,1],type="n",xlab=plot.opt$xlab, ylab=laby,ylim=limy[ipar,],main=plot.opt$main)
# for(isamp in 1:nsamp)
# lines(1:k,ebar[ipar,1:k,isamp],col=plot.opt$col,lty=plot.opt$lty, lwd=plot.opt$lwd)
# abline(h=apply(ebar[,k,],1,mean), col=plot.opt$ablinecol, lty=plot.opt$ablinelty,lwd=plot.opt$ablinelwd)
# abline(h=apply(ekmax,1,mean)[ipar], col=plot.opt$ablinecol, lty=plot.opt$ablinelty,lwd=plot.opt$ablinelwd)
# abline(h=apply(ekmin,1,mean)[ipar], col=plot.opt$ablinecol, lty=plot.opt$ablinelty,lwd=plot.opt$ablinelwd)
# }
# }
}
}
# cat("k=",k,"\n")
if(saemixObject["options"]$warnings) cat("\n")
if(saemixObject["options"]$warnings) {
if(k>=kmax){cat("Computing the empirical conditional mean and variance: maximum number of iterations reached without meeting convergence criterion (max.iter=",kmax,")\n") }
else cat("Convergence achieved in",k,"iterations\n")
}
if(displayPlot) {
ibeg<-1
nitr<-(k-ibeg)
if(FALSE) { # phi-scale
ebarbar<-apply(ebar[,ibeg:(k-1),],c(1,2),mean) # mean across the samples => npar x nitr
sdbarbar<-apply(sdbar[,ibeg:(k-1),],c(1,2),mean) # mean across the samples => npar x nitr
# sdbarbar<-apply(sdbar,3,mean)
if(displayPlot) {
ypl1<-data.frame(Iteration=rep((ibeg:(k-1)),nb.parameters),condMean=c(t(ebarbar)),par=rep(saemixObject["model"]["name.modpar"],each=nitr),ymin=rep(ekmin,each=nitr),ymax=rep(ekmax,each=nitr),type="mean")
ypl2<-data.frame(Iteration=rep((ibeg:(k-1)),nb.parameters),condMean=c(t(sdbarbar)),par=rep(saemixObject["model"]["name.modpar"],each=nitr),ymin=rep(sdkmin,each=nitr),ymax=rep(sdkmax,each=nitr),type="SD")
ypl<-rbind(ypl1,ypl2)
gpl<-ggplot(data=ypl,aes(x=.data$Iteration,y=.data$condMean)) +geom_line() + facet_wrap(as.factor(.data$type)~as.factor(.data$par), scales="free") + geom_ribbon(aes(ymin=.data$ymin,ymax=.data$ymax),alpha=0.2,fill="darkcyan")
print(gpl)
}
} else { # psi-scale
nitr<-(k-ibeg+1)
for(i in 1:nsamp) {
for(itr in ibeg:k) {
# sdbar[,itr,i]<-dtransphi(t(ebar[,itr,i]),saemixObject["model"]["transform.par"])*sdbar[,itr,i]
ebar[,itr,i]<-transphi(t(ebar[,itr,i]),saemixObject["model"]["transform.par"])
}
}
# tsdkmin<-rowMeans(sdkmin)*dtransphi(t(rowMeans(ekmin)),saemixObject["model"]["transform.par"])
# tsdkmax<-rowMeans(sdkmax)*dtransphi(t(rowMeans(ekmax)),saemixObject["model"]["transform.par"])
if(nsamp==1) {
tsdkmin<-sdkmin
tsdkmax<-sdkmax
tekmin<-transphi(t(ekmin),saemixObject["model"]["transform.par"])
tekmax<-transphi(t(ekmax),saemixObject["model"]["transform.par"])
} else {
tsdkmin<-rowMeans(sdkmin)
tsdkmax<-rowMeans(sdkmax)
tekmin<-transphi(t(rowMeans(ekmin)),saemixObject["model"]["transform.par"])
tekmax<-transphi(t(rowMeans(ekmax)),saemixObject["model"]["transform.par"])
}
ebarbar<-apply(ebar[,ibeg:k,],c(1,2),mean) # mean across the samples => npar x nitr
sdbarbar<-apply(sdbar[,ibeg:k,],c(1,2),mean) # mean across the samples => npar x nitr
# sdbarbar<-apply(sdbar,3,mean)
if(displayPlot) {
ypl1<-data.frame(Iteration=rep(c(ibeg:k),nb.parameters),condMean=c(t(ebarbar)),par=rep(saemixObject["model"]["name.modpar"],each=nitr),ymin=rep(tekmin,each=nitr),ymax=rep(tekmax,each=nitr),type="mean")
ypl2<-data.frame(Iteration=rep(c(ibeg:k),nb.parameters),condMean=c(t(sdbarbar)),par=rep(saemixObject["model"]["name.modpar"],each=nitr),ymin=rep(tsdkmin,each=nitr),ymax=rep(tsdkmax,each=nitr),type="SD")
ypl<-rbind(ypl1,ypl2)
gpl<-ggplot(data=ypl,aes(x=.data$Iteration,y=.data$condMean)) +geom_line() + facet_wrap(as.factor(.data$type)~as.factor(.data$par), scales="free") + geom_ribbon(aes(ymin=.data$ymin,ymax=.data$ymax),alpha=0.2,fill="darkcyan")
print(gpl)
}
}
}
eta.cond<-matrix(0,nrow=dim(etaM)[1],ncol=nb.parameters)
eta.cond[,ind.eta]<-etaM
eta.cond<-array(t(eta.cond),dim=c(nb.parameters,N,nsamp))
eta.cond<-t(apply(eta.cond,c(1,2),mean))
colnames(eta.cond)<-paste("eta.",saemixObject["model"]["name.modpar"],sep="")
cond.shrinkage<-100*(1-apply(eta.cond,2, var)/mydiag(saemixObject["results"]["omega"]))
names(cond.shrinkage)<-paste("Sh.",saemixObject["model"]["name.modpar"],".%",sep="")
resh.eik<-resh.varik<-array(0,dim=c(nrow=N,ncol=nb.parameters,nsamp))
phi1<-array(t(phiM),dim=c(nb.parameters,N,nsamp))
phi.samp<-psi.samp<-array(0,dim=c(N,nb.parameters,nsamp))
for(isamp in 1:nsamp) {
resh.eik[,,isamp]<-t(eik[,,isamp])
resh.varik[,,isamp]<-t(varik[,,isamp])
phi.samp[,,isamp]<-t(phi1[,,isamp])
for(isuj in 1:N)
psi.samp[isuj,,isamp]<-transphi(t(phi.samp[isuj,,isamp]),saemixObject["model"]["transform.par"])
}
cond.mean.psi<-cond.mean.phi<-apply(resh.eik,c(1,2),mean) # Mean across samples of means of phiM at last iteration
for(isuj in 1:N)
cond.mean.psi[isuj,]<-transphi(t(cond.mean.psi[isuj,]),saemixObject["model"]["transform.par"])
saemixObject["results"]["cond.mean.phi"]<-cond.mean.phi
saemixObject["results"]["cond.mean.psi"]<-cond.mean.psi
saemixObject["results"]["cond.var.phi"]<-apply(resh.varik,c(1,2),mean)
saemixObject["results"]["cond.mean.eta"]<-eta.cond
saemixObject["results"]["cond.shrinkage"]<-cond.shrinkage
saemixObject["results"]["phi.samp"]<-phi.samp
saemixObject["results"]["phi.samp.var"]<-resh.varik
saemixObject["results"]["psi.samp"]<-psi.samp
# ECO TODO: verifier ce qui se passe quand un parametre n'est pas estime => faut-il prendre la moyenne des valeurs ou grace a chol.omega qui vaut 0, ca marche ? Verifier d'ailleurs si on a chol.omega quand il manque un parametre (normalement oui)
return(saemixObject)
}
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.