2017-09-07 5 views
0

以下のロジスティック回帰モデルでは、n(およびy)の非整数値を使用して事後からサンプリングすることができます。これは、部分的なデータが利用可能である場合、または体重を下げることが望ましい場合に、この種のモデルで起こり得る。非整数加重を使用したJAGSロジスティック回帰モデル

model <- function() { 
    ## Specify likelihood 
    for (i in 1:N1) { 
     y[i] ~ dbin(p[i], n[i]) 
     logit(p[i]) <- log.alpha[1] + alpha[2] * d[i] 
    } 
    ## Specify priors 
    alpha[1] <- exp(log.alpha[1]) 
    alpha[2] <- exp(log.alpha[2]) 
    Omega[1:2, 1:2] <- inverse(p2[, ]) 
    log.alpha[1:2] ~ dmnorm(p1[], Omega[, ]) 
} 

DBINはnの整数値を必要とするので、非整数Nの場合にエラーを返します。

私は、これを1つのトリックで行うことは可能であるが、正しく動作させることができないと読んだ。ヘルプは高く評価しました。

答えて

1

あなたが言ったように、あなたは1つのトリックでこれを行うことができるはずです。 JAGSには二項係数関数がないため、困難な部分は二項尤度を正しくコーディングしています。しかし、there are ways to do this。以下のモデルは、あなたが望むことをすることができるはずです。それらのトリックのより具体的な説明については、see my answer hereを参照してください。

data{ 
    C <- 10000 
    for(i in 1:N1){ 
    ones[i] <- 1 
    } 
} 
model{ 
for(i in 1:N1){ 
# calculate a binomial coefficient 
bin_co[i] <- exp(logfact(n[i]) - (logfact(y[i]) + logfact(n[i] - y[i]))) 
# logit p 
logit(p[i]) <- log.alpha[1] + alpha[2] * d[i] 
# calculate a binomial likelihood using ones trick 
prob[i] <- (bin_co[i]*(p[i]^y[i])) * ((1-p[i])^(n[i] - y[i])) 
# put prob in Bernoulli trial and divide by large constant 
ones[i] ~ dbern(prob[i]/C) 
} 
## Specify priors 
alpha[1] <- exp(log.alpha[1]) 
alpha[2] <- exp(log.alpha[2]) 
Omega[1:2, 1:2] <- inverse(p2[, ]) 
log.alpha[1:2] ~ dmnorm(p1[], Omega[, ]) 
} 
+0

これは素晴らしいです!どうもありがとうございました。私はbin_coがアルファで定数なので不要だと考えるのは正しいですか? –

+0

いいえ、 'bin_co'は各データポイント(二項係数)とともに変化するので、非常に必要です。 'n'と' y'は各データ点で変化するので、 'bin_co'も同様です。 2番目の理由は、2項係数が2項尤度の一部であるためです。それを削除すると、分析では二項尤度を使用しなくなります。 –

関連する問題