foceiFit.Rd

% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/foceiFit.R
\name{foceiFit}
\alias{foceiFit}
\alias{focei.fit}
\alias{foceiFit.data.frame}
\alias{foceiFit.data.frame0}
\title{FOCEi fit}
\usage{
foceiFit(data, ...)

focei.fit(data, ...)

\method{foceiFit}{data.frame}(data, ...)

\method{foceiFit}{data.frame0}(
  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{
\item{data}{Data to fit; Needs to be RxODE compatible and have
\code{DV}, \code{AMT}, \code{EVID} in the dataset.}

\item{...}{Ignored parameters}

\item{inits}{Initialization list}

\item{PKpars}{Pk Parameters function}

\item{model}{The RxODE model to use}

\item{pred}{The Prediction function}

\item{err}{The Error function}

\item{lower}{Lower bounds}

\item{upper}{Upper Bounds}

\item{fixed}{Boolean vector indicating what parameters should be
fixed.}

\item{skipCov}{Boolean vector indicating what parameters should be
fixed when calculating covariances}

\item{control}{FOCEi options Control list.  See
\code{\link{foceiControl}}}

\item{thetaNames}{Names of the thetas to be used in the final
object.}

\item{etaNames}{Eta names to be used in the final object.}

\item{etaMat}{Eta matrix for initial estimates or final estimates
of the ETAs.}

\item{env}{An environment used to build the FOCEi or nlmixr
object.}

\item{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.}

\item{drop}{Columns to drop from the output}
}
\value{
A focei fit or nlmixr fit object

FOCEi fit object
}
\description{
FOCEi fit
}
\examples{

\donttest{

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

}
}
\author{
Matthew L. Fidler and Wenping Wang
}