# NOT RUN {
## 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 <- log(90) #log V (L)
# Bounds may be specified by c(lower, est, upper), like NONMEM:
# Residuals errors are assumed to be population parameters
prop.err <- c(0, 0.2, 1)
# 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.err)
})}
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 <- log(90) #log V (L)
## Bounds may be specified by c(lower, est, upper), like NONMEM:
## Residuals errors are assumed to be population parameters
prop.err <- c(0, 0.2, 1)
add.err <- c(0, 0.001)
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)
})}
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 <- log(90) #log V (L)
## Bounds may be specified by c(lower, est, upper), like NONMEM:
## Residuals errors are assumed to be population parameters
prop.err <- c(0, 0.2, 1)
delta <- c(0.1, 1, 10)
add.err <- c(0, 0.001)
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 <- .5 # log Ka
tcl <- -3.2 # log Cl
tv <- -1 # log V
eta.ka ~ 1
eta.cl ~ 2
eta.v ~ 1
add.err <- 0.1
})
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))
# }
Run the code above in your browser using DataLab