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

Question

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

Благодарю.

Кристиан

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, который называется Бернулли - простой "бросок монеты".

Вот вывод:

> 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

Как видите, параметры соответствуют параметрам, используемым для генерации данных. Кроме того, если вы сравните с частым решением, вы увидите, что оно соответствует.

РЕДАКТИРОВАТЬ: но типичная логистическая (~ биномиальная) регрессия будет измерять некоторые подсчеты с некоторым верхним значением 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)

(конец редактирования)

Дополнительные примеры из экологии населения см. В книгах Марка Кери ("Введение в 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)