0

这是我的代码

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 格式的数据有关：

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，：节点 y[1] 中的错误节点与父节点不一致

4

1 回答 1

0

您已经在先验边界处开始了参数的初始值。这本身不是问题，但这些 logit 标度值可能会偏极端，因此会产生 Pr == 0 或 Pr == 1 的初始估计值。

只是处理您的数据，假设我们用beta0 = -100和初始化模型beta1 = 100。

对于您的第一个数据点x = 30，因此您的 logit 线性预测器开始为：

theta = -100 + 100 * 30

theta = 2900
plogis(theta) = 1

所以我们一开始的成功概率是 1，但是y=1对于n=20试验，成功概率不能是 1。你可以尝试一些事情来鼓励模型开始采样。

  1. 改变你的初始值。使它们更接近于 0（例如，介于 -4 和 4 之间）。
  2. 执行第 1 步，但还要重新调整您的x协变量，使其具有均值 = 0 和 sd = 1（即，使用.中的scale函数不是直接输出的粉丝，因此您最终会这样做。这很有帮助，因为这意味着您可以使用标准先验进行回归，但这意味着您需要以不同的方式解释您的系数。网上有很多资源可以查看与均值居中变量相关的资源。RJAGSscalex = 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)

当然，这只有在您没有先验信息并且对于非常简单的模型（例如这个模型）时才真正起作用。

于 2020-04-23T13:54:10.357 回答