From: Greg Tarpinian <sasprog474474_at_yahoo.com>

Date: Fri 21 Apr 2006 - 02:47:55 EST

}

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Fri Apr 21 02:59:06 2006

Date: Fri 21 Apr 2006 - 02:47:55 EST

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@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 Received on Fri Apr 21 02:59:06 2006

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.1.8, at Sat 27 May 2006 - 06:10:27 EST.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*