foceiFit: FOCEi fit

In nlmixr: Nonlinear Mixed Effects Models in Population PK/PD

FOCEi fit

Description

FOCEi fit

Usage

foceiFit(data, ...)

focei.fit(data, ...)

## S3 method for class 'data.frame'
foceiFit(data, ...)

## S3 method for class 'data.frame0'
foceiFit(
  data,
  inits,
  PKpars,
  model = NULL,
  pred = NULL,
  err = NULL,
  lower = -Inf,
  upper = Inf,
  fixed = NULL,
  skipCov = NULL,
  control = foceiControl(),
  thetaNames = NULL,
  etaNames = NULL,
  etaMat = NULL,
  ...,
  env = NULL,
  keep = NULL,
  drop = NULL
)

Arguments

data	Data to fit; Needs to be RxODE compatible and have DV, AMT, EVID in the dataset.
...	Ignored parameters
inits	Initialization list
PKpars	Pk Parameters function
model	The RxODE model to use
pred	The Prediction function
err	The Error function
lower	Lower bounds
upper	Upper Bounds
fixed	Boolean vector indicating what parameters should be fixed.
skipCov	Boolean vector indicating what parameters should be fixed when calculating covariances
control	FOCEi options Control list. See foceiControl
thetaNames	Names of the thetas to be used in the final object.
etaNames	Eta names to be used in the final object.
etaMat	Eta matrix for initial estimates or final estimates of the ETAs.
env	An environment used to build the FOCEi or nlmixr object.
keep	Columns to keep from either the input dataset. For the input dataset, if any records are added to the data LOCF (Last Observation Carried forward) imputation is performed.
drop	Columns to drop from the output

Value

A focei fit or nlmixr fit object

FOCEi fit object

Author(s)

Matthew L. Fidler and Wenping Wang

Examples

## Comparison to Wang2007 objective functions

mypar2 = function ()
{
    k = theta[1] * exp(eta[1]);
}

mod <- RxODE({
    ipre = 10 * exp(-k * t)
})
pred <- function() ipre

errProp <- function(){
  return(prop(0.1))
}

inits <- list(THTA=c(0.5),
              OMGA=list(ETA[1] ~ 0.04));

w7 <- Wang2007

w7$DV <- w7$Y
w7$EVID <- 0
w7$AMT <- 0

## Wang2007 prop error OBF 39.458 for NONMEM FOCEi, nlmixr matches.
fitPi <- foceiFit(w7, inits, mypar2,mod,pred,errProp,
     control=foceiControl(maxOuterIterations=0,covMethod=""))

print(fitPi$objective)

## Wang2007 prop error OBF 39.207 for NONMEM FOCE; nlmixr matches.
fitP <- foceiFit(w7, inits, mypar2,mod,pred,errProp,
     control=foceiControl(maxOuterIterations=0,covMethod="",
     interaction=FALSE))

print(fitP$objective)

## Wang 2007 prop error OBF 39.213 for NONMEM FO; nlmixr matches
fitPfo <- foceiFit(w7, inits, mypar2,mod,pred,errProp,
     control=foceiControl(maxOuterIterations=0,covMethod="",
     fo=TRUE))

print(fitPfo$objective)

## Note if you have the etas you can evaluate the likelihood
## of an arbitrary model.  It doesn't have to be solved by
## FOCEi

etaMat <- matrix(fitPi$eta[,-1])

fitP2 <- foceiFit(w7, inits, mypar2,mod,pred,errProp, etaMat=etaMat,
      control=foceiControl(maxOuterIterations=0,maxInnerIterations=0,
      covMethod=""))


errAdd <- function(){
  return(add(0.1))
}

## Wang2007 add error of -2.059 for NONMEM FOCE=NONMEM FOCEi;
## nlmixr matches.
fitA <- foceiFit(w7, inits, mypar2,mod,pred,errAdd,
     control=foceiControl(maxOuterIterations=0,covMethod=""))

## Wang2007 add error of 0.026 for NONMEM FO; nlmixr matches

fitAfo <- foceiFit(w7, inits, mypar2,mod,pred,errAdd,
     control=foceiControl(maxOuterIterations=0,fo=TRUE,covMethod=""))

## Extending Wang2007 to add+prop with same dataset
errAddProp <- function(){
  return(add(0.1) + prop(0.1));
}

