Package 'saemix'

Description

Access slots of an SaemixData object using the object["slot"] format

Usage

## S4 method for signature 'SaemixData'
x[i, j, drop]

## S4 method for signature 'SaemixRepData'
x[i, j, drop]

## S4 replacement method for signature 'SaemixRepData'
x[i, j] <- value

## S4 method for signature 'SaemixSimData'
x[i, j, drop]

## S4 replacement method for signature 'SaemixSimData'
x[i, j] <- value

## S4 replacement method for signature 'SaemixRes'
x[i, j] <- value

Arguments

Description

Access slots of an SaemixModel object using the object["slot"] format

Usage

## S4 method for signature 'SaemixModel'
x[i, j, drop]

Arguments

Description

Access slots of an SaemixObject object using the object["slot"] format

Usage

## S4 method for signature 'SaemixObject'
x[i, j, drop]

Arguments

Description

Access slots of a SaemixRes object using the object["slot"] format

Usage

## S4 method for signature 'SaemixRes'
x[i, j, drop]

Arguments

Description

Joint selection of covariates and random effects in a nonlinear mixed effects model by a backward-type algorithm based on two different versions of BIC for covariate selection and random effects selection respectively. Selection is made among the covariates as such specified in the SaemixData object. Only uncorrelated random effects structures are considered.

Usage

backward.procedure(saemixObject, trace = TRUE)

Arguments

Value

An object of the SaemixObject class storing the covariate model and the covariance structure of random effects of the final model.

Author(s)

Maud Delattre

References

M Delattre, M Lavielle, MA Poursat (2014) A note on BIC in mixed effects models. Electronic Journal of Statistics 8(1) p. 456-475 M Delattre, MA Poursat (2017) BIC strategies for model choice in a population approach. (arXiv:1612.02405)

Description

Check initial fixed effects for an SaemixModel object applied to an SaemixData object

Usage

checkInitialFixedEffects(model, data, psi = c(), id = c(), ...)

Arguments

Details

The function uses the model slot of the SaemixModel object to obtain predictions, using the predictors object. The user is responsible for giving all the predictors needed by the model function. if psi is not given, the predictions will be computed for the population parameters (first line of the psi0 slot) of the object.

The predictions correspond to the structure of the model. For models defined as a structural model, individual plots for the subjects in id overlaying the predictions for the parameters psi and the individual data are shown, and the predictions correspond to f(t_ij, psi). For models defined in terms of their likelihood, the predictions returned correspond to the log-likelihood. No individual graphs are currently available for discrete data models.

Warning: this function is currently under development and the output may change in future versions of the package

Value

the predictions corresponding to the values for each observation in the predictors of either the model f or log-likelihood. For Gaussian data models, the function also plots the data overlayed with the model predictions for each subject in id (where id is the index in the N subjects).

Examples

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,modeltype="structural",
  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")

checkInitialFixedEffects(saemix.model, saemix.data, id=c(1:6))
checkInitialFixedEffects(saemix.model, saemix.data, id=c(1:6), psi=c(0.5, 30, 2)) # better fit

Description

Extract coefficients from an saemix fit

Usage

## S3 method for class 'SaemixObject'
coef(object, ...)

Arguments

Value

a list with 3 components:

fixed

fixed effects

population

a list of population parameters with two elements, a matrix containing the untransformed parameters psi and a matrix containing the transformed parameters phi

individual

a list of individual parameters with two elements, a matrix containing the untransformed parameters psi and a matrix containing the transformed parameters phi

Description

A specific penalty is used for BIC (BIC.cov) when the compared models have in common the structural model and the covariance structure for the random effects (see Delattre et al., 2014).

Usage

compare.saemix(..., method = c("is", "lin", "gq"))

Arguments

Details

Note that the comparison between two or more models will only be valid if they are fitted to the same dataset.

Value

A matrix of information criteria is returned, with at least two columns containing respectively AIC and BIC values for each of the compared models. When the models have in common the structural model and the covariance structure for the random effects, the matrix includes an additional column with BIC.cov values that are more appropriate when the comparison only concerns the covariates.

Author(s)

Maud Delattre

References

M Delattre, M Lavielle, MA Poursat (2014) A note on BIC in mixed effects models. Electronic Journal of Statistics 8(1) p. 456-475

Examples

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

# Definition of models to be compared
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)
}
# Model with one covariate
saemix.model1<-saemixModel(model=model1cpt,modeltype="structural",
  description="One-compartment model, clearance dependent on weight",
  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,0,1,0,0,0),ncol=3,byrow=TRUE))
# Model with two covariates
saemix.model2<-saemixModel(model=model1cpt,modeltype="structural",
   description="One-compartment model, clearance dependent on weight, volume dependent on sex",
   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,0,1,0,1,0),ncol=3,byrow=TRUE))
# Model with three covariates
saemix.model3<-saemixModel(model=model1cpt,modeltype="structural",
  description="One-cpt model, clearance, absorption dependent on weight, volume dependent on sex",
  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(1,0,1,0,1,0),ncol=3,byrow=TRUE))


# Running the main algorithm to estimate the population parameters
saemix.options<-list(seed=632545,save=FALSE,save.graphs=FALSE)
saemix.fit1<-saemix(saemix.model1,saemix.data,saemix.options)
saemix.fit2<-saemix(saemix.model2,saemix.data,saemix.options)
saemix.fit3<-saemix(saemix.model3,saemix.data,saemix.options)

compare.saemix(saemix.fit1, saemix.fit2, saemix.fit3)
compare.saemix(saemix.fit1, saemix.fit2, saemix.fit3, method = "lin")

# We need to explicitly run Gaussian Quadrature if we want to use it in
# the comparisons
saemix.fit1 <- llgq.saemix(saemix.fit1)
saemix.fit2 <- llgq.saemix(saemix.fit2)
saemix.fit3 <- llgq.saemix(saemix.fit3)
compare.saemix(saemix.fit1, saemix.fit2, saemix.fit3, method = "gq")

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

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 [email protected]

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

Description

The cow.saemix data contains records of the weight of 560 cows on 9 or 10 occasions.

Usage

cow.saemix

Format

This data frame contains the following columns:

cow

the unique identifier for each cow

time

time (days)

weight

a numeric vector giving the weight of the cow (kg)

birthyear

year of birth (between 1988 and 1998)

