回帰係数のばらつき

• この続き
• 理想値からある程度のばらつきを入れてデータをシミュレーションし、回帰式推定をする
• 推定係数がどれくらいばらつくかも見てみよう

# シミュレーション
par(ask=TRUE)
Niter<-100
# 係数格納行列２つ

Km1<-Km2<-matrix(0,Niter,4)

for(ii in 1:Niter){

# 指数的に濃度を振る

minC<--2
maxC<-2
x<-c(2^(seq(from=minC,to=maxC,length=10)))
# Hill's eq のパラメタ
# a EC50
# b n
# U upper limit
# L lower limit
a<-2
b<-2
U<-5
L<-2

# データを作る
y<-U/(1+(a/x)^b)+L

# ばらつきを入れる
# 理想データの分散を基準にして、そのどれくらいかで指定する
s<-0.1
y<-y+rnorm(length(y),sd=sqrt(var(y))*s)

# 濃度の対数も作っておく
X<-log(x)

# optim()関数での推定をするための関数を２つ
# y=U*1/(1+(a/x)^b)+L
fr <- function(ks) {   
    k1 <- ks[1]
    k2 <- ks[2]
    k3 <- ks[3]
    k4 <- ks[4]
    sum((y-(k3/(1+(k1/x)^k2)+k4))^2)
    #sum((y-(x^k2/(k1^k2+x^k2)))^2)
}
# y=U*1/(1+exp(-b*(log(x)-log(a))))+L
fr2 <- function(ks) {   ## Rosenbrock Banana function
    k1 <- ks[1]
    k2 <- ks[2]
    k3 <- ks[3]
    k4 <- ks[4]
    sum((y-(k3/(1+exp(-k2*(log(x)-log(k1))))+k4))^2)
}
# optim()関数による推定
# 推定パラメタの初期値は適当にして推定する
opt.out<-optim(runif(4)*10, fr)
opt.out2<-optim(runif(4)*10, fr2)

Km1[ii,]<-opt.out$par
Km2[ii,]<-opt.out2$par

}

par(mfcol=c(2,2))
boxplot(Km1)
boxplot(Km2)
plot(Km1[,1],Km2[,1])
plot(Km1[,2],Km2[,2])


# plotする
par(mfcol=c(2,2))

plot(x,y)
t<-2^(seq(from=min(X),to=max(X),length=100))
points(t,opt.out$par[3]/(1+(opt.out$par[1]/t)^opt.out$par[2])+opt.out$par[4],col=2,cex=0.01)

plot(log(x),y)
#t<-seq(from=min(x),to=max(x),length=100)
points(log(t),opt.out$par[3]/(1+(opt.out$par[1]/t)^opt.out$par[2])+opt.out$par[4],col=2,cex=0.01)

plot(X,y)

points(log(t),opt.out2$par[3]/(1+exp(-opt.out2$par[2]*(log(t)-log(opt.out2$par[1]))))+opt.out2$par[4],col=2,cex=0.01)

plot(exp(X),y)

points(t,opt.out2$par[3]/(1+exp(-opt.out2$par[2]*(log(t)-log(opt.out2$par[1]))))+opt.out2$par[4],col=2,cex=0.01)
# 推定係数を見る
opt.out
opt.out2

hf1<-function(t){opt.out$par[3]/(1+(opt.out$par[1]/t)^opt.out$par[2])+opt.out$par[4]}
hf2<-function(t){opt.out2$par[3]/(1+exp(-opt.out2$par[2]*(log(t)-log(opt.out2$par[1]))))+opt.out2$par[4]}

integrate(hf1,0,4)
integrate(hf2,0,4)