こちらの続き。
データはこちらから。
曲線で囲まれる斜線部の面積を求めてみる。
という関数なので、与えられた濃度で積分してみる。
ただし、x軸は常用対数であることに注意して
conc<- c(0.01,0.03,0.1, 0.3,1, 3,10, 30,100,300,1000) resp<- c(1,1,1,7,18,22,31,33,33,34,33)*(80/22) library(drc) res<- drm(resp~conc,fct=LL.4()) parm1<- res$fit$par[1] # ヒル係数 parm2<- res$fit$par[2] # 最小値 parm3<- res$fit$par[3] # 最大値 parm4<- res$fit$par[4] # EC50 # 積分する関数 int_f<- function(x){ return((1/log(10))*(parm3-parm2)*(parm4^parm1)/((parm4^parm1+x^parm1)*x)) } # プロット用の関数 f<- function(x){return(parm2+(parm3-parm2)*(parm4^parm1)/(parm4^parm1+x^parm1))} # 積分する。 integrate(int_f,min(conc),max(conc)) dose<- c(seq(conc[1],conc[3],by=0.0001),seq(conc[3],max(conc),by=0.01)) plot(log(conc,10),f(conc),log="",type="n",xlab="",ylab="") par(new=TRUE) xval<- log(conc,10) dval<- f(conc) polygon(c(xval,rev(xval)),c(rep(parm2,length(xval)),rev(dval)),col=5,density=c(20,20)) par(new=TRUE) plot(log(conc,10),f(conc),log="",type="o",xlab="",ylab="") par(new=TRUE) plot(log(dose,10),f(dose),log="",type="l",col=2,lty=2,xlab="concentration",ylab="response") # 検算のため台形近似で面積計算をする。 S<- 0 for(i in 2:length(conc)){ S<- c(S,(resp[i-1]+resp[i])*(log(conc[i],10)-log(conc[i-1],10))/2) }
integrate(int_f,min(conc),max(conc))$value [1] 363.8412 sum(S) [1] 359.5129