10月1日（水）

ロジスティック回帰モデル

３Ｄプロットは僕の理解の範疇を超えているようである。

意味が分からなさすぎて発狂しそうになったので、諦めてperspでグラフをつくってみた。

age <- c(18,21,22,25,26,28,33,34,35,37,42,47,55,56,58,61,65,68,75,77)
out <- c(0,0,0,0,0,0,0,0,0,0,1,1,0,1,0,1,1,1,1,1)

d <- cbind(age,out)

age.n <- d[which(d[,2]==0),1]
age.o <- d[which(d[,2]==1),1]

p.n <- function(a,x0){
1-1/(1+exp(-a*(age.n-x0)))
}

p.o <- function(a,x0){
1/(1+exp(-a*(age.o-x0)))
}

p <- function(a,x0){
prod(p.n(a,x0))*prod(p.o(a,x0))
}

a. <- seq(0,0.5,0.01)
x0. <- 20:70

z <- c()
for(i in a.)
for(j in x0.)
z <- c(z,p(i,j))
z. <- matrix(z,ncol=length(x0.),byrow=T)

persp(a.,x0.,z.,theta = 30)

また、max(z.)を調べたところ、a=0.3,x0=34のときに最も尤度が高くなることが分かった。
実際に係数をこれにしてグラフを描くと、

a <- 0.3
x0 <- 34

f <- function(x){
1/(1+exp(-a*(x-x0)))
}

plot(age,out,xlim=c(15,77),ylim=c(0,1),ylab="リスク",main="年齢とリスクの関係")

par(new=T)
curve(f,15,77,xlim=c(15,77),ylim=c(0,1),ann=F)

なるほど確かに、それらしい曲線が引けている。

なお、対数をとって微分する方法もやってみたのだが、うまく解が求められなかったので挫折した。
一応そういう方法もある、ということだけ頭に入れておこうと思う。

回帰モデルによる交絡の制御

回帰モデルの強みは、例えば曝露群と非曝露群とで年齢層が全く異なっていて、層化が使えないような場合でも統計量の推定が行えることである。
（ただし、モデルが合っている保証はないが）

曝露群、非曝露群で全く年齢層が異なっているデータを作って、曝露による効果の推定を行ってみた。

ここで、アウトカム＝　k0 + k1*（年齢）+ k2*（曝露の有無1or0）
（k0,k1,k2は定数）

というモデル式を設定した。

k0 <-　0.2
k1 <-　0.3
k2 <-　2

N <- 10000
n <- 100

out. <- function(x){
k0+k1*x[1]+k2*x[2]+rnorm(1)
}

a1.d <- sample(20:80,N,replace=T)
a2.d <- sample(0:1,N,replace=T)
a3.d <- apply(cbind(a1.d,a2.d),1,out.)

D <- cbind(a1.d,a2.d,a3.d)

d <- data.frame(D)

d.e <- d[sample(which(d[,1]<=50&d[,2]==1),n),]
d.n <- d[sample(which(d[,1]>50&d[,2]==0),n),]

plot(d.e[,1],d.e[,3],xlim=c(20,80),ylim=c(0,30))

par(new=T)
plot(d.n[,1],d.n[,3],xlim=c(20,80),ylim=c(0,30))

result.e <- lm(d.e[,3]~d.e[,1],data=d)
result.n <- lm(d.n[,3]~d.n[,1],data=d)

E <- result.e$coefficients[1]-result.n$coefficients[1]
E

> E
(Intercept) 
   1.736189 

このように、Rに装備されているlm()の利用によって推定値を出すことは出来た。

ただ、この場合年齢による効果（k2）は曝露によらず一定とするモデルで考えているので、その条件を加えて推定を行った方がよいと思われる。
そこで、自分で最小二乗法を利用してやってみた。
（もしかしたらRに装備されている関数で出来るのかもしれないが）

k0 <-　0.2
k1 <-　0.3
k2 <-　3

N <- 10000
n <- 100

out. <- function(x){
k0+k1*x[1]+k2*x[2]+rnorm(1)
}

a1.d <- sample(20:80,N,replace=T)
a2.d <- sample(0:1,N,replace=T)
a3.d <- apply(cbind(a1.d,a2.d),1,out.)

d <- cbind(a1.d,a2.d,a3.d)

x <- sample(which(d[,1]<=50&d[,2]==1),n)
x. <- sample(which(d[,1]>50&d[,2]==0),n)

age <- d[x,1]
age. <- d[x.,1]

f1 <- function(k1,b1){
k1*age+b1
}

f2 <- function(k1,b2){
k1*age.+b2
}

P <- function(X){
sum((d[x,3]-f1(X[1],X[2]))^2)+sum((d[x.,3]-f1(X[1],X[3]))^2)
}

P(c(0.3,0.2,3.2))

op <- optim(c(0.3,0.2,3.2),P)

E <- op$par[3]-op$par[2]
E

> E
[1] 6.1551

理論的には合ってるような気がするのだが、値がうまく合わない。
ちょっと理由が分からないので、明日また検討することにする。

以上