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

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

6 голосов
/ 24 ноября 2011

Мне просто интересно, есть ли у кого-нибудь 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)

1 Ответ

5 голосов
/ 28 ноября 2011

Твой сценарий - это как раз то, что нужно.Он почти работает, для его работы требуется всего одно простое изменение: 

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

, которое определяет, какие данные отправляются в WinBugs.Переменная C должна быть заполнена значением true.presence, N должно быть 1 в соответствии с сгенерированными вами данными - обратите внимание, что это особый случай биномиального распределения для N = 1, который называется Бернулли - простой"coin flip".

Вот вывод: 

> 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)

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

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