r-statistics-fanの日記

統計好き人間の覚書のようなもの

RでTDM,AUC, C-max,Time above MIC

R

RでTDM,AUC, C-max,Time above MIC

 

一応rで薬剤の血中濃度を描くことができた。 C-maxやAUCやTime above MICも計算できる。 もっと正確に計算することは出来るだろうが、 正直この辺の薬剤濃度関連は誤差だらけなので、 正確性はあまり必要ないので、こんなもんで 良いかな。

あとは、各薬剤のパラメーターの設定だが、 難しい。なんというか、誤差が大きい。 添付文書のパラメータで自前で計算したのと、 専用のソフトで計算したのとで結果が 異なる。そもそも出力される曲線自体が 異なる。専用ソフトでも、それ、添付文書と 違うんじゃね？って感じの曲線。 しかしまあ、誤差といえば誤差なんだよね。 この辺のパラメター設定も将来の課題として おこう。しっかし、現実の測定データがない 素の段階でこれなので、測定データから、パラ メーターをベイズ推定するのも、誤差に飲み 込まれるんじゃね？ッて思うけど、いずれ やってみたい。

nmax <- 10
d <- c(rep(2000, nmax))  #dose mg
Ti <- c(rep(0.5, nmax))  #infusion time hour
cl <- c(rep(11.71, nmax))    #clrearance L/hour
ke <- c(rep(0.693 , nmax))  #ke  /hour
int <- 8   #infusion interval
time2 <- c(seq(from = 0, to = int * (10 - 1), by = int))

cpiv <- function(x, d = d[1], Ti = Ti[1], cl = cl[1], ke = ke[1]){
      d/Ti/cl * (1-exp(-ke*x))*(x < Ti) + (x >= Ti)*d/Ti/cl * (1-exp(-ke*Ti))*exp(-ke*(x-Ti))
}

cpn <- function(x, time2 = time2[1], d = d[1], Ti = Ti[1], cl = cl[1], ke = ke[1]){
      cpiv(x - time2, d = d, Ti = Ti, cl = cl, ke = ke)* (x >= time2)
} 

cp.sum <- function(x, nmax){
      res <- 0
      for (i in 1:nmax){
            res <- res + cpn(x, time2[i], d = d[i], Ti = Ti[i], cl = cl[i], ke = ke[i])  
      }
      return(res)
}

curve(cp.sum(x, nmax = nmax), 0.1, 80)#, log = "y", ylim = c(1, 100))

 

curve(cp.sum(x, nmax = nmax), 0.1, 80, log = "y", ylim = c(1, 100))

 

TAM <- function(mic = 8){
ct <- 0
num <- 0
for (i in seq(0.01, time2[nmax], by = 0.01)){
      ct <- ct + as.numeric(cp.sum(i, nmax = nmax) > mic)
      num <- num + 1
}
return(ct / num)
}

TAM(32)  #Time above MIC
## [1] 0.2522
integ <- function(x){cp.sum(x, nmax = nmax)}

AUC.calc <- function(){
temp <- numeric(nmax-1)
for (i in 1:(nmax-1)){
temp[i] <- integrate(integ, lower = time2[i], upper = time2[i+1])$value  #AUC
}
return(temp)
}
AUC.calc()
## [1] 170.0 170.8 170.8 170.8 170.8 170.8 170.8 170.8 170.8
xx <- seq(0.01, time2[nmax], by = 0.01)
plot(xx, cp.sum(xx, nmax = nmax),ylim = c(0, 110), xlim = c(0, time2[nmax]), type="l")

 

cp.sum(Ti + time2[1:nmax], nmax = nmax) #Cmax
##  [1] 100.0 100.4 100.4 100.4 100.4 100.4 100.4 100.4 100.4 100.4

今日はここまで