R2WinBUGS - 邏輯迴歸與模擬數據

Question

我只是想知道是否有人有一些R代碼使用包R2WinBUGS來運行邏輯迴歸 - 理想情況下用模擬數據來生成「真值」和兩個連續的協變量。

謝謝。

基督教

PS：

潛在的代碼來生成通過r2winbugs人工數據（一維的情況下），然後運行WinBUGS軟件（它不工作還沒有）。

library(MASS) 
library(R2WinBUGS) 

setwd("d:/BayesianLogisticRegression") 

n.site <- 150 

X1<- sort(runif(n = n.site, min = -1, max =1)) 

xb <- 0.0 + 3.0*X1 

occ.prob <- 1/(1+exp(-xb)) 

plot(X1, occ.prob,xlab="X1",ylab="occ.prob") 

true.presence <- rbinom(n = n.site, size = 1, prob = occ.prob) 

plot(X1, true.presence,xlab="X1",ylab="true.presence") 

# combine data as data frame and save 
data <- data.frame(X1, true.presence) 
write.matrix(data, file = "data.txt", sep = "\t") 

sink("model.txt") 
cat(" 
model { 

# Priors 
alpha ~ dnorm(0,0.01) 
beta ~ dnorm(0,0.01) 

# Likelihood 
for (i in 1:n) { 
    C[i] ~ dbin(p[i], N)  # Note p before N 
    logit(p[i]) <- alpha + beta *X1[i] 
} 
} 
",fill=TRUE) 
sink() 

# Bundle data 
win.data <- list(mass = X1, n = length(X1)) 

# Inits function 
inits <- function(){ list(alpha=rlnorm(1), beta=rlnorm(1))} 

# Parameters to estimate 
params <- c("alpha", "beta") 

# MCMC settings 
nc <- 3 #Number of Chains 
ni <- 1200 #Number of draws from posterior 
nb <- 200 #Number of draws to discard as burn-in 
nt <- 2 Thinning rate 

# Start Gibbs sampling 
out <- bugs(data=win.data, inits=inits, parameters.to.save=params, 
model.file="model.txt", n.thin=nt, n.chains=nc, n.burnin=nb, 
n.iter=ni, debug = TRUE) 

Answer 1

你的腳本就是這樣做的。這幾乎是工作，它只是需要一個簡單的改變，使其工作：

win.data <- list(X1 = X1, n = length(X1), C = true.presence, N = 1) 

其中定義的數據去WinBUGS軟件。變量C必須填寫true.presence，根據您生成的數據，N必須爲1 - 請注意，這是N = 1的二項分佈的特例，稱爲Bernoulli - 一個簡單的「硬幣翻轉」。

這裏是輸出：

> print(out, dig = 3) 
Inference for Bugs model at "model.txt", fit using WinBUGS, 
3 chains, each with 1200 iterations (first 200 discarded), n.thin = 2 
n.sims = 1500 iterations saved 
      mean sd 2.5%  25%  50%  75% 97.5% Rhat n.eff 
alpha  -0.040 0.221 -0.465 -0.187 -0.037 0.114 0.390 1.001 1500 
beta  3.177 0.478 2.302 2.845 3.159 3.481 4.165 1.000 1500 
deviance 136.438 2.059 134.500 135.000 135.800 137.200 141.852 1.001 1500 

你看，所述參數對應於用於生成所述數據的參數。另外，如果您與頻率主義者解決方案進行比較，您會發現它相對應。

EDIT：但典型物流（〜二項式）迴歸將測量與某些上限值N [1]一些計數，它會允許不同的N [1]對於每個觀測。例如說少年對全體人口的比例（N）。這就要求只是指數模型中加入N：

C[i] ~ dbin(p[i], N[i]) 

數據生成看起來是這樣的：

N = rpois(n = n.site, lambda = 50) 
juveniles <- rbinom(n = n.site, size = N, prob = occ.prob) 
win.data <- list(X1 = X1, n = length(X1), C = juveniles, N = N) 

（編輯結束）

對於來自更多的例子人口生態學見books of Marc Kéry（生態學家介紹WinBUGS，特別是使用WinBUGS的貝葉斯種羣分析：層次分析法，這是一本很棒的書）。

完整的劇本我用 - 你的修正後的腳本是列在這裏（末尾與頻率論的解決方案比較）：

#library(MASS) 
library(R2WinBUGS) 

#setwd("d:/BayesianLogisticRegression") 

n.site <- 150 

X1<- sort(runif(n = n.site, min = -1, max =1)) 

xb <- 0.0 + 3.0*X1 

occ.prob <- 1/(1+exp(-xb)) # inverse logit 

plot(X1, occ.prob,xlab="X1",ylab="occ.prob") 

true.presence <- rbinom(n = n.site, size = 1, prob = occ.prob) 

plot(X1, true.presence,xlab="X1",ylab="true.presence") 

# combine data as data frame and save 
data <- data.frame(X1, true.presence) 
write.matrix(data, file = "data.txt", sep = "\t") 

sink("tmp_bugs/model.txt") 
cat(" 
model { 

# Priors 
alpha ~ dnorm(0,0.01) 
beta ~ dnorm(0,0.01) 

# Likelihood 
for (i in 1:n) { 
    C[i] ~ dbin(p[i], N)  # Note p before N 
    logit(p[i]) <- alpha + beta *X1[i] 
} 
} 
",fill=TRUE) 
sink() 

# Bundle data 
win.data <- list(X1 = X1, n = length(X1), C = true.presence, N = 1) 

# Inits function 
inits <- function(){ list(alpha=rlnorm(1), beta=rlnorm(1))} 

# Parameters to estimate 
params <- c("alpha", "beta") 

# MCMC settings 
nc <- 3 #Number of Chains 
ni <- 1200 #Number of draws from posterior 
nb <- 200 #Number of draws to discard as burn-in 
nt <- 2 #Thinning rate 

# Start Gibbs sampling 
out <- bugs(data=win.data, inits=inits, parameters.to.save=params, 
model.file="model.txt", n.thin=nt, n.chains=nc, n.burnin=nb, 
n.iter=ni, 
working.directory = paste(getwd(), "/tmp_bugs/", sep = ""), 
debug = TRUE) 

print(out, dig = 3) 

# Frequentist approach for comparison 
m1 = glm(true.presence ~ X1, family = binomial) 
summary(m1) 

# normally, you should do it this way, but the above also works: 
#m2 = glm(cbind(true.presence, 1 - true.presence) ~ X1, family = binomial)