这是我的代码
library(R2jags) #library(rjags)
library(bayesplot)
library(coda)
# set working directory
setwd("/Users/isa/Desktop/logreg")
# BUGS model code
cat("model {
for( i in 1 : 8 ) {
y[i] ~ dbin(theta[i],n[i])
logit(theta[i]) <- beta0 + beta1 * x[i]
}
beta0 ~ dunif(-100, 100)
beta1 ~ dunif(-100, 100)
}",
file = "model_log.txt")
data <- read.delim("data.txt",
sep = "",
header = TRUE,
check.names = "FALSE",
stringsAsFactors = FALSE)
initsone <- list(beta0 = -100, beta1 = 100)
initstwo <- list(beta0 = 100, beta1 = -100)
initslog <- list(initsone, initstwo)
paramslog <- c("beta0", "beta1", "theta[6]")
outputlog <-
jags(data = data,
inits = initslog,
parameters.to.save = paramslog,
model.file = "model_log.txt",
n.chains = 2,
n.iter = 1000,
n.burnin = 1000,
n.thin = 1,
DIC = TRUE#,
# bugs.directory = getwd(),
# working.directory = getwd()
)
一切都很好，直到我尝试编译输出。我得到一个错误。
Error in jags.model(model.file, data = data, inits = init.values, n.chains = n.chains, :
Error in node y[1]
Node inconsistent with parents
我相信这和我的数据有关，我的数据是OpenBugs格式的，但我把它转换成了R格式：
list(y = c(1, 3, 6, 8, 11, 15, 17, 19),
n = c(20, 20, 20, 20, 20, 20, 20, 20),
x = c(30, 32, 34, 36, 38, 40, 42, 44),
N = 8 )
但我把它转换成了R格式
y n x
1 20 30
3 20 32
6 20 34
8 20 36
11 20 38
15 20 40
17 20 42
19 20 44
我是不是把数据转换错了？数据中哪里出了问题？一切都很好，直到我尝试编译输出。我得到一个错误说明。在jags.model(model.file, data = data, inits = init.values, n.chains = n.chains, : Error in node y[1]Node inconsistent with parents.
你已经在你的前值的边界开始了你的参数的初始值。这本身并不是一个问题，但这些对数比例值很可能太过极端，因此会产生Pr == 0或Pr == 1的初始估计。
仅仅通过您的数据，让我们假设我们初始化模型为 beta0 = -100 和 beta1 = 100.
对于您的第一个数据点 x = 30 所以你的对数线性预测器开始是。
theta = -100 + 100 * 30
theta = 2900
plogis(theta) = 1
所以我们开始时的成功概率是1，但是... ... y=1 对于 n=20 试验，所以成功概率不能为1。你可以尝试一些事情来鼓励模型开始采样。
x 协变量的均值=0，sd=1（即使用 scale 作用于 R. JAGS 是不喜欢直接从 scale 所以你最终会做 x = as.numeric(scale(c(30, 32, 34, 36, 38, 40, 42, 44))). 这很有帮助，因为这意味着你可以使用标准的前值来进行回归，但意味着你需要以不同的方式解释你的系数。网上有很多资源可以看与均值中心变量相关的资料。另一种在正确区域的某个地方获得初始值的方法（假设你使用的是模糊前值）是只需拟合一个频繁主义逻辑回归模型，并将这些估计值作为每个参数的一些随机正态分布的平均值。毕竟，在模糊前值的情况下，频繁主义估计应该非常接近贝叶斯估计，因为可能性会大大超过前值。
dat <- list(y = c(1, 3, 6, 8, 11, 15, 17, 19),
n = c(20, 20, 20, 20, 20, 20, 20, 20),
x = c(30, 32, 34, 36, 38, 40, 42, 44),
N = 8 )
# set up as matrix of successes and failures
y <- matrix(NA, ncol = 2, nrow = 8)
y[,1] <- dat$y
y[,2] <- dat$n - dat$y
m1 <- glm(y ~ dat$x, family = binomial)
summary(m1)
Call:
glm(formula = y ~ dat$x, family = binomial)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.39289 -0.20654 -0.04323 0.21294 0.50657
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -13.55295 2.05832 -6.584 4.57e-11 ***
dat$x 0.36630 0.05536 6.616 3.68e-11 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 70.7352 on 7 degrees of freedom
Residual deviance: 0.7544 on 6 degrees of freedom
AIC: 27.71
Number of Fisher Scoring iterations: 4
在这里你可以看到，数值在 abs(100) 与这个特定模型的参数估计值相差甚远。
所以，如果你愿意，你可以这样设置一些初始值。
initsone <- list(
beta0 = rnorm(1, m1$coefficients[1], 2),
beta1 = rnorm(1, m1$coefficients[2], 2)
)
initstwo <- list(
beta0 = rnorm(1, m1$coefficients[1], 2),
beta1 = rnorm(1, m1$coefficients[2], 2)
)
initslog <- list(initsone, initstwo)
当然，这只有在没有预先信息的情况下 才会真正起作用，而且也只适用于非常简单的模型，比如这个模型。