Rで積分その2

MikuHatsune2011-11-06

こちらの続き。
データはこちらから。
 
曲線で囲まれる斜線部の面積を求めてみる。
\begin{align}f(x)&=E_{lower}+(E_{upper}-E_{lower})\frac{1}{1+e^{n(log{x}-log{EC50})}}\\&=E_{lower}+(E_{upper}-E_{lower})\frac{EC50^n}{EC50^n+x^n}\end{aling}
という関数なので、与えられた濃度x_{min} \leq x \leq x_{max}積分してみる。
ただし、x軸は常用対数であることに注意して
\begin{align}\int_{\hspace{20}x_{min}}^{\hspace{45}x_{max}} \{f(x)-E_{lower}\}d(log_{10}x)&=\int_{\hspace{20}x_{min}}^{\hspace{45}x_{max}} \{f(x)-E_{lower}\} \frac{d\frac{logx}{log10}}{dx}dx\\&=\frac{1}{log10}\int_{\hspace{20}x_{min}}^{\hspace{45}x_{max}} \{f(x)-E_{lower}\} (\frac{1}{x})dx\\&=\frac{1}{log10}\int_{\hspace{20}x_{min}}^{\hspace{45}x_{max}}(E_{upper}-E_{lower})\frac{EC50^n}{(EC50^n+x^n)x}dx\end{align}

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