[R] nlminb( ) : one compartment open PK model

[Spencer Graves]

	  I've had pretty good luck solving problems like this when the post 
included a simple, reproducible example, as suggested in the posting 
guide (www.R-project.org/posting-guide.html).  Without that, I'm close 
to clueless.  I've encountered this kind of error myself before, and I 
think the way I solved it was just simplifying the example in steps 
until the problem went away.  By doing this, I could isolate the source 
of the problem.

	  Have you tried "trace=1", as described in the nlminb help page?  With 
this, nlminb prints "The value of the objective function and the 
parameters [at] every trace'th iteration."  Also, have you worked all 
the examples in the nlminb help page?  If that didn't show me how to 
solve my problem, I think I'd next try having nlminb estimate the 
gradient and hessian, rather than supplying them myself.

	  hope this helps,
	  spencer graves

Greg Tarpinian wrote:

> 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))
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html