conddist.saemix: Estimate conditional mean and variance of individual...

In saemix: Stochastic Approximation Expectation Maximization (SAEM) Algorithm

Estimate conditional mean and variance of individual parameters using the MCMC algorithm

Description

When the parameters of the model have been estimated, we can estimate the individual parameters (psi_i).

Usage

conddist.saemix(saemixObject, nsamp = 1, max.iter = NULL, plot = FALSE, ...)

Arguments

saemixObject	an object returned by the saemix function
nsamp	Number of samples to be drawn in the conditional distribution for each subject. Defaults to 1
max.iter	Maximum number of iterations for the computation of the conditional estimates. Defaults to twice the total number of iterations (see above)
plot	a boolean indicating whether to display convergence plots (defaults to FALSE)
...	optional arguments passed to the plots. Plots will appear if the argument plot is set to TRUE

Details

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:

cond.mean.phi

Conditional mean of the individual distribution of the parameters (obtained as the mean of the samples)

cond.var.phi

Conditional variance of the individual distribution of the parameters (obtained as the mean of the estimated variance of the samples)

cond.shrinkage

Estimate of the shrinkage for the conditional estimates

cond.mean.eta

Conditional mean of the individual distribution of the parameters (obtained as the mean of the samples)

phi.samp

An array with 3 dimensions, giving nsamp samples from the conditional distributions of the individual parameters

phi.samp.var

The estimated individual variances for the sampled parameters phi.samp

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).

Author(s)

Emmanuelle Comets [e-mail]

Audrey Lavenu

Marc Lavielle

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.

See Also

SaemixData,SaemixModel, SaemixObject,saemixControl,saemix

Examples

# 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]))

saemix documentation built on July 9, 2023, 7:43 p.m.