R2WinBUGS - regressione logistica con dati simulati

Question

Mi chiedo se qualcuno abbia qualche codice R che utilizza il pacchetto R2WinBUGS per eseguire la regressione logistica, idealmente con dati simulati per generare la "verità" e due continue co-variate.

Grazie.

Christian

PS:

codice potenziale per generare i dati artificiali (un caso) e dimensionale winbugs eseguito tramite r2winbugs (non funziona ancora).

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

Il tuo script è esattamente il modo per farlo. E 'quasi funzionando, è solo necessario un cambiamento semplice per farlo funzionare:

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

che definisce quali dati andare WinBugs. La variabile C deve essere riempita con true.presence, N deve essere 1 in base ai dati generati - si noti che questo è un caso speciale di distribuzione binomiale per N = 1, che è chiamato Bernoulli - un semplice "coin flip".

Ecco l'output:

> 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 

come si vede, i parametri corrispondono ai parametri utilizzati per generare i dati. Inoltre, se si confronta con la soluzione frequentista, si vede che corrisponde.

EDIT: ma la regressione logistica (~ binomiale) tipica misurerebbe alcuni conteggi con un valore superiore N [i] e consentirebbe N differenti per ciascuna osservazione. Ad esempio, diciamo la proporzione di giovani per l'intera popolazione (N). Ciò richiederebbe solo per aggiungere indice N nel modello:

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

La generazione di dati sarebbe qualcosa di simile:

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) 

(fine della modifica)

Per ulteriori esempi da ecologia della popolazione vedi books of Marc Kéry (Introduzione a WinBUGS per ecologisti, e soprattutto Analisi della popolazione bayesiana usando WinBUGS: una prospettiva gerarchica, che è un grande libro).

Lo script completo che ho usato - lo script corretto dei tuoi elencati qui (confronto con la soluzione frequentista alla fine):

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