Невозможно выполнить байесовский анализ данных с помощью зазубрин в R

Question

Вот мой код

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] Узел несовместим с родителями

Answer 1

Вы начали свои начальные значения для параметров на границах ваших приоров. По сути, это не проблема, но эти lo git масштабированные значения, вероятно, являются крайними и поэтому создадут начальные оценки, которые являются Pr == 0 или Pr == 1.

. данные, давайте предположим, что мы инициализируем модель с beta0 = -100 и beta1 = 100.

Для вашей первой точки данных x = 30, так что ваш git -линейный предиктор начинается как:

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 в R. JAGS не является вентилятором прямого вывода из scale, чтобы вы в конечном итоге выполнили x = as.numeric(scale(c(30, 32, 34, 36, 38, 40, 42, 44))). Это полезно, поскольку это означает, что вы можете использовать стандартные приоры для своей регрессии, но означает, что вам нужно интерпретировать ваши коэффициенты по-разному. В Интернете много ресурсов взглянуть на относящиеся к среднему центру переменные.

Еще один способ получения начальных значений где-то в правильной области (при условии, что вы используете неопределенные априорные значения) - это просто подгонка логистов для частых c регрессионной модели и используйте эти оценки в качестве среднего некоторого случайного нормального распределения для каждого параметра. В конце концов, при неопределенных априорных оценках частота должна быть действительно близка к байесовской оценке, так как вероятность значительно превзойдет предыдущую.

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) довольно далеки от оценок параметров это специфицирует c модель.

Итак, если вы хотите, вы можете установить некоторые начальные значения следующим образом:

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)

Это, конечно, будет действительно работать, только если у вас нет предварительная информация и для очень простых моделей, таких как эта.