R2WinBUGS - अनुरूपित डेटा के साथ लॉजिस्टिक रिग्रेशन

Question

मैं सोच रहा हूं कि किसी के पास कुछ आर कोड है जो लॉजिस्टिक रिग्रेशन चलाने के लिए पैकेज R2WinBUGS का उपयोग करता है - आदर्श रूप से 'सत्य' और दो निरंतर सह-भिन्नता उत्पन्न करने के लिए अनुरूपित डेटा के साथ।

धन्यवाद।

ईसाई

पुनश्च:

संभावित कोड कृत्रिम डेटा (एक आयामी मामले) और r2winbugs के माध्यम से चलाने के 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)) 

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 में जाते हैं। परिवर्तनीय सी को सत्य.presence से भरा जाना चाहिए, एन आपके द्वारा जेनरेट किए गए डेटा के अनुसार 1 होना चाहिए - ध्यान दें कि यह एन = 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 

जैसा कि आप देख, मानकों डेटा उत्पन्न करने के लिए इस्तेमाल किया मानकों के अनुरूप:

यहाँ उत्पादन होता है। इसके अलावा, यदि आप लगातार समाधान के साथ तुलना करते हैं, तो आप देखते हैं कि यह मेल खाता है।

संपादित: लेकिन ठेठ रसद (~ द्विपद) प्रतिगमन के साथ कुछ ऊपरी मूल्य एन [i] कुछ मायने रखता है मापने होता है और यह विभिन्न एन [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) 

(संपादित के अंत)

से अधिक उदाहरण के लिए:

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

डेटा पीढ़ी कुछ ऐसा दिखाई देगा: यह सिर्फ अपने मॉडल में एन के सूचकांक से जोड़ने के लिए की आवश्यकता होगी जनसंख्या पारिस्थितिकता books of Marc Kéry (पारिस्थितिक विज्ञानी के लिए विनब्यूजीएस का परिचय, और विशेष रूप से बेयसियन जनसंख्या विश्लेषण WinBUGS का उपयोग करके: एक पदानुक्रमित परिप्रेक्ष्य, जो एक महान पुस्तक है)।

पूरा स्क्रिप्ट मैं इस्तेमाल किया - तुम्हारा को सही स्क्रिप्ट यहां सूचीबद्ध है (अंत में frequentist समाधान के साथ तुलना): http://books.google.ca/books?id की

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