twin

existence of a twin (no=1, yes=2)

birthrank

the rank of birth (beetween 3 and 7)

Details

An exponential model was assumed to describe the weight gain with time: y_ij = A_i (1- B_i exp( - K_i t_ij)) +epsilon_ij

References

JC Pinheiro, DM Bates (2000) Mixed-effects Models in S and S-PLUS, Springer, New York (Appendix A.19)

Examples

data(cow.saemix)
saemix.data<-saemixData(name.data=cow.saemix,header=TRUE,name.group=c("cow"), 
      name.predictors=c("time"),name.response=c("weight"), 
      name.covariates=c("birthyear","twin","birthrank"), 
      units=list(x="days",y="kg",covariates=c("yr","-","-")))

growthcow<-function(psi,id,xidep) {
  x<-xidep[,1]
  a<-psi[id,1]
  b<-psi[id,2]
  k<-psi[id,3]
  f<-a*(1-b*exp(-k*x))
  return(f)
}
saemix.model<-saemixModel(model=growthcow,
      description="Exponential growth model", 
      psi0=matrix(c(700,0.9,0.02,0,0,0),ncol=3,byrow=TRUE, 
        dimnames=list(NULL,c("A","B","k"))),transform.par=c(1,1,1),fixed.estim=c(1,1,1), 
      covariate.model=matrix(c(0,0,0),ncol=3,byrow=TRUE), 
      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(algorithms=c(1,1,1),nb.chains=1,nbiter.saemix=c(200,100), 
             seed=4526,save=FALSE,save.graphs=FALSE,displayProgress=FALSE)

# Plotting the data
plot(saemix.data,xlab="Time (day)",ylab="Weight of the cow (kg)")

saemix.fit<-saemix(saemix.model,saemix.data,saemix.options)

Description

Create saemix objects either with empty results or with parameters set by the user. This is an internal function used as a preliminary step to obtain predictions for new data and is not intended to be used directly.

Usage

createSaemixObject.empty(model, data, control = list())

Arguments

Details

with createSaemixObject.empty, the data component is set to the data object, the model component is set to the model object, and the result component is empty

with createSaemixObject.initial, the data and model are set as with createSaemixObject.empty, but the population parameter estimates are used to initialise the result component as in the initialisation step of the algorithm (initialiseMainAlgo)

Value

an object of class "SaemixObject".

Examples

# TODO

Description

These functions create bootstrapped datasets using the bootstrap methods described in saemix.bootstrap

Usage

dataGen.case(saemixObject)

dataGen.NP(saemixObject, nsamp = 0, eta.sampc = NULL, conditional = FALSE)

dataGen.Par(saemixObject)

Arguments

Description

These functions produce default sets of plots, corresponding to diagnostic or individual fits.

Usage

default.saemix.plots(saemixObject, ...)

Arguments

Details

These functions are wrapper functions designed to produce default sets of plots to help the user assess their model fits.

Value

Depending on the type argument, the following plots are produced:

• default.saemix.plots by default, the following plots are produced: a plot of the data, convergence plots, plot of the likelihood by importance sampling (if computed), plots of observations versus predictions, scatterplots and distribution of residuals, boxplot of the random effects, correlations between random effects, distribution of the parameters, VPC

• basic.gof basic goodness-of-fit plots: convergence plots, plot of the likelihood by importance sampling (if computed), plots of observations versus predictions

• advanced.gof advanced goodness-of-fit plots: scatterplots and distribution of residuals, VPC,...

• covariate.fits plots of all estimated parameters versus all covariates in the dataset

• individual.fits plots of individual predictions (line) overlayed on individual observations (dots) for all subjects in the dataset

Author(s)

Emmanuelle Comets [email protected], 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

saemix, saemix.plot.data, saemix.plot.setoptions, plot.saemix

Examples

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

# Reducing the number of iterations due to time constraints for CRAN
saemix.options<-list(seed=632545,save=FALSE,save.graphs=FALSE,nbiter.saemix=c(100,100))

saemix.fit<-saemix(saemix.model,saemix.data,saemix.options)

default.saemix.plots(saemix.fit)

# Not run (time constraints for CRAN)
# basic.gof(saemix.fit)

# Not run (time constraints for CRAN)
# advanced.gof(saemix.fit)

individual.fits(saemix.fit)

Description

This function provides VPC plots for non Gaussian data models (work in progress)

Usage

discreteVPC(object, outcome = "categorical", verbose = FALSE, ...)

Arguments

Details

This function is a very rough first attempt at automatically creating VPC plots for models defined through their log-likelihood (categorical, count, or TTE data). It makes use of the new element simulate.function in the model component of the object

• for TTE data, a KM-VPC plot will be produced

• for count, categorical and binary data, a plot showing the proportion of each score/category across time will be shown along with the corresponding prediction intervals from the model These plots can be stratified over a covariate in the data set (currently only categorical covariates) by passing an argument which.cov='name' to the call

Author(s)

Emmanuelle Comets [email protected]

References

Brendel, K, Comets, E, Laffont, C, Laveille, C, Mentre, F. Metrics for external model evaluation with an application to the population pharmacokinetics of gliclazide, Pharmaceutical Research 23 (2006), 2036-2049.

Holford, N. The Visual Predictive Check: superiority to standard diagnostic (Rorschach) plots (Abstract 738), in: 14th Meeting of the Population Approach Group in Europe, Pamplona, Spain, 2005.

Ron Keizer, tutorials on VPC TODO

See Also

SaemixObject, saemix, saemix.plot.vpc, simulateDiscreteSaemix

Examples

data(lung.saemix)

saemix.data<-saemixData(name.data=lung.saemix,header=TRUE,name.group=c("id"),
name.predictors=c("time","status","cens"),name.response=c("status"),
name.covariates=c("age", "sex", "ph.ecog", "ph.karno", "pat.karno", "wt.loss","meal.cal"),
units=list(x="days",y="",covariates=c("yr","","-","%","%","cal","pounds")))

weibulltte.model<-function(psi,id,xidep) {
  T<-xidep[,1]
  y<-xidep[,2] # events (1=event, 0=no event)
  cens<-which(xidep[,3]==1) # censoring times (subject specific)
  init <- which(T==0)
  lambda <- psi[id,1] # Parameters of the Weibull model
  beta <- psi[id,2]
  Nj <- length(T)
  ind <- setdiff(1:Nj, append(init,cens)) # indices of events
  hazard <- (beta/lambda)*(T/lambda)^(beta-1) # H'
  H <- (T/lambda)^beta # H
  logpdf <- rep(0,Nj) # ln(l(T=0))=0
  logpdf[cens] <- -H[cens] + H[cens-1] # ln(l(T=censoring time))
  logpdf[ind] <- -H[ind] + H[ind-1] + log(hazard[ind]) # ln(l(T=event time))
  return(logpdf)
}

simulateWeibullTTE <- function(psi,id,xidep) {
  T<-xidep[,1]
  y<-xidep[,2] # events (1=event, 0=no event)
  cens<-which(xidep[,3]==1) # censoring times (subject specific)
  init <- which(T==0)
  lambda <- psi[,1] # Parameters of the Weibull model
  beta <- psi[,2]
  tevent<-T
  Vj<-runif(dim(psi)[1])
  tsim<-lambda*(-log(Vj))^(1/beta) # nsuj events
  tevent[T>0]<-tsim
  tevent[tevent[cens]>T[cens]] <- T[tevent[cens]>T[cens]]
  return(tevent)
  }
saemix.model<-saemixModel(model=weibulltte.model,description="time model",modeltype="likelihood",
       simulate.function = simulateWeibullTTE,
                psi0=matrix(c(1,2),ncol=2,byrow=TRUE,dimnames=list(NULL,  c("lambda","beta"))),
                transform.par=c(1,1),covariance.model=matrix(c(1,0,0,0),ncol=2, byrow=TRUE))
saemix.options<-list(seed=632545,save=FALSE,save.graphs=FALSE, displayProgress=FALSE)


tte.fit<-saemix(saemix.model,saemix.data,saemix.options)
simtte.fit <- simulateDiscreteSaemix(tte.fit, nsim=100)
gpl <- discreteVPC(simtte.fit, outcome="TTE")
plot(gpl)

Description

This function provides VPC plots for time-to-event data models (work in progress)

Usage

discreteVPCTTE(
  object,
  ngrid = 200,
  interpolation.method = "step",
  verbose = FALSE,
  ...
)

Arguments

Details

add details on TTE VPC, RTTE VPC, etc...

Author(s)

Emmanuelle Comets [email protected]

See Also

SaemixObject, saemix, saemix.plot.vpc, simulateDiscreteSaemix

Tutorials on TTE-VPC TODO

Description

The epilepsy data from Thall and Vail (1990), available from the MASS package, records two-week seizure counts for 59 epileptics. The number of seizures was recorded for a baseline period of 8 weeks, and then patients were randomly assigned to a treatment group or a control group. Counts were then recorded for four successive two-week periods. The subject's age is the only covariate. See the documentation for epil in the MASS package for details on the dataset.

Source

MASS package in R

References

P Thall, S Vail (1990). Some covariance models for longitudinal count data with overdispersion. Biometrics 46(3):657-71.

Examples

# You need to have MASS installed to successfully run this example
if (requireNamespace("MASS")) {
  
  epilepsy<-MASS::epil
  saemix.data<-saemixData(name.data=epilepsy, name.group=c("subject"),
     name.predictors=c("period","y"),name.response=c("y"),
     name.covariates=c("trt","base", "age"), units=list(x="2-week",y="",covariates=c("","","yr")))
  ## Poisson model with one parameter
  countPoi<-function(psi,id,xidep) { 
    y<-xidep[,2]
    lambda<-psi[id,1]
    logp <- -lambda + y*log(lambda) - log(factorial(y))
    return(logp)
    }
 saemix.model<-saemixModel(model=countPoi,description="Count model Poisson",modeltype="likelihood", 
   psi0=matrix(c(0.5),ncol=1,byrow=TRUE,dimnames=list(NULL, c("lambda"))), transform.par=c(1))
 
  saemix.options<-list(seed=632545,save=FALSE,save.graphs=FALSE, displayProgress=FALSE)
  poisson.fit<-saemix(saemix.model,saemix.data,saemix.options)
  
}

Description

Estimate by linearisation the Fisher Information Matrix and the standard error of the estimated parameters.

Usage

fim.saemix(saemixObject)

Arguments

Details

The inverse of the Fisher Information Matrix provides an estimate of the variance of the estimated parameters theta. This matrix cannot be computed in closed-form for nonlinear mixed-effect models; instead, an approximation is obtained as the Fisher Information Matrix of the Gaussian model deduced from the nonlinear mixed effects model after linearisation of the function f around the conditional expectation of the individual Gaussian parameters. This matrix is a block matrix (no correlations between the estimated fixed effects and the estimated variances).

Value

The function returns an updated version of the object saemix.fit in which the following elements have been added:

se.fixed:

standard error of fixed effects, obtained as part of the diagonal of the inverse of the Fisher Information Matrix (only when fim.saemix has been run, or when saemix.options$algorithms[2] is 1)

se.omega:

standard error of the variance of random effects, obtained as part of the diagonal of the inverse of the Fisher Information Matrix (only when fim.saemix has been run, or when the saemix.options$algorithms[2] is 1)

se.res:

standard error of the parameters of the residual error model, obtained as part of the diagonal of the inverse of the Fisher Information Matrix (only when fim.saemix has been run, or when the saemix.options$algorithms[2] is 1)

fim:

Fisher Information Matrix

ll.lin:

likelihood calculated by linearisation

Author(s)

Emmanuelle Comets [email protected], 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

SaemixObject,saemix

Examples

# Running the main algorithm to estimate the population parameters
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)

