具有非整数权重的 JAGS 逻辑回归模型
JAGS logistic regression models with non-integer weights
对于下面的逻辑回归模型,我希望能够使用 n(和 y)的非整数值从后验中抽样。当部分数据可用或希望减轻重量时,这种模型可能会发生。
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 的整数值,因此 returns 在非整数 n 的情况下会出错。
我读到应该可以用这些技巧来做到这一点,但未能使其正常工作。感谢帮助。
正如你所说,你应该可以用那些技巧来做到这一点。困难的部分是正确编码二项式似然,因为 JAGS 没有二项式系数函数。然而,there are ways to do this. The model below should be able to do what you would like to. For a more specific explanation of the ones trick, 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[, ])
}
对于下面的逻辑回归模型,我希望能够使用 n(和 y)的非整数值从后验中抽样。当部分数据可用或希望减轻重量时,这种模型可能会发生。
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 的整数值,因此 returns 在非整数 n 的情况下会出错。
我读到应该可以用这些技巧来做到这一点,但未能使其正常工作。感谢帮助。
正如你所说,你应该可以用那些技巧来做到这一点。困难的部分是正确编码二项式似然,因为 JAGS 没有二项式系数函数。然而,there are ways to do this. The model below should be able to do what you would like to. For a more specific explanation of the ones trick, 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[, ])
}