漸近線を考える

今やっていることで、用量反応曲線の漸近線の有無を評価したい。
毎度おなじみ、こちらのデータを使って漸近線について考えてみる。
　
用量反応曲線を書くために、R上では以下のロジスティック曲線を推定している。
　

　
微分係数が0のとき、傾きが0という高校数学の知識をもとに、与えられた濃度の両端で微分係数がどうなっているかをRでやってみる。
　
まずは関数を定義し、一階微分して数値を代入する……
　
その前に、用量反応曲線を推定するときに、濃度x軸を対数でなんやかんやしているので、微分をいじる（サンクスNK）。
　

　
ということで、をxで微分したあとにxをかけておけばよさそう。
　

# リンクのデータの最初だけやってみる
a<- 80/22  
conc    <- c(0.01,0.03,0.1, 0.3,1,  3,10, 30,100,300,1000) 
pro0   <- c(1,1,1,7,18,22,31,33,33,34,33)*a 
library(drc)
res<- drm(pro0~conc,fct=LL.4())
parm1<- res$parm[1,] # ヒル係数
parm2<- res$parm[2,] # 最小値
parm3<- res$parm[3,] # 最大値
parm4<- res$parm[4,] # EC50
df_drcddose<- deriv(~ parm2 + (parm3 - parm2)/((1 + exp(parm1*(log(dose/parm4))))),"dose")
# "dose"を変数と見て一階微分する。
dose<- c(min(conc),max(conc))
eval(df_drcddose) # doseを代入した時の値と微分した値を返す。

すべてについてやってみる。

a<- 80/22  
conc    <- c(0.01,0.03,0.1, 0.3,1,  3,10, 30,100,300,1000)  # Dose of isoproterenol
pro0   <- c(1,1,1,7,18,22,31,33,33,34,33)*a 
pro10  <- c(1,1,1,6,17,26,31,33,33,33,33)*a        
pro30  <- c(0,0,0,4,15,24,30,32,32,32,32)*a                 
pro100 <- c(0,1,1,1,11,21,29,32,33,33,33)*a                 
pro300 <- c(1,1,1,1,4,15,25,30,32,33,33)*a
pro1000<- c(1,1,1,1,1,5,16,26,31,32,33)*a
pro3000<- c(1,1,1,1,1,1,7,18,27,31,32)*a
pro.cont<- c(0,10,30,100,300,1000,3000)
EC50<- c(0.8442,0.9240,1.116,1.968,4.064,11.79,30.89)  

library(drc)
data<- rbind(pro0,pro10,pro30,pro100,pro300,pro1000,pro3000)
ldata<- matrix(0,nrow(data),4)
gradlist<- matrix(0,nrow(data),2) # 端点での傾きを収める。
par(mfrow=c(nrow(data),2))
for(p in 1:nrow(data)){
	res<- drm(data[p,]~conc,fct = LL.4())
	ldata[p,]<- res$start
	parm1<- res$parm[1,] # ヒル係数
	parm2<- res$parm[2,] # 最小値
	parm3<- res$parm[3,] # 最大値
	parm4<- res$parm[4,] # EC50
	df_drcddose<- deriv(~ parm2 + (parm3 - parm2)/((1 + exp(parm1*(log(dose/parm4))))),"dose")
	dose<- seq(min(conc),max(conc),by=0.05) # 細かくプロットしてみる。
	gradients<- attributes(eval(df_drcddose))$gradient*dose #attributes()を使って傾きを取り出す。変換した微分のためdoseをかけ忘れない。
	gradlist[p,1]<- head(gradients,n=1) # 左端での傾き
	gradlist[p,2]<- tail(gradients,n=1) # 右端での傾き
	plot(res,col=p,ylab="",main="log-logistic estimation")
	abline(v=parm4) # EC50の垂線
	plot(dose,attributes(eval(df_drcddose))$gradient*dose,type="l",log="x",col=p,
		 xlab="",ylab="gradient",main="gradient")
	abline(v=c(min(conc),parm4,max(conc))) # EC50と両端の垂線
}

gradlist

右が微分したもののプロット。垂線は濃度の端点とEC50である。
EC50で最大となっているのでよさそう。