# Not run (strict time constraints for CRAN)
# saemix.fit<-saemix(saemix.model,saemix.data,saemix.options)

# Estimating the Fisher Information Matrix using the result of saemix 
# & returning the result in the same object
# fim.saemix(saemix.fit)

Description

fitted is a generic function which extracts model predictions from objects returned by modelling functions

Usage

## S3 method for class 'SaemixRes'
fitted(object, type = c("ipred", "ypred", "ppred", "icpred"), ...)

## S3 method for class 'SaemixObject'
fitted(object, type = c("ipred", "ypred", "ppred", "icpred"), ...)

Arguments

Value

Model predictions

Description

Joint selection of covariates and random effects in a nonlinear mixed effects model by a forward-type algorithm based on two different versions of BIC for covariate selection and random effects selection respectively. Selection is made among the covariates as such specified in the SaemixData object. Only uncorrelated random effects structures are considered.

Usage

forward.procedure(saemixObject, trace = TRUE)

Arguments

Value

An object of the SaemixObject class storing the covariate model and the covariance structure of random effects of the final model.

Author(s)

Maud Delattre

References

M Delattre, M Lavielle, MA Poursat (2014) A note on BIC in mixed effects models. Electronic Journal of Statistics 8(1) p. 456-475 M Delattre, MA Poursat (2017) BIC strategies for model choice in a population approach. (arXiv:1612.02405)

