このようなデータセットで
time C
0.1 2.6
0.25 4.817
0.5 6.596
0.75 6.471
1 6.049
1.5 5.314
2 4.611
2.5 4.5
3 4.392
4 4.013
5 3.698
6 3.505
8 3.382
12 2.844
14 2.383
24 1.287
このデータを以下のように定義されたモデルに当てはめたい
twocpt <- function(t, Cc, parms){
with(as.list(parms),{
dC0 <- -k01*C0
dCc <- k01*C0 + k21*Cp -(k12+ke)*Cc
dCp <- k12*Cc - k21*Cp
list(dCc)
})
}
このページ ( http://www.inside-r.org/packages/cran/FME/docs/modCost ) のリファレンスを参考にして、次のコードを作成しました。
#two compartment model, oral dosing
require(ggplot2)
require(FME)
require(XLConnect)
#Read Data from xlsx file, draw a scatter plot of the plasma-concentration profile
conc <- readWorksheetFromFile("E:/R/Book1.xlsx", sheet=1, header=TRUE)
pprofile <- ggplot(conc, aes(time, C))
pprofile <- pprofile + scale_x_continuous("Time (hr)")+scale_y_continuous("Concentration (ng/mL)")
(pprofile <- pprofile + geom_point()+geom_line())
#Create a matrix of the data frame.
concm <- as.matrix(conc)
#Define the parameters in the current simulation
k01 <- 1
k12 <- 10
k21 <- 0.5
ke <- 4
# wrap them up in the parms
parms <- c(k01=k01, k12=k12, k21=k21, ke=ke)
#Define the differential function
twocpt <- function(t, Cc, parms){
with(as.list(parms),{
dC0 <- -k01*C0
dCc <- k01*C0 + k21*Cp -(k12+ke)*Cc
dCp <- k12*Cc - k21*Cp
list(dCc)
})
}
#Define Cost function
Cost <- function(P) {
parms["k01"] <- P[1]
parms["k12"] <- P[2]
parms["k21"] <- P[3]
parms["ke"] <- P[4]
time <- conc[,1]
out <- ode(c(C=0), time, twocpt, parms)
return(modCost(out, concm))
}
Fit <- modFit(p=c(k01=10, k12=0.1, k21=0.4, ke=2), f=Cost)
summary(Fit)
ただし、次の警告メッセージが表示されました。
illegal input detected before taking any integration steps - see written message
どこに問題があるのか誰か教えてください。それともどうすれば手っ取り早く?ありがとう。