fitAP <- foceiFit(w7, inits, mypar2,mod,pred,errAddProp,
     control=foceiControl(maxOuterIterations=0,covMethod=""))

## Checking lognormal

errLogn <- function(){
   return(lnorm(0.1));
}

## First run the fit with the nlmixr lnorm error

fitLN <- foceiFit(w7, inits, mypar2,mod,pred,errLogn,
     control=foceiControl(maxOuterIterations=0,covMethod=""))


## Next run on the log-transformed space
w72 <- w7; w72$DV <- log(w72$DV)

predL <- function() log(ipre)

fitLN2 <- foceiFit(w72, inits, mypar2,mod,predL,errAdd,
     control=foceiControl(maxOuterIterations=0,covMethod=""))

## Correct the fitLN2's objective function to be on the normal scale
print(fitLN2$objective + 2*sum(w72$DV))

## Note the objective function of the lognormal error is on the normal scale.
print(fitLN$objective)

mypar2 <- function ()
{
    ka <- exp(THETA[1] + ETA[1])
    cl <- exp(THETA[2] + ETA[2])
    v  <- exp(THETA[3] + ETA[3])
}

mod <- RxODE({
    d/dt(depot) <- -ka * depot
    d/dt(center) <- ka * depot - cl / v * center
    cp <- center / v
})

pred <- function() cp

err <- function(){
    return(add(0.1))
}

inits <- list(THTA=c(0.5, -3.2, -1),
              OMGA=list(ETA[1] ~ 1, ETA[2] ~ 2, ETA[3] ~ 1));

## Remove 0 concentrations (should be lloq)

d <- theo_sd[theo_sd$EVID==0 & theo_sd$DV>0 | theo_sd$EVID>0,];

fit1 <- foceiFit(d, inits, mypar2,mod,pred,err)

## you can also fit lognormal data with the objective function on the same scale

errl <- function(){
    return(lnorm(0.1))
}

fit2 <- foceiFit(d, inits, mypar2,mod,pred,errl)

## You can also use the standard nlmixr functions to run FOCEi

library(data.table);
datr <- Infusion_1CPT;
datr$EVID<-ifelse(datr$EVID==1,10101,datr$EVID)
datr<-data.table(datr)
datr<-datr[EVID!=2]
datro<-copy(datr)
datIV<-datr[AMT>0][,TIME:=TIME+AMT/RATE][,AMT:=-1*AMT]
datr<-rbind(datr,datIV)

one.compartment.IV.model <- function(){
  ini({ # Where initial conditions/variables are specified
    # '<-' or '=' defines population parameters
    # Simple numeric expressions are supported
    lCl <- 1.6      #log Cl (L/hr)
    lVc <- 4.5  #log V (L)
    # Bounds may be specified by c(lower, est, upper), like NONMEM:
    # Residuals errors are assumed to be population parameters
    prop.sd <- 0.3
    # Between subject variability estimates are specified by '~'
    # Semicolons are optional
    eta.Vc ~ 0.1   #IIV V
    eta.Cl ~ 0.1; #IIV Cl
  })
  model({ # Where the model is specified
    # The model uses the ini-defined variable names
    Vc <- exp(lVc + eta.Vc)
    Cl <- exp(lCl + eta.Cl)
    # RxODE-style differential equations are supported
    d / dt(centr) = -(Cl / Vc) * centr;
    ## Concentration is calculated
    cp = centr / Vc;
    # And is assumed to follow proportional error estimated by prop.err
    cp ~ prop(prop.sd)
   })}

fitIVp <- nlmixr(one.compartment.IV.model, datr, "focei");

## You can also use the Box-Cox Transform of both sides with
## proportional error (Donse 2016)

one.compartment.IV.model <- function(){
ini({ # Where initial conditions/variables are specified
    ## '<-' or '=' defines population parameters
    ## Simple numeric expressions are supported
    lCl <- 1.6      #log Cl (L/hr)
    lVc <- 4.5  #log V (L)
    ## Bounds may be specified by c(lower, est, upper), like NONMEM:
    ## Residuals errors are assumed to be population parameters
    prop.err <- 0.3
    add.err <- 0.01
    lambda <- c(-2, 1, 2)
    zeta <- c(0.1, 1, 10)
    ## Between subject variability estimates are specified by '~'
    ## Semicolons are optional
    eta.Vc ~ 0.1   #IIV V
    eta.Cl ~ 0.1; #IIV Cl
})
model({ ## Where the model is specified
    ## The model uses the ini-defined variable names
    Vc <- exp(lVc + eta.Vc)
    Cl <- exp(lCl + eta.Cl)
    ## RxODE-style differential equations are supported
    d / dt(centr) = -(Cl / Vc) * centr;
    ## Concentration is calculated
    cp = centr / Vc;
    ## And is assumed to follow proportional error estimated by prop.err
    cp ~ pow(prop.err, zeta) + add(add.err) + boxCox(lambda)
    ## This is proportional to the untransformed f; You can use the transformed f by using powT()
})}