Description

Constructor functions for Classes in the saemix package

Usage

## S4 method for signature 'SaemixData'
initialize(
  .Object,
  name.data,
  header,
  sep,
  na,
  name.group,
  name.predictors,
  name.response,
  name.covariates,
  name.X,
  units,
  name.mdv,
  name.cens,
  name.occ,
  name.ytype,
  verbose = TRUE,
  automatic = TRUE
)

## S4 method for signature 'SaemixRepData'
initialize(.Object, data = NULL, nb.chains = 1)

## S4 method for signature 'SaemixSimData'
initialize(.Object, data = NULL, datasim = NULL)

## S4 method for signature 'SaemixModel'
initialize(
  .Object,
  model,
  description,
  modeltype,
  psi0,
  name.response,
  name.sigma,
  transform.par,
  fixed.estim,
  error.model,
  covariate.model,
  covariance.model,
  omega.init,
  error.init,
  name.modpar,
  verbose = TRUE
)

## S4 method for signature 'SaemixRes'
initialize(
  .Object,
  status = "empty",
  modeltype,
  name.fixed,
  name.random,
  name.sigma,
  fixed.effects,
  fixed.psi,
  betaC,
  betas,
  omega,
  respar,
  cond.mean.phi,
  cond.var.phi,
  mean.phi,
  phi,
  phi.samp,
  parpop,
  allpar,
  MCOV
)

## S4 method for signature 'SaemixObject'
initialize(.Object, data, model, options = list())

Arguments

Methods

list("signature(.Object = \"SaemixData\")")

create a SaemixData object. Please use the saemixData function.

list("signature(.Object = \"SaemixModel\")")

create a SaemixModel object Please use the saemixModel function.

list("signature(.Object = \"SaemixObject\")")

create a SaemixObject object. This object is obtained after a successful call to saemix

list("signature(.Object = \"SaemixRepData\")")

create a SaemixRepData object

list("signature(.Object = \"SaemixRes\")")

create a SaemixRes object

list("signature(.Object = \"SaemixSimData\")")

create a SaemixSimData object

Description

The knee.saemix data represents pain scores recorded in a clinical study in 127 patients with sport related injuries treated with two different therapies. After 3,7 and 10 days of treatment the pain occuring during knee movement was observed.

Usage

knee.saemix

Format

This data frame contains the following columns:

id

subject index in file

time

time of measurement (in days)

y

knee pain (0=none to 4=severe)

Age

patient age (scaled and centered)

Sex

patient gender (0=male, 1=female)

RD

moderate knee pain (defined as pain score 2 or more)

treatment

treatment indicator (0=placebo, 1=treatment)

Age2

patient age, squared (Age^2)

#'

Details

The data in the knee.saemix was reformatted from the knee dataset provided by the catdata package (see data(knee, package="catdata")). A time column was added representing the day of the measurement (with 0 being the baseline value) and each observation corresponds to a different line in the dataset. Treatment was recoded as 0/1 (placebo/treatment), gender as 0/1 (male/female) and Age2 represents the squared of centered Age.

Please refer to the PDF documentation (chapter 4, section 4.6) for more details on the analysis, including examples of diagnostic plots.

Source

catdata package in R

References

G Tutz (2012), Regression for Categorical Data, Cambridge University Press.

#' @examples data(knee.saemix)

#' @keywords datasets

Description

Estimate the log-likelihood using Gaussian Quadrature (multidimensional grid)

Usage

llgq.saemix(saemixObject)

Arguments

Details

The likelihood of the observations is estimated using Gaussian Quadrature (see documentation).

Value

the log-likelihood estimated by Gaussian Quadrature

Author(s)

Emmanuelle Comets [email protected], 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

SaemixObject,saemix,llis.saemix

Examples

# Running the main algorithm to estimate the population parameters
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)

# Not run (strict time constraints for CRAN)
# saemix.fit<-saemix(saemix.model,saemix.data,saemix.options)

# Estimating the likelihood by Gaussian Quadrature using the result of saemix 
# & returning the result in the same object
# saemix.fit<-llgq.saemix(saemix.fit)

Description

Estimate the log-likelihood using Importance Sampling

Usage

llis.saemix(saemixObject)

Arguments

Details

The likelihood of the observations is estimated without any approximation using a Monte-Carlo approach (see documentation).

Value

the log-likelihood estimated by Importance Sampling

Author(s)

Emmanuelle Comets [email protected], 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

SaemixObject,saemix,llgq.saemix

Examples

# Running the main algorithm to estimate the population parameters
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", modeltype="structural",
  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)

# Not run (strict time constraints for CRAN)
# saemix.fit<-saemix(saemix.model,saemix.data,saemix.options)

# Estimating the likelihood by importance sampling using the result of saemix 
# & returning the result in the same object
# saemix.fit<-llis.saemix(saemix.fit)

Description

The likelihood in saemix can be computed by one of three methods: linearisation (linearisation of the model), importance sampling (stochastic integration) and gaussian quadrature (numerical integration). The linearised likelihood is obtained as a byproduct of the computation of the Fisher Information Matrix (argument FIM=TRUE in the options given to the saemix function). If no method argument is given, this function will attempt to extract the likelihood computed by importance sampling (method="is"), unless the object contains the likelihood computed by linearisation, in which case the function will extract this component instead. If the requested likelihood is not present in the object, it will be computed and aded to the object before returning.

Usage

## S3 method for class 'SaemixObject'
logLik(object, method = c("is", "lin", "gq"), ...)

## S3 method for class 'SaemixObject'
AIC(object, method = c("is", "lin", "gq"), ..., k = 2)

## S3 method for class 'SaemixObject'
BIC(object, method = c("is", "lin", "gq"), ...)

