I have done more research and I am trying to pass values of parameters from lsoda into an objective function that I want to
minimize using optim function.
I am not having any sucess. My code is below annotated. oh also please excuse the terrbile coding.
My main question is how to pass parameter values back and forth from lsoda and optim?? ??
thanks z
#LOAD LIBRARIES #######################################################################
library(deSolve)
## A simple two compartment PK model#######################
# say liver and fat
#4 different doses
#input times, x and parms.
#x <-state varibables (need small starting values)
#parms <-parms? are these starting values? how to update?
#output list
twocompmodel <- function(t, x, parms) {
with(as.list(c(parms, x)), {
#3.3
Dose <- 3.3 #sigimp(t)
dCl <- Dose -(Ke+Klf)*Cl + Kfl*Cf #Liver
dCf <- Klf*Cl - Kfl*Cf #Fat
#7.3
Dose2 <- 7.3 #sigimp(t)
dCl2 <- Dose2 -(Ke+Klf)*Cl2 + Kfl*Cf2 #Liver
dCf2 <- Klf*Cl2 - Kfl*Cf2 #Fat
#15.2
Dose3 <- 15.2 #sigimp(t)
dCl3 <- Dose3 -(Ke+Klf)*Cl3 + Kfl*Cf3 #Liver
dCf3 <- Klf*Cl3 - Kfl*Cf3 #Fat
#33
Dose4 <- 33 #sigimp(t)
dCl4 <- Dose4 -(Ke+Klf)*Cl4 + Kfl*Cf4 #Liver
dCf4 <- Klf*Cl4 - Kfl*Cf4 #Fat
res <- c(dCl,dCf,dCl2,dCf2,dCl3,dCf3,dCl4,dCf4)
list(res)
})
}
## Parameters
parms <- c(Ke = 0.00194, Kfl = 0.015, Klf = .005)
## vector of timesteps
times <-c(0,91, 210, 364, 728 )
## Start values for steady state
y <- xstart <- c(Cl = .000001, Cf = 0.0001,Cl2 = .000001, Cf2 = 0.0001,Cl3 = .000001, Cf3 = 0.0001,Cl4 = .000001, Cf4 =
0.0001)
#FUNCTION TO MINIMIZE ##################################################################
#FUNCTION TO MINIMIZE ##################################################################
#FUNCTION TO MINIMIZE ##################################################################
#use output from lsoda and get RESIDUAL SUM OF SQUARES USING OBSERVED CONCENTRATIONS.
#function name "fitPBPK" input:: xstart,times, the twocompmodel function, parms, and observed conc in fat and liver.
#output SSR
fitPBPK<-function(xstart, times, twocompmodel, parms,obsfatconc,obsliverconc){
#HOW TO SPECIFY twocompmodel as a function # with the inputs times, xstart and parms???
OBSTcddFatConc<-obsfatconc
OBSTcddLiverConc<-obsliverconc
out<- as.data.frame(lsoda(xstart, times, twocompmodel, parms)) #OUTPUT OF LSODA
LiverConc<-c(rep(out$Cl[2:5],c(10,10,8,10)),rep(out$Cl2[2:5],c(10,10,8,10)),rep(out$Cl3[2:5],c(10,10,8,10)),rep(out$Cl4
[2:5],c(10,10,8,8)))
FatConc<-c(rep(out$Cf[2:5],c(10,10,8,10)),rep(out$Cf2[2:5],c(10,10,8,10)),rep(out$Cf3[2:5],c(10,10,8,10)),rep(out$Cf4[2:5],c
(10,10,8,8)))
SSR.fat<-sum((OBSTcddFatConc - FatConc)^2)
SSR.liver<-sum((OBSTcddLiverConc - LiverConc)^2)
sumofsqs<-SSR.fat+SSR.liver
return(sumofsqs)
}
#should I return PARMS as well?
####tests fitPBPK
fitPBPK(xstart, times, twocompmodel,c(Ke=.02, Kfl=.015 , Klf=.002 ),testdf$TcddFatConc,testdf$TcddLiverConc)
#[1] 2798343793 #works
##TRY TO USE WITH OPTIM#######################################################################
##TRY TO USE WITH OPTIM#######################################################################
##TRY TO USE WITH OPTIM#######################################################################
#ERROR!!!!
optim(fn=fitPBPK(xstart, times, twocompmodel,parms,testdf$TcddFatConc,testdf$TcddLiverConc),
par=c(Ke = 0.00194, Kfl = 0.015, Klf = .005), method = "BFGS", control=list(trace=6),hessian = T)
# Error in function (par) : could not find function "fn" ERROR!!!
####I know there are many things to fix but...
# how do I get the optim function to update the parms in the lsoda?
# how to pass parameter values back and forth from lsoda and optim?? ??
# S