fitIVtbs <- nlmixr(one.compartment.IV.model, datr, "focei")

## If you want to use a variance normalizing distribution with
## negative/positive data you can use the Yeo-Johnson transformation
## as well.  This is implemented by the yeoJohnson(lambda) function.
one.compartment.IV.model <- function(){
ini({ # Where initial conditions/variables are specified
    ## '<-' or '=' defines population parameters
    ## Simple numeric expressions are supported
    lCl <- 1.6      #log Cl (L/hr)
    lVc <- 4.5  #log V (L)
    ## Bounds may be specified by c(lower, est, upper), like NONMEM:
    ## Residuals errors are assumed to be population parameters
    prop.err <- 0.3
    delta <- c(0.1, 1, 10)
    add.err <- 0.01
    lambda <- c(-2, 1, 2)
    ## Between subject variability estimates are specified by '~'
    ## Semicolons are optional
    eta.Vc ~ 0.1   #IIV V
    eta.Cl ~ 0.1; #IIV Cl
})
model({ ## Where the model is specified
    ## The model uses the ini-defined variable names
    Vc <- exp(lVc + eta.Vc)
    Cl <- exp(lCl + eta.Cl)
    ## RxODE-style differential equations are supported
    d / dt(centr) = -(Cl / Vc) * centr;
    ## Concentration is calculated
    cp = centr / Vc;
    ## And is assumed to follow proportional error estimated by prop.err
    cp ~ pow(prop.err, delta) + add(add.err) + yeoJohnson(lambda)
})}

fitIVyj <- nlmixr(one.compartment.IV.model, datr, "focei")

## In addition to using L-BFGS-B for FOCEi (outer problem) you may
## use other optimizers.  An example is below

one.cmt <- function() {
  ini({
      tka <- .44   # log Ka
      tcl <- log(c(0, 2.7, 100)) # log Cl
      tv <- 3.45    # log V
      eta.ka ~ 0.6
      eta.cl ~ 0.3
      eta.v ~ 0.1
      add.err <- 0.7
  })
  model({
      ka <- exp(tka + eta.ka)
      cl <- exp(tcl + eta.cl)
      v <- exp(tv + eta.v)
      linCmt() ~ add(add.err)
  })
}

fit <- nlmixr(one.cmt, theo_sd, "focei", foceiControl(outerOpt="bobyqa"))

## You may also make an arbitrary optimizer work by adding a wrapper function:

newuoa0 <- function(par, fn, gr, lower = -Inf, upper = Inf, control = list(), ...){
  ## The function requires par, fn, gr, lower, upper and control
  ##
  ## The par, fn, gr, lower and upper and sent to the function from nlmixr's focei.
  ## The  control is the foceiControl list
  ##
  ##  The following code modifies the list control list for no warnings.
  .ctl <- control;
  if (is.null(.ctl$npt)) .ctl$npt <- length(par) * 2 + 1
  ## nlmixr will print information this is to suppress the printing from the
  ## optimizer
  .ctl$iprint <- 0L;
  .ctl <- .ctl[names(.ctl) %in% c("npt", "rhobeg", "rhoend", "iprint", "maxfun")];
  ## This does not require gradient and is an unbounded optimization:
  .ret <- minqa::newuoa(par, fn, control=.ctl);
  ## The return statement must be a list with:
  ##    - x for the final parameter message
  ##    - message for a minimization message
  ##    - convergence for a convergence code
  .ret$x <- .ret$par;
  .ret$message <- .ret$msg;
  .ret$convergence <- .ret$ierr
  ## you can access the final list from the optimization by fit$optReturn
  return(.ret);
}

fit <- nlmixr(one.cmt, theo_sd, "focei", foceiControl(outerOpt=newuoa0))

nlmixr documentation built on March 27, 2022, 5:05 p.m.