## S3 method for class 'covariate'
BIC(object, method = c("is", "lin", "gq"), ...)

Arguments

Details

BIC.covariate implements the computation of the BIC from Delattre et al. 2014.

Value

Returns the selected statistical criterion (log-likelihood, AIC, BIC) extracted from the SaemixObject, computed with the 'method' argument if given (defaults to IS).

Author(s)

Emmanuelle Comets [email protected]

Audrey Lavenu

Marc Lavielle

Maud Delattre

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

AIC,BIC, saemixControl, saemix

Description

The lung.saemix contains survival data in patients with advanced lung cancer from the North Central Cancer Treatment Group. Performance scores rate how well the patient can perform usual daily activities. This data is available in the survival library for R and has been reformatted here for use in saemix (see details).

Usage

lung.saemix

Format

This data frame contains the following columns:

id

subject index in file

inst

institution code

time

observation time since the beginning of follow-up

status

0=alive, 1=dead

cens

0=observed, 1=censored

sex

patient gender (0=male, 1=female)

age

age in years

ph.ecog

ECOG performance score as rated by the physician. 0=asymptomatic, 1= symptomatic but completely ambulatory, 2= in bed <50% of the day, 3= in bed > 50% of the day but not bedbound, 4 = bedbound

ph.karno

Karnofsky performance score (bad=0-good=100) rated by physician (%)

pat.karno

Karnofsky performance score (bad=0-good=100) rated by patient (%)

meal.cal

calories consumed at meals (cal)

wt.loss

weight loss in last six months (pounds)

Details

The data in the lung.saemix was reformatted from the lung cancer dataset (see data(cancer, package="survival")). Patients with missing age, sex, institution or physician assessments were removed from the dataset. Status was recoded as 1 for death and 0 for a censored event, and a censoring column was added to denote whether the patient was dead or alive at the time of the last observation. For saemix, a line at time=0 was added for all subjects. Finally, subjects were numbered consecutively from 0 to 1.

Source

Terry Therneau from the survival package in R

References

CL Loprinzi, JA Laurie, HS Wieand, JE Krook, PJ Novotny, JW Kugler, et al. (1994). Prospective evaluation of prognostic variables from patient-completed questionnaires. North Central Cancer Treatment Group. Journal of Clinical Oncology. 12(3):601-7.

Examples

data(lung.saemix)

saemix.data<-saemixData(name.data=lung.saemix,header=TRUE,name.group=c("id"),
name.predictors=c("time","status","cens"),name.response=c("status"),
name.covariates=c("age", "sex", "ph.ecog", "ph.karno", "pat.karno", "wt.loss","meal.cal"),
units=list(x="days",y="",covariates=c("yr","","-","%","%","cal","pounds")))
weibulltte.model<-function(psi,id,xidep) {
  T<-xidep[,1]
  y<-xidep[,2] # events (1=event, 0=no event)
  cens<-which(xidep[,3]==1) # censoring times (subject specific)
  init <- which(T==0)
  lambda <- psi[id,1] # Parameters of the Weibull model
  beta <- psi[id,2]
  Nj <- length(T)
  ind <- setdiff(1:Nj, append(init,cens)) # indices of events
  hazard <- (beta/lambda)*(T/lambda)^(beta-1) # H'
  H <- (T/lambda)^beta # H
  logpdf <- rep(0,Nj) # ln(l(T=0))=0
  logpdf[cens] <- -H[cens] + H[cens-1] # ln(l(T=censoring time))
  logpdf[ind] <- -H[ind] + H[ind-1] + log(hazard[ind]) # ln(l(T=event time))
  return(logpdf)
}
saemix.model<-saemixModel(model=weibulltte.model,description="time model",modeltype="likelihood",
                psi0=matrix(c(1,2),ncol=2,byrow=TRUE,dimnames=list(NULL,  c("lambda","beta"))),
                transform.par=c(1,1),covariance.model=matrix(c(1,0,0,0),ncol=2, byrow=TRUE))
saemix.options<-list(seed=632545,save=FALSE,save.graphs=FALSE, displayProgress=FALSE)

tte.fit<-saemix(saemix.model,saemix.data,saemix.options)

# The fit from saemix using the above Weibull model may be compared 
# to the non-parametric KM estimate
## Not run: 
library(survival)
  lung.surv<-lung.saemix[lung.saemix$time>0,]
  lung.surv$status<-lung.surv$status+1
  Surv(lung.surv$time, lung.surv$status) # 1=censored, 2=dead
  f1 <- survfit(Surv(time, status) ~ 1, data = lung.surv)
  xtim<-seq(0,max(lung.saemix$time), length.out=200)
  estpar<-tte.fit@results@fixed.effects
  ypred<-exp(-(xtim/estpar[1])^(estpar[2]))
  plot(f1, xlab = "Days", ylab = "Overall survival probability")
  lines(xtim,ypred, col="red",lwd=2)

## End(Not run)

Description

Compute the estimates of the individual parameters PSI_i (conditional mode - Maximum A Posteriori)

Usage

map.saemix(saemixObject)

Arguments

Details

The MCMC procedure is used to estimate the conditional mode (or Maximum A Posteriori) m(phi_i |yi ; hattheta) = Argmax_phi_i p(phi_i |yi ; hattheta)

Value

Author(s)

Emmanuelle Comets [email protected], 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

SaemixObject,saemix

Examples

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)

# Not run (strict time constraints for CRAN)
# saemix.fit<-saemix(saemix.model,saemix.data,saemix.options)

# Estimating the individual parameters using the result of saemix 
# & returning the result in the same object
# saemix.fit<-map.saemix(saemix.fit)

Description

Extract or replace the diagonal of a matrix, or construct a diagonal matrix (replace diag function from R-base)

Usage

mydiag(x = 1, nrow, ncol)

Arguments

Value

If x is a matrix then diag(x) returns the diagonal of x. The resulting vector will have names if the matrix x has matching column and rownames.

Author(s)

Emmanuelle Comets [email protected], Audrey Lavenu, Marc Lavielle.

See Also

diag

Examples

mydiag(1)
mydiag(c(1,2))

Description

This function uses the npde library to compute normalised prediction distribution errors (npde) and normalised prediction discrepancies (npd). The simulations can also be used to plot VPCs. As of version 3.0, the plot functions for diagnostics involving VPC, npde or npd are deprecated and the user will be redirected to the current proposed method (creating an NpdeObject and using the plot functions from the library npde).

