Conditional logistic regression

雖然說WinBUGS中有三種寫法, 但是實際上測試結果只有第一種寫法最穩定..

參考資料 Volumn 2. Endo: conditional inference in case-control studies

直接用例子來說明(這裡的例子不是WinBUGS的, 而是R裡頭的School例子)

model
{
for (i in 1:N) {
y[i,1] <- 1 y[i,2] <- 0 }

for (i in 1:N) {

# model 1. 用一般的logistic regression去做, 缺點是只能處理1-1 matching
# 結果(update 10K, 取後面5K, 每隔10取一點)
# (b1,b2,b3)=(1.079 0.968 , -0.471 0.7271, -1.795 0.7117)

y[i,1] ~ dbin(p[i,1],1)
logit(p[i,1]) <- b1*(smoking[i,1]-smoking[i,2]) + b2*(rubber[i,1]-rubber[i,2]) + b3*(alcohol[i,1]-alcohol[i,2])

# model 2. 用multinomial distribution去做, 缺點是容易發散..
# y[i,1:2] ~ dmulti(p[i,1:2],1)
# for (j in 1:2) {
# p[i,j] <- e[i,j]/sum(e[i,]) # e[i,j] <- max(0.00001, min( exp(b1*smoking[i,j]+b2*rubber[i,j]+b3*alcohol[i,j]) , 0.99999)) # log(e[i,j]) <- b1*smoking[2*(i-1)+j]+b2*rubber[2*(i-1)+j]+b3*alcohol[2*(i-1)+j] # }

# model 3. 使用poisson regression去做, 但是我始終做不出來...(sigh)
# for (j in 1:2) {
# Y[i,j] ~ dpois(mu2[i,j])
# log(mu[i,j]) <- b0[i] + b1*smoking[2*(i-1)+j]+b2*rubber[2*(i-1)+j]+b3*alcohol[2*(i-1)+j] # mu[i,j] <- exp(b0[i] + b1*smoking[i,j]+b2*rubber[i,j]+b3*alcohol[i,j]) # mu2[i,j] <- max(mu[i,j], 0.0001) # }
# b0[i] ~ dnorm(0,1)

}

b1 ~ dnorm(0,1.0E-6)
b2 ~ dnorm(0,1.0E-6)
b3 ~ dnorm(0,1.0E-6)
}

list(b1=0,b2=0,b3=0)

list(N=26,
smoking=structure(.Data=c(1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,0,1,1,1,1,0,0,1,1,1,1,0,1,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,1,0,1,1,1,1,0,0),.Dim=c(26,2)),
alcohol=structure(.Data=c(0,1,0,1,1,1,0,0,0,1,0,1,0,1,1,1,1,1,0,1,1,0,0,0,1,1,0,0,0,1,0,1,1,1,0,1,0,1,0,1,1,1,1,0,0,0,1,0,0,1,0,0),.Dim=c(26,2)),
rubber=structure(.Data=c(0,0,0,0,0,0,1,0,1,0,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,1,1,0,1,1,0,1,1,1,0,1,0,0,0,0,0,1,1,0,0,0),.Dim=c(26,2)))

model

{

for (i in 1:N) {

y[i,1] <- 1

y[i,2] <- 0

}

for (i in 1:N) {

y[i,1] ~ dbin(p[i,1],1)

logit(p[i,1]) <- b1*(smoking[i,1]-smoking[i,2]) + b2*(rubber[i,1]-rubber[i,2]) + b3*(alcohol[i,1]-alcohol[i,2])

# Y[i,1:2] ~ dmulti(p[i,1:2],1)

# for (j in 1:2) {

# p[i,j] <- e[i,j]/sum(e[i,])

# e[i,j] <- max(0.00001, min( exp(b1*smoking[i,j]+b2*rubber[i,j]+b3*alcohol[i,j]) , 0.99999))

# log(e[i,j]) <- b1*smoking[2*(i-1)+j]+b2*rubber[2*(i-1)+j]+b3*alcohol[2*(i-1)+j]

# }

# for (j in 1:2) {

# Y[i,j] ~ dpois(mu2[i,j])

# log(mu[i,j]) <- b0[i] + b1*smoking[2*(i-1)+j]+b2*rubber[2*(i-1)+j]+b3*alcohol[2*(i-1)+j]

# mu[i,j] <- exp(b0[i] + b1*smoking[i,j]+b2*rubber[i,j]+b3*alcohol[i,j])

# mu2[i,j] <- max(mu[i,j], 0.0001)

# }

# b0[i] ~ dnorm(0,1)

}

b1 ~ dnorm(0,1.0E-6)

b2 ~ dnorm(0,1.0E-6)

b3 ~ dnorm(0,1.0E-6)

}

list(b1=0,b2=0,b3=0)

list(N=26,

smoking=structure(.Data=c(1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,0,1,1,1,1,0,0,1,1,1,1,0,1,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,1,0,1,1,1,1,0,0),.Dim=c(26,2)),

alcohol=structure(.Data=c(0,1,0,1,1,1,0,0,0,1,0,1,0,1,1,1,1,1,0,1,1,0,0,0,1,1,0,0,0,1,0,1,1,1,0,1,0,1,0,1,1,1,1,0,0,0,1,0,0,1,0,0),.Dim=c(26,2)),

rubber=structure(.Data=c(0,0,0,0,0,0,1,0,1,0,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,1,1,0,1,1,0,1,1,1,0,1,0,0,0,0,0,1,1,0,0,0),.Dim=c(26,2)))