[R] nlminb( ) : one compartment open PK model
Greg Tarpinian
sasprog474474 at yahoo.com
Thu Apr 20 18:47:55 CEST 2006
All,
I have been able to successfully use the optim( ) function with
"L-BFGS-B" to find reasonable parameters for a one-compartment
open pharmacokinetic model. My loss function in this case was
squared error, and I made no assumptions about the distribution
of the plasma values. The model appeared to fit pretty well.
Out of curiosity, I decided to try to use nlminb( ) applied to
a likelihood function that assumes the plasma values are normally
distributed on the semilog scale (ie, modeling log(conc) over
time). nlminb( ) keeps telling me that it has converged, but
the estimated parameters are always identical to the initial
values.... I am certain that I have committed "ein dummheit"
somewhere in the following code, but not sure what... Any help
would be greatly appreciated.
Kind regards,
Greg
model2 <- function(parms, dose, time, log.conc)
{
exp1 <- exp(-parms[1]*time)
exp2 <- exp(-parms[2]*time)
right.hand <- log(exp1 - exp2)
numerator <- dose*parms[1]*parms[2]
denominator <- parms[3]*(parms[2] - parms[1])
left.hand <- log(numerator/(denominator))
pred <- left.hand + right.hand
# defining the distribution of the values
const <- 1/(sqrt(2*pi)*parms[4])
exponent <- (-1/(2*(parms[4]^2)))*(log.conc - pred)^2
likelihood <- const*exp(exponent)
#defining the merit function
-sum(log(likelihood))
}
deriv2
<- deriv( expr = ~ -log(1/(sqrt(2*pi)*S)*exp((-1/(2*(S^2)))*
(log.conc-(log(dose*Ke*Ka/(Cl*(Ka-Ke)))
+log(exp(-Ke*time)-exp(-Ka*time))))^2)),
namevec = c("Ke","Ka","Cl","S"),
function.arg = function(Ke, Ka, Cl, S, dose, time, log.conc) NULL )
gradient2.1compart <- function(parms, dose, time, log.conc)
{
Ke <- parms[1]; Ka <- parms[2]; Cl <- parms[3]; S <- parms[4]
colSums(attr(deriv2.1compart(Ke, Ka, Cl, S, dose, time, log.conc), "gradient"))
}
attach(foo.frame)
inits <- c(Ke = .5,
Ka = .5,
Cl = 1,
S = 1)
#Trying out the code
nlminb(start = inits,
objective = model2,
gradient = gradient2,
control = list(eval.max = 5000, iter.max = 5000),
lower = rep(0,4),
dose = DOSE,
time = TIME,
log.conc = log(RESPONSE))
More information about the R-help
mailing list