Usage

npdeSaemix(saemixObject, nsim = 1000)

Arguments

Details

Since version 3.0, the saemix package depends on the npde package, which computes the npd/npde and produces graphs. See the documentation for npde for details on the computation methods See the PDF documentation and the bookdown https://iame-researchcenter.github.io/npde_bookdown/ for details on the different plots available.

Value

An object of class NpdeObject

Author(s)

Emmanuelle Comets [email protected]

References

E Comets, K Brendel, F Mentre (2008). Computing normalised prediction distribution errors to evaluate nonlinear mixed-effect models: the npde add-on package for R. Computer Methods and Programs in Biomedicine, 90:154-66.

K Brendel, E Comets, C Laffont, C Laveille, F Mentre (2006). Metrics for external model evaluation with an application to the population pharmacokinetics of gliclazide. Pharmaceutical Research, 23:2036–49

PDF documentation for npde 3.0: https://github.com/ecomets/npde30/blob/main/userguide_npde_3.1.pdf

See Also

npde.graphs, gof.test

NpdeObject npde.plot.select autonpde

npde.plot.scatterplot npde.plot.dist

Examples

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
# Works interactively but not in the contained environment of CRAN (it looks for a datafile 
# instesad of finding the dataset in the environment)
# saemix.fit<-saemix(saemix.model,saemix.data,saemix.options)
# npde.obj<-npdeSaemix(saemix.fit)
# plot(npde.obj)
# plot(npde.obj, plot.type="vpc")
# plot(npde.obj, plot.type="covariates")
# plot(npde.obj, plot.type="cov.x.scatter")
# plot(npde.obj, plot.type="cov.ecdf")

Description

The oxboys.saemix data collects the height and age of students in Oxford measured over time. The data frame has 234 rows and 4 columns.

Usage

oxboys.saemix

Format

This data frame contains the following columns:

Subject

an ordered factor giving a unique identifier for each boy in the experiment

age

a numeric vector giving the standardized age (dimensionless)

height

a numeric vector giving the height of the boy (cm)

Occasion

an ordered factor - the result of converting 'age' from a continuous variable to a count so these slightly unbalanced data can be analyzed as balanced

Details

These data are described in Goldstein (1987) as data on the height of a selection of boys from Oxford, England versus a standardized age. The dataset can be found in the package nlme. We use an linear model for this data: y_ij = Base_i + slope_i x_ij +epsilon_ij

References

JC Pinheiro, DM Bates (2000) Mixed-effects Models in S and S-PLUS, Springer, New York (Appendix A.19)

Examples

data(oxboys.saemix)
saemix.data<-saemixData(name.data=oxboys.saemix,header=TRUE,
      name.group=c("Subject"),name.predictors=c("age"),name.response=c("height"),
      units=list(x="yr",y="cm"))

# plot the data
plot(saemix.data)

growth.linear<-function(psi,id,xidep) {
  x<-xidep[,1]
  base<-psi[id,1]
  slope<-psi[id,2]
  f<-base+slope*x
  return(f)
}
saemix.model<-saemixModel(model=growth.linear,description="Linear model",
      psi0=matrix(c(140,1),ncol=2,byrow=TRUE,dimnames=list(NULL,c("base","slope"))),
      transform.par=c(1,0),covariance.model=matrix(c(1,1,1,1),ncol=2,byrow=TRUE), 
      error.model="constant")

saemix.options<-list(algorithms=c(1,1,1),nb.chains=1,seed=201004,
      save=FALSE,save.graphs=FALSE,displayProgress=FALSE)

saemix.fit<-saemix(saemix.model,saemix.data,saemix.options)

Description

The PD1.saemix and PD2.saemix data frames were simulated according to an Emax dose-response model.

Usage

PD1.saemix

PD2.saemix

Format

This data frame contains the following columns:

subject

an variable identifying the subject on whom the observation was made. The ordering is by the dose at which the observation was made.

dose

simulated dose.

response

simulated response

gender

gender (0 for male, 1 for female)

Details

These examples were used by P. Girard and F. Mentre for the symposium dedicated to Comparison of Algorithms Using Simulated Data Sets and Blind Analysis, that took place in Lyon, France, September 2004. The datasets contain 100 individuals, each receiving 3 different doses:(0, 10, 90), (5, 25, 65) or (0, 20, 30). It was assumed that doses were given in a cross-over study with sufficient wash out period to avoid carry over. Responses (y_ij) were simulated with the following pharmacodynamic model: y_ij = E0_i + D_ij Emax_i/(D_ij + ED50_i) +epsilon_ij The individual parameters were simulated according to log (E0_i) = log (E0) + eta_i1 log (Emax_i) = log (Emax) + eta_i2 log (E50_i) = log (E50) + beta w_i + eta_i3

PD1.saemix contains the data simulated with a gender effect, beta=0.3. PD2.saemix contains the data simulated without a gender effect, beta=0.

References

P Girard, F Mentre (2004). Comparison of Algorithms Using Simulated Data Sets and Blind Analysis workshop, Lyon, France.

Examples

data(PD1.saemix)
saemix.data<-saemixData(name.data=PD1.saemix,header=TRUE,name.group=c("subject"),
      name.predictors=c("dose"),name.response=c("response"),
      name.covariates=c("gender"), units=list(x="mg",y="-",covariates=c("-")))

modelemax<-function(psi,id,xidep) {
# input:
#   psi : matrix of parameters (3 columns, E0, Emax, EC50)
#   id : vector of indices 
#   xidep : dependent variables (same nb of rows as length of id)
# returns:
#   a vector of predictions of length equal to length of id
  dose<-xidep[,1]
  e0<-psi[id,1]
  emax<-psi[id,2]
  e50<-psi[id,3]
  f<-e0+emax*dose/(e50+dose)
  return(f)
}

# Plotting the data
plot(saemix.data,main="Simulated data PD1")

# Compare models with and without covariates with LL by Importance Sampling
model1<-saemixModel(model=modelemax,description="Emax growth model", 
       psi0=matrix(c(20,300,20,0,0,0),ncol=3,byrow=TRUE,dimnames=list(NULL,
       c("E0","Emax","EC50"))), transform.par=c(1,1,1),
       covariate.model=matrix(c(0,0,0), ncol=3,byrow=TRUE),fixed.estim=c(1,1,1))

