雖然說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)))
