6

誰かがパッケージR2WinBUGSを使用してロジスティック回帰を実行するRコードを持っているかどうか疑問に思っています。理想的には、シミュレートされたデータを使用して「真実」と2つの連続共変量を生成します。

ありがとう。

キリスト教徒

PS：

人工データ（1次元の場合）を生成し、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)
4

1 に答える 1

5

あなたのスクリプトはまさにそれを行う方法です。それはほとんど機能しています、それを機能させるためにただ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)

（編集終了）

人口生態学のその他の例については、MarcKéryの本を参照してください （生態学者のための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)
于 2011-11-28T18:39:14.673 に答える