R2WinBUGS - regresión logística con datos simulados

Question

Me pregunto si alguien tiene algún código R que use el paquete R2WinBUGS para ejecutar la regresión logística, idealmente con datos simulados para generar la 'verdad' y dos covariables continuas.

Gracias.

Cristiano

PS:

código potencial para generar datos artificiales (un caso) y dimensiones Winbugs ejecutar a través de r2winbugs (que no funciona todavía).

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

Su secuencia de comandos es exactamente la manera de hacerlo. Casi se está trabajando, sólo se requiere un cambio sencillo de hacer que funcione:

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

que define qué datos van a WinBugs. La variable C debe llenarse con verdadera. La presencia, N debe ser 1 de acuerdo con los datos generados. Tenga en cuenta que este es un caso especial de distribución binomial para N = 1, que se llama Bernoulli, un simple "lanzamiento de moneda".

Aquí está la salida:

> 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 

como ves, los parámetros se corresponden con los parámetros utilizados para generar los datos. Además, si lo compara con la solución frecuentadora, verá que corresponde.

EDIT: pero la típica regresión logística (~ binomial) mediría algunos conteos con algún valor superior N [i], y permitiría diferentes N [i] para cada observación. Por ejemplo, diga la proporción de juveniles a toda la población (N). Para ello sería necesario sólo para añadir índice a N en su modelo:

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

La generación de los datos sería algo como:

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) 

(final de edición)

Para más ejemplos de ecología de la población, ver books of Marc Kéry (Introducción a WinBUGS para ecologistas, y especialmente a Bayesian Population Analysis usando WinBUGS: una perspectiva jerárquica, que es un gran libro).

El guión completo que utiliza - el guión corregido de los suyos está aquí (la comparación con soluciones frequentist al final):

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