model2<-saemixModel(model=modelemax,description="Emax growth model", 
       psi0=matrix(c(20,300,20,0,0,0),ncol=3,byrow=TRUE,dimnames=list(NULL, 
       c("E0","Emax","EC50"))), transform.par=c(1,1,1),
       covariate.model=matrix(c(0,0,1), ncol=3,byrow=TRUE),fixed.estim=c(1,1,1))

# SE not computed as not needed for the test
saemix.options<-list(algorithms=c(0,1,1),nb.chains=3,seed=765754, 
       nbiter.saemix=c(500,300),save=FALSE,save.graphs=FALSE,displayProgress=FALSE)

fit1<-saemix(model1,saemix.data,saemix.options)
fit2<-saemix(model2,saemix.data,saemix.options)
wstat<-(-2)*(fit1["results"]["ll.is"]-fit2["results"]["ll.is"])

cat("LRT test for covariate effect on EC50: p-value=",1-pchisq(wstat,1),"\n")

Description

Methods for function plot

Methods

list("signature(x = \"ANY\")")

default plot function ?

list("signature(x = \"SaemixData\")")

Plots the data. Defaults to a spaghetti plot of response versus predictor, with lines joining the data for one individual.

list("signature(x = \"SaemixModel\")")

Plots prediction of the model

list("signature(x = \"SaemixObject\")")

This method gives access to a number of plots that can be performed on a SaemixObject

list("signature(x = \"SaemixSimData\")")

Plots simulated datasets

Description

This function will plot predictions obtained from an SaemixModel object over a given range of X. Additional predictors may be passed on to the function using the predictors argument.

Usage

## S4 method for signature 'SaemixModel,ANY'
plot(x, y, range = c(0, 1), psi = NULL, predictors = NULL, ...)

Arguments

Examples

# Note that we have written the PK model so that time is the first predictor (xidep[,1]) 
# and dose the second
model1cpt<-function(psi,id,xidep) { 
     tim<-xidep[,1]  
     dose<-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)
     }
 x<-saemixModel(model=model1cpt,description="One-compartment model with first-order absorption", 
               psi0=matrix(c(1.5,30,1), ncol=3,byrow=TRUE, dimnames=list(NULL, c("ka","V","CL"))))
 # Plot the model over 0-24h, using the parameters given in psi0 and a dose of 300
 plot(x, range=c(0,24), predictors=300, verbose=TRUE)
 # Plot the model over 0-24h, using another set of parameters and a dose of 350
 plot(x, range=c(0,24), psi=c(1.5,20,2), predictors=350, verbose=TRUE)

Description

Plot model predictions for a new dataset. If the dataset is large, only the first 20 subjects (id's) will be shown.

Usage

## S4 method for signature 'SaemixModel,SaemixData'
plot(x, y, ...)

Arguments

Details

The function uses the model slot of the SaemixModel object to obtain predictions, using the dataset contained in the SaemixData object. The user is responsible for making sure data and model match. If psi is not given, the predictions will be computed for the population parameters (first line of the psi0 slot) of the object. If psi is given, the number of columns in psi (or the number of elements of psi, if psi is given as a vector) should match the number of parameters in the model, otherwise an error message will be shown and the function will return empty. If psi is a dataframe, each line will be used for a separate subject of the smx.data object. Elements of psi will be recycled if psi has less lines than the number of subjects in the dataset.

Currently this function only works for models defined as 'structural'.

Value

a ggplot object

Examples

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,modeltype="structural",
  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")

plot(saemix.model, saemix.data)
plot(saemix.model, saemix.data, psi=c(2, 40, 3))
indpsi<-data.frame(ka=2, V=seq(25,47,2), CL=seq(2.5,4.7, 0.2))
plot(saemix.model, saemix.data, psi=indpsi)

Description

Several plots (selectable by the type argument) are currently available: convergence plot, individual plots, predictions versus observations, distribution plots, VPC, residual plots, mirror.

Usage

## S4 method for signature 'SaemixObject,ANY'
plot(x, y, ...)

Arguments

Details

This is the generic plot function for an SaemixObject object, which implements different graphs related to the algorithm (convergence plots, likelihood estimation) as well as diagnostic graphs. A description is provided in the PDF documentation. Arguments such as main, xlab, etc... that can be given to the generic plot function may be used, and will be interpreted according to the type of plot that is to be drawn.

A special argument plot.type can be set to determine the type of plot; it can be one of:

data:

A spaghetti plot of the data,displaying the observed data y as a function of the regression variable (time for a PK application)

convergence:

For each parameter in the model, this plot shows the evolution of the parameter estimate versus the iteration number

likelihood:

Graph showing the evolution of the log-likelihood during the estimation by importance sampling

observations.vs.predictions:

Plot of the predictions computed with the population parameters versus the observations (left), and plot of the predictions computed with the individual parameters versus the observations (right)

residuals.scatter:

Scatterplot of the residuals versus the predictor (top) and versus predictions (bottom), for weighted residuals (population residuals, left), individual weighted residuals (middle) and npde (right).

residuals.distribution:

Distribution of the residuals, plotted as histogram (top) and as a QQ-plot (bottom), for weighted residuals (population residuals, left), individual weighted residuals (middle) and npde (right).

individual.fit:

Individual fits are obtained using the individual parameters with the individual covariates

population.fit:

Population fits are obtained using the population parameters with the individual covariates

both.fit:

Individual fits, superposing fits obtained using the population parameters with the individual covariates (red) and using the individual parameters with the individual covariates (green)

mirror:

Mirror plots assessing the compatibility of simulated data compared to the original

marginal.distribution:

Distribution of the parameters (conditional on covariates when some are included in the model). A histogram of individual parameter estimates can be overlayed on the plot, but it should be noted that the histogram does not make sense when there are covariates influencing the parameters and a warning will be displayed

random.effects:

Boxplot of the random effects

correlations:

Correlation between the random effects

parameters.vs.covariates:

Plots of the estimates of the individual parameters versus the covariates, using scatterplot for continuous covariates, boxplot for categorical covariates

randeff.vs.covariates:

Plots of the estimates of the random effects versus the covariates, using scatterplot for continuous covariates, boxplot for categorical covariates

npde:

Plots 4 graphs to evaluate the shape of the distribution of the normalised prediction distribution errors (npde)

vpc:

Visual Predictive Check, with options to include the prediction intervals around the boundaries of the selected interval as well as around the median (50th percentile of the simulated data).

In addition, the following values for plot.type produce a series of plots:

reduced:

produces the following plots: plot of the data, convergence plots, plot of the likelihood by importance sampling (if computed), plots of observations versus predictions. This is the default behaviour of the plot function applied to an SaemixObject object

full:

produces the following plots: plot of the data, convergence plots, plot of the likelihood by importance sampling (if computed), plots of observations versus predictions, scatterplots and distribution of residuals, VPC, npde, boxplot of the random effects, distribution of the parameters, correlations between random effects, plots of the relationships between individually estimated parameters and covariates, plots of the relationships between individually estimated random effects and covariates

Each plot can be customised by modifying options, either through a list of options set by the saemix.plot.setoptions function, or on the fly by passing an option in the call to the plot (see examples).

Value

None

Author(s)

Emmanuelle Comets [email protected], 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

SaemixObject,saemix, saemix.plot.setoptions, saemix.plot.select, saemix.plot.data

Examples

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)

