saemix.bootstrap: Bootstrap for saemix fits

In saemix: Stochastic Approximation Expectation Maximization (SAEM) Algorithm

Bootstrap for saemix fits

Description

This function provides bootstrap estimates for a saemixObject run. Different bootstrap approaches have been implemented (see below for details), with the default method being the conditional non-parametric bootstrap ()

Usage

saemix.bootstrap(
  saemixObject,
  method = "conditional",
  nboot = 200,
  nsamp = 100,
  saemix.options = NULL
)

Arguments

saemixObject	an object returned by the saemix function
method	the name of the bootstrap algorithm to use (one of: "case", "residual", "parametric" or "conditional") (defaults to "conditional")
nboot	number of bootstrap samples
nsamp	number of samples from the conditional distribution (for method="conditional")
saemix.options	list of options to run the saemix algorithm. Defaults to the options in the object, but suppressing the estimation of individual parameters (map=FALSE) and likelihood (ll.is=FALSE), and the options to print out intermediate results and graphs (displayProgress=FALSE,save.graphs=FALSE,print=FALSE)

Details

Different bootstrap algorithms have been proposed for non-linear mixed effect models, to account for the hierarchical nature of these models (see review in Thai et al. 2013, 2014) In saemix, we implemented the following bootstrap algorithms, which can be selected using the "method" argument:

• case refers to case bootstrap, where resampling is performed at the level of the individual and the bootstrap sample consists in sampling with replacement from the observed individuals

• residual refers to non-parametric residual bootstrap, where for each individual the residuals for the random effects and for the residual error are resampled from the individual estimates after the fit

• parametric refers to parametric residual bootstrap, where these residuals are sampled from their theoretical distributions. In saemix, random effects are drawn from a normal distribution using the estimated variance in saemixObject, and transformed to the distribution specified by the transform.par component of the model, and the residual error is drawn from a normal distribution using the estimated residual variance.

• conditional refers to the conditional non-parametric bootstrap, where the random effects are resampled from the individual conditional distributions as in (Comets et al. 2021) and the residual errors resampled from the corresponding residuals

Important note: for discrete data models, all residual-based bootstraps (residual, parametric and conditional) need a simulate.function slot to be included in the model object, as the algorithm will need to generate predictions with the resampled individual parameters in order to generate bootstrap samples. See SaemixModel and the examples provided as notebooks for details.

Author(s)

Emmanuelle Comets [e-mail]

References

Thai H, Mentré F, Holford NH, Veyrat-Follet C, Comets E. A comparison of bootstrap approaches for estimating uncertainty of parameters in linear mixed-effects models. Pharmaceutical Statistics, 2013 ;12:129–40.

Thai H, Mentré F, Holford NH, Veyrat-Follet C, Comets E. Evaluation of bootstrap methods for estimating uncertainty of parameters in nonlinear mixed-effects models : a simulation study in population pharmacokinetics. Journal of Pharmacokinetics and Pharmacodynamics, 2014; 41:15–33.

Comets E, Rodrigues C, Jullien V, Moreno U. Conditional non-parametric bootstrap for non-linear mixed effect models. Pharmaceutical Research, 2021; 38, 1057–66.

Examples

# Bootstrap for the theophylline data 
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(algorithm=c(1,0,0),seed=632545,save=FALSE,save.graphs=FALSE, 
displayProgress=FALSE)

# Not run (strict time constraints for CRAN)

saemix.fit<-saemix(saemix.model,saemix.data,saemix.options)
# Only 10 bootstrap samples here for speed, please increase to at least 200
theo.case <- saemix.bootstrap(saemix.fit, method="case", nboot=10)
theo.cond <- saemix.bootstrap(saemix.fit, nboot=10)


# Bootstrap for the toenail data
data(toenail.saemix)
saemix.data<-saemixData(name.data=toenail.saemix,name.group=c("id"), name.predictors=c("time","y"), 
 name.response="y", name.covariates=c("treatment"),name.X=c("time"))
 
binary.model<-function(psi,id,xidep) {
  tim<-xidep[,1]
  y<-xidep[,2]
  inter<-psi[id,1]
  slope<-psi[id,2]
  logit<-inter+slope*tim
  pevent<-exp(logit)/(1+exp(logit))
  pobs = (y==0)*(1-pevent)+(y==1)*pevent
  logpdf <- log(pobs)
  return(logpdf)
}
simulBinary<-function(psi,id,xidep) {
    tim<-xidep[,1]
    y<-xidep[,2]
    inter<-psi[id,1]
    slope<-psi[id,2]
    logit<-inter+slope*tim
    pevent<-1/(1+exp(-logit))
    ysim<-rbinom(length(tim),size=1, prob=pevent)
    return(ysim)
    }

saemix.model<-saemixModel(model=binary.model,description="Binary model",
     modeltype="likelihood", simulate.function=simulBinary,
     psi0=matrix(c(-5,-.1,0,0),ncol=2,byrow=TRUE,dimnames=list(NULL,c("inter","slope"))),
     transform.par=c(0,0))
     

saemix.options<-list(seed=1234567,save=FALSE,save.graphs=FALSE, displayProgress=FALSE, 
   nb.chains=10, fim=FALSE)
binary.fit<-saemix(saemix.model,saemix.data,saemix.options)
# Only 10 bootstrap samples here for speed, please increase to at least 200
toenail.case <- saemix.bootstrap(binary.fit, method="case", nboot=10)
toenail.cond <- saemix.bootstrap(binary.fit, nboot=10)

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