- こちらで、lm()関数を使って線形回帰をしている
- 最後に以下の一行を付けると、回帰直線がついた図も得られる
plot(lm(log(pro.cont[c(2:length(pro.cont))],10)~log((r-1),10)))
- 線形回帰は、一番単純な回帰
- これをもう少し面倒くさくするべく、変化させて行くと、glm()となる
- 近似関数を得るにはoptim()関数もある。以下の図は、ryamada本"R13-2.3R"
- glm()はこちらやこちら
- 「回帰」はこちらの記事の筆者とその後続ができるようになろうとしている処理
- 全部を通したコマンドを、便利のために末尾につけておくこととする
a<- 80/22 # 印刷されたものを測って手で変換するというなんというハイテク
#cont<- c(-2, -1.5, -1,-0.5,0,0.5, 1,1.5, 2,2.5, 3)
cont<- 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 # Score meant mean blood pressure difference
pro10 <- c(1,1,1,6,17,26,31,33,33,33,33)*a # on printed paper.
pro30 <- c(0,0,0,4,15,24,30,32,32,32,32)*a # a converted score into real blood pressure difference.
pro100 <- c(0,1,1,1,11,21,29,32,33,33,33)*a # Dose of propranolol was from 0 to 3000.
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)
M<- max(pro0,pro10,pro30,pro100,pro300,pro1000,pro3000)
plot(cont,pro0 ,log="x",type="l",col=2,ylim=c(0,M),ylab="",xlab="")
par(new=T)
plot(cont,pro10 ,log="x",type="l",col=3,ylim=c(0,M),ylab="",xlab="")
par(new=T)
plot(cont,pro30 ,log="x",type="l",col=4,ylim=c(0,M),ylab="",xlab="")
par(new=T)
plot(cont,pro100 ,log="x",type="l",col=5,ylim=c(0,M),ylab="",xlab="")
par(new=T)
plot(cont,pro300 ,log="x",type="l",col=6,ylim=c(0,M),ylab="",xlab="")
par(new=T)
plot(cont,pro1000,log="x",type="l",col=7,ylim=c(0,M),ylab="",xlab="")
par(new=T)
plot(cont,pro3000,log="x",type="l",col=8,ylim=c(0,M),
ylab="?Blood pressure [mmHg]",xlab="Isoproterenol concentration [nM/kg]")
legend(100,40,c(0,10,30,100,300,1000,3000),col=c(2:8),lty=1)
text(270,43,"propranolol")
title("Isoproterenol dose response curve")
EC50<- c(0.8442,0.9240,1.116,1.968,4.064,11.79,30.89) # 反応が50%となるアゴニスト濃度
r<- EC50/EC50[1] # 容量比というものらしい
r<- r[-c(1)]
plot(pro.cont[c(2:length(pro.cont))],(r-1),log="xy",
xlab="log([propranolol])",
ylab="log(r-1)")
abline(h=0)
# 直線近似をします。
lm(log(pro.cont[c(2:length(pro.cont))],10)~log((r-1),10))
plot(lm(log(pro.cont[c(2:length(pro.cont))],10)~log((r-1),10)))