# Not run (strict time constraints for CRAN)
# saemix.fit<-saemix(saemix.model,saemix.data,saemix.options)

# Set of default plots
# plot(saemix.fit)

# Data
# plot(saemix.fit,plot.type="data")

# Convergence
# plot(saemix.fit,plot.type="convergence")

# Individual plot for subject 1, smoothed
# plot(saemix.fit,plot.type="individual.fit",ilist=1,smooth=TRUE)

# Individual plot for subject 1 to 12, with ask set to TRUE 
# (the system will pause before a new graph is produced)
# plot(saemix.fit,plot.type="individual.fit",ilist=c(1:12),ask=TRUE)

# Diagnostic plot: observations versus population predictions
# par(mfrow=c(1,1))
# plot(saemix.fit,plot.type="observations.vs.predictions",level=0,new=FALSE)

# LL by Importance Sampling
# plot(saemix.fit,plot.type="likelihood")

# Scatter plot of residuals
# Data will be simulated to compute weighted residuals and npde
# the results shall be silently added to the object saemix.fit
# plot(saemix.fit,plot.type="residuals.scatter")

# Boxplot of random effects
# plot(saemix.fit,plot.type="random.effects")

# Relationships between parameters and covariates
# plot(saemix.fit,plot.type="parameters.vs.covariates")

# Relationships between parameters and covariates, on the same page
# par(mfrow=c(3,2))
# plot(saemix.fit,plot.type="parameters.vs.covariates",new=FALSE)

# VPC
# Not run (time constraints for CRAN)
# plot(saemix.fit,plot.type="vpc")

Description

This function will plot a longitudinal dataframe contained in an SaemixData object. By default it produces a spaghetti plot, but arguments can be passed on to modify this behaviour.

Usage

## S3 method for class 'SaemixData'
plot(x, y, ...)

## S3 method for class 'SaemixSimData'
plot(x, y, irep = -1, prediction = FALSE, ...)

Arguments

Details

this function can also be used to visualise the predictions for simulated values of the individual parameters, using the ypred element instead of the ysim element normally used here

Description

This function provides exploration plots for non Gaussian longitudinal data (work in progress, doesn't work yet for RTTE)

Usage

plotDiscreteData(object, outcome = "continuous", verbose = FALSE, ...)

plotDiscreteDataElement(
  object,
  outcome = "categorical",
  mirror = FALSE,
  irep = 1,
  verbose = FALSE,
  ...
)

Arguments

Details

This function is a very rough first attempt at automatically creating plots to explore discrete longitudinal data.

• for TTE data, a KM plot will be produced

• for count, categorical and binary data, a plot showing the proportion of each score/category across time will be shown These plots can be stratified over a covariate in the data set (currently only categorical covariates) by passing an argument which.cov='name' to the call #'

Author(s)

Emmanuelle Comets [email protected]

References

Brendel, K, Comets, E, Laffont, C, Laveille, C, Mentre, F. Metrics for external model evaluation with an application to the population pharmacokinetics of gliclazide, Pharmaceutical Research 23 (2006), 2036-2049.

Holford, N. The Visual Predictive Check: superiority to standard diagnostic (Rorschach) plots (Abstract 738), in: 14th Meeting of the Population Approach Group in Europe, Pamplona, Spain, 2005.

Ron Keizer, tutorials on VPC TODO

See Also

SaemixObject, saemix, saemix.plot.vpc, simulateDiscreteSaemix

Examples

# Time-to-event data
data(lung.saemix)

saemix.data<-saemixData(name.data=lung.saemix,header=TRUE,name.group=c("id"),
name.predictors=c("time","status","cens"),name.response=c("status"),
name.covariates=c("age", "sex", "ph.ecog", "ph.karno", "pat.karno", "wt.loss","meal.cal"),
units=list(x="days",y="",covariates=c("yr","","-","%","%","cal","pounds")))

# Plots a KM survival plot
plotDiscreteData(saemix.data, outcome="TTE")
# Plots a KM survival plot, stratified by sex
plotDiscreteData(saemix.data, outcome="TTE", which.cov="sex")

# Count data
data(rapi.saemix)
saemix.data<-saemixData(name.data=rapi.saemix, name.group=c("id"),
                 name.predictors=c("time","rapi"),name.response=c("rapi"),
                 name.covariates=c("gender"),units=list(x="months",y="",covariates=c("")))

# Plots a histogram of the counts
plotDiscreteData(saemix.data, outcome="count")

Description

Methods for function predict

Usage

## S4 method for signature 'SaemixObject'
predict(
  object,
  newdata = NULL,
  type = c("ipred", "ypred", "ppred", "icpred"),
  se.fit = FALSE,
  ...
)

Arguments

Value

a vector or a dataframe (if more than one type) with the corresponding predictions for each observation in the dataframe

Methods

list("signature(object = \"ANY\")")

Default predict functions

list("signature(object = \"SaemixObject\")")

Computes predictions using the results of an SAEM fit

Description

Predictions for a new dataset

Usage

## S3 method for class 'SaemixModel'
predict(object, predictors, psi = c(), id = c(), ...)

Arguments