R2WinBUGS - логистическая регрессия с имитируемыми данными

Question

Мне просто интересно, есть ли у кого-нибудь код R, который использует пакет R2WinBUGS для запуска логистической регрессии - в идеале с имитируемыми данными для генерации «истины» и двух непрерывных совместных вариаций.

Спасибо.

Christian

PS:

Потенциальный код для создания искусственных данных (один одномерный случай) и запустить winbugs через r2winbugs (он пока не работает).

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 [i], и это позволило бы для разных N [i] для каждого наблюдения. Например, укажите долю несовершеннолетних для всего населения (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)