9月30日（火）

ロジスティックモデル

「ロスマンの疫学」p.303のデータにおいて、モデル式を

リスク＝

とし、最尤法による係数推定を行った。

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)

plot(age,out)

age.n <- deta[which(d[,2]==0),1]
age.o <- deta[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)

plot3d(a.,x0.,z.)

しかし、なんだか変な図が出来た上に、ブログへのデータの貼り方がよくわからない。

明日ちゃんとまとめようと思う。

以上