saemix: Stochastic Approximation Expectation Maximization (SAEM)...

Description Usage Arguments Details Value Author(s) References See Also Examples

View source: R/main.R

Description

SAEM algorithm perform parameter estimation for nonlinear mixed effects models without any approximation of the model (linearization, quadrature approximation, . . . )

Usage

1
saemix(model, data, control = list())

Arguments

model

an object of class SaemixModel, created by a call to the function saemixModel

data

an object of class SaemixData, created by a call to the function saemixData

control

a list of options, see saemixControl

Details

The SAEM algorithm is a stochastic approximation version of the standard EM algorithm proposed by Kuhn and Lavielle (see reference). Details of the algorithm can be found in the pdf file accompanying the package.

Value

An object of class SaemixObject containing the results of the fit of the data by the non-linear mixed effect model. A summary of the results is printed out to the terminal, and, provided the appropriate options have not been changed, numerical and graphical outputs are saved in a directory.

Author(s)

Emmanuelle Comets <emmanuelle.comet[email protected]>, Audrey Lavenu, Marc Lavielle.

References

Kuhn E, Lavielle M. Maximum likelihood estimation in nonlinear mixed effects models. Computational Statistics and Data Analysis 49, 4 (2005), 1020-1038.

Comets E, Lavenu A, Lavielle M. SAEMIX, an R version of the SAEM algorithm. 20th meeting of the Population Approach Group in Europe, Athens, Greece (2011), Abstr 2173.

See Also

SaemixData,SaemixModel, SaemixObject, saemixControl, plot.saemix

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
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")


# Not run (strict time constraints for CRAN)
# saemix.fit<-saemix(saemix.model,saemix.data,list(seed=632545,directory="newtheo",
# save=FALSE,save.graphs=FALSE))

# Prints a summary of the results
# print(saemix.fit)

# Outputs the estimates of individual parameters
# psi(saemix.fit)

# Shows some diagnostic plots to evaluate the fit
# plot(saemix.fit)

belhal/saemix documentation built on March 22, 2018, 2:34 p.m.