Почему выходные данные модели при использовании rjags и R2Jags отличаются?

Я работаю над подгонкой модели многоуровневой логистической регрессии с предикторами на уровне группы. Я использую JAGS через R. Я получаю другое поведение, когда подбираю модель с пакетами runjags и R2Jags.

Я попытался написать воспроизводимый пример, показывающий проблему. Ниже я моделирую данные из биномиальной модели, индексирую данные по 8 графикам и 2 блокам, а затем использую многоуровневую логистическую регрессию для восстановления вероятностей успеха (b1 и b2) в приведенном ниже коде. Прокрутите вниз, чтобы увидеть краткое описание двух вариантов.

У меня вопрос:

  1. Почему задние части этих двух посадок отличаются? Я использую одни и те же данные, одну спецификацию модели и устанавливаю генератор случайных чисел перед каждым из них. Почему различается среднее значение апостериоров и почему значения Rhat так различаются?
# -------------------------------------------------------------------
# Loading required packages
# -------------------------------------------------------------------
library(rjags) 
library(R2jags)
library(MCMCvis)

Информация о версии пакета:

jags.version()
[1] ‘4.3.0’

R2jags_0.5-7   MCMCvis_0.13.5 rjags_4-10

# -------------------------------------------------------------------
# Simulate data
# -------------------------------------------------------------------
set.seed(10)

N.plots = 8
N.blocks = 2
trials=400

n = rep(100,trials)
N=length(n)
plotReps=N/N.plots
blockReps=N/N.blocks

# Block 1
b1<-rep(c(.25,.75,.9,.1),each=plotReps)-.05
# Block 2
b2<-rep(c(.25,.75,.9,.1),each=plotReps)+.05

y = rbinom(trials, 100, p = c(b1,b2))

# vectors indexing plots and blocks
plot = rep(1:8,each=plotReps)
block = rep(1:2,each=blockReps)

# pass data to list for JAGS
data = list(
  y = y,
  n = n,
  N = length(n),
  plot = plot,
  block= block,
  N.plots = N.plots,
  N.blocks = N.blocks
)

# -------------------------------------------------------------------
# Code for JAGS model
# -------------------------------------------------------------------

modelString <- "model { 
  ## Priors

  # hyperpriors
  mu.alpha ~ dnorm(0, 0.0001)

  sigma.plot ~ dunif(0,100) 
  tau.plot <- 1 / sigma.plot^2

  sigma.block ~ dunif(0,100) 
  tau.block <- 1 / sigma.block^2

  # priors 
  for(i in 1:N.plots){     
    eps.plot[i]~dnorm(0,tau.plot)
  }

  for(i in 1:N.blocks){
    eps.block[i]~dnorm(0,tau.block)
  }

  # Likelihood
  for(i in 1:N){
    logit(p[i]) <- mu.alpha + eps.plot[plot[i]] + eps.block[block[i]]
    y[i] ~ dbin(p[i], n[i])

  }
}"

# -------------------------------------------------------------------
# Initial values
# -------------------------------------------------------------------
# set inits for rjags
inits = list(list(mu.alpha = 0,sigma.plot=2,sigma.block=2),
             list(mu.alpha = 0,sigma.plot=2,sigma.block=2),
             list(mu.alpha = 0,sigma.plot=2,sigma.block=2)) 

# set inits function for R2jags
initsFun<-function(){list(
  mu.alpha=0,
  sigma.plot=2,
  sigma.block=2
)}

# -------------------------------------------------------------------
# Set JAGS parameters and random seed
# -------------------------------------------------------------------
# scalars that specify the 
# number of iterations in the chain for adaptation
# number of iterations for burn-in
# number of samples in the final chain
n.adapt = 500
n.update = 5000
n.iterations = 1000
n.thin = 1
parsToMonitor = c("mu.alpha","sigma.plot","sigma.block","eps.plot","eps.block")

# -------------------------------------------------------------------
# Call to JAGS via rjags
# -------------------------------------------------------------------
set.seed(2)
# tuning (n.adapt)
jm = jags.model(textConnection(modelString), data = data, inits = inits,
                n.chains = length(inits), n.adapt = n.adapt)

# burn-in (n.update)
update(jm, n.iterations = n.update)

# chain (n.iter)
samples.rjags = coda.samples(jm, variable.names = c(parsToMonitor), n.iter = n.iterations, thin = n.thin)

# -------------------------------------------------------------------
# Call to JAGS via R2jags
# -------------------------------------------------------------------
set.seed(2)
samples.R2jags <-jags(data=data,inits=initsFun,parameters.to.save=parsToMonitor,model.file=textConnection(modelString),
                      n.thin=n.thin,n.chains=length(inits),n.burnin=n.adapt,n.iter=n.iterations,DIC=T)

# -------------------------------------------------------------------
# Summarize posteriors using MCMCvis
# -------------------------------------------------------------------
sum.rjags <- MCMCvis::MCMCsummary(samples.rjags,params=c("mu.alpha","eps.plot","sigma.plot","sigma.block","eps.block"))
sum.rjags

sum.R2jags2 <- MCMCvis::MCMCsummary(samples.R2jags,params=c("mu.alpha","eps.plot","sigma.plot","sigma.block","eps.block"))
sum.R2jags2

Вот результат подгонки rjags:

                     mean         sd         2.5%         50%       97.5% Rhat n.eff
mu.alpha      0.07858079 21.2186737 -48.99286669 -0.04046538 45.16440893 1.11  4063
eps.plot[1]  -1.77570813  0.8605892  -3.45736942 -1.77762035 -0.02258692 1.00  2857
eps.plot[2]  -0.37359614  0.8614370  -2.07913650 -0.37581522  1.36611635 1.00  2846
eps.plot[3]   0.43387001  0.8612820  -1.24273657  0.42332033  2.20253810 1.00  2833
eps.plot[4]   1.31279883  0.8615840  -0.38750596  1.31179143  3.06307745 1.00  2673
eps.plot[5]  -1.34317034  0.8749558  -3.06843578 -1.34747145  0.44451006 1.00  2664
eps.plot[6]  -0.40064738  0.8749104  -2.13233876 -0.41530587  1.37910977 1.00  2677
eps.plot[7]   0.36515253  0.8738092  -1.35364716  0.35784379  2.15597251 1.00  2692
eps.plot[8]   1.71826293  0.8765952  -0.01057452  1.70627507  3.50314147 1.00  2650
sigma.plot    1.67540914  0.6244529   0.88895789  1.53080631  3.27418094 1.01   741
sigma.block  19.54287007 26.1348353   0.14556791  6.68959552 93.21927035 1.22    94
eps.block[1] -0.55924545 21.2126905 -46.34099332 -0.24261169 48.81435107 1.11  4009
eps.block[2]  0.35658731 21.2177540 -44.65998407  0.25801739 49.31921639 1.11  4457

и вот результат подгонки R2jags:

                   mean         sd         2.5%         50%       97.5% Rhat n.eff
mu.alpha     -0.09358847 19.9972601 -45.81215297 -0.03905447 47.32288503 1.04  1785
eps.plot[1]  -1.70448172  0.8954054  -3.41749845 -1.70817566  0.08187877 1.00  1141
eps.plot[2]  -0.30070570  0.8940527  -2.01982416 -0.30458798  1.46954632 1.00  1125
eps.plot[3]   0.50295713  0.8932038  -1.20985348  0.50458106  2.29271214 1.01  1156
eps.plot[4]   1.37862742  0.8950657  -0.34965321  1.37627777  3.19545411 1.01  1142
eps.plot[5]  -1.40421696  0.8496819  -3.10743244 -1.41880218  0.25843323 1.01  1400
eps.plot[6]  -0.45810643  0.8504694  -2.16755579 -0.47087931  1.20827684 1.01  1406
eps.plot[7]   0.30319019  0.8492508  -1.39045509  0.28668886  1.96325582 1.01  1500
eps.plot[8]   1.65474420  0.8500635  -0.03632306  1.63399429  3.29585024 1.01  1395
sigma.plot    1.66375532  0.6681285   0.88231891  1.49564854  3.45544415 1.04   304
sigma.block  20.64694333 23.0418085   0.41071589 11.10308188 85.56459886 1.09    78
eps.block[1] -0.45810120 19.9981027 -46.85060339 -0.33090743 46.27709625 1.04  1795
eps.block[2]  0.58896195 19.9552211 -46.39310677  0.28183123 46.57874408 1.04  1769

Вот графики трассировки для mu.alpha из двух подборов. Во-первых, из rjags подойдут:

График трассировки для mu.alpha из rjags fit

Во-вторых, из R2jags подходят:

График трассировки для mu.alpha от R2Jags fit


r bayesian jags rjags r2jags
person gregor-fausto    schedule 19.02.2020    source источник
comment
Вы пробовали увеличить количество повторений в задней части? Вы проверили, сошлись ли они оба?   -  person Dason    schedule 19.02.2020
comment
Похоже, ваш параметр mu.alpha очень плохо ограничен данными. В итоге вы получите 95% центрального распределения в диапазоне от [-49,3, 51,7] в случае rjags и [-36,3, 46,1] в случае R2jags. Две разные аппроксимации выглядят так, как будто они сходятся в одном и том же апостериорном распределении. Я думаю, что любая разница связана с вариацией выборки в MCMC. Установка одного и того же случайного начального числа для каждого из них не гарантирует одинакового поведения, поскольку два пакета могут устанавливать вызов JAGS по-разному.   -  person qdread    schedule 19.02.2020
comment
Я только что добавил графики трассировки для mu.alpha, которые, как мне кажется, предполагают, что они не сходятся (?).   -  person gregor-fausto    schedule 19.02.2020
comment
Да, эти графики трассировки определенно предполагают, что они не сходятся   -  person qdread    schedule 19.02.2020
comment
Я увеличил количество повторений до 10000, и результаты задних упражнений более похожи между ними. Однако график трассировки для mu.alpha по-прежнему отбирает много пространства параметров, так что, возможно, проблема в том, что mu.alpha плохо ограничивается данными. Я пытался понять, как должна выглядеть модель случайных эффектов для многоуровневой логистики с помощью JAGS; Я ожидал, что это будет похоже на glmer (cbind (y, n) ~ 1 + (1 | plot) + (1 | block), family = binomial).   -  person gregor-fausto    schedule 19.02.2020
comment
Спасибо вам обоим. Я думаю, что моя проблема здесь на самом деле не связана с rjags или R2jags ... Как предположил @qdread, мой mu.alpha плохо ограничен данными. Думаю, это потому, что у меня всего 2 блока. Когда я запускаю это без block в модели, mu.alpha сходится.   -  person gregor-fausto    schedule 19.02.2020

Ответы (2)


arrow_upward
0
arrow_downward

Хотя часть проблемы связана с отсутствием конвергенции для mu.alpha, другая проблема заключается в том, как оба пакета определяют количество выборок для сбора из апостериорного распределения. Кроме того, вызов update после jags.model должен быть:

update(jm, n.iter = n.update)

вместо

update(jm, n.iterations = n.update)

Для rjags вы можете довольно легко указать количество шагов адаптации, шагов обновления и шагов итерации. Глядя на samples.rjags, совершенно ясно, что каждая цепочка имеет апостериор длиной n.iterations, всего (в этом примере) 3000 выборок (n.iterations * n.chains). И наоборот, R2jags::jags будет отбирать апостериорное число раз, равное аргументу n.iter минус аргумент n.burnin. Итак, как вы это указали, у вас 1) не включены n.update шагов в R2jags::jags и 2) производится выборка только апостериорных в общей сложности 1500 раз (каждая цепочка хранит только 500 выборок) по сравнению с 3000 раз из rjags.

Если вы хотите сделать аналогичное прожигание и сэмплировать такое же количество раз, вы могли бы вместо этого запустить:

samples.R2jags <-jags(
  data=data,
  inits=inits,
  parameters.to.save=parsToMonitor,
  model.file=textConnection(modelString),
  n.thin=n.thin,
  n.chains=length(inits),
  n.burnin=n.adapt + n.update ,
  n.iter=n.iterations +n.update + n.adapt,
  DIC=T
)

Наконец, R2jags по умолчанию загружает модуль glm, а rjags - нет. Это потенциально может вызвать некоторые расхождения, поскольку используемые сэмплеры, вероятно, разные (по крайней мере, в этом случае, потому что вы устанавливаете glm). Вы можете загрузить модуль glm с помощью rjags::load.module('glm') вызова перед вызовом jags.model.

И хотя это не связано с вопросом, по сути, я бы избегал использования i в ваших циклах for внутри модели для каждого цикла (используйте другую букву, если количество итераций варьируется между циклами):

modelString <- "model { 
  ## Priors

  # hyperpriors
  mu.alpha ~ dnorm(0, 0.0001)

  sigma.plot ~ dunif(0,100) 
  tau.plot <- 1 / sigma.plot^2

  sigma.block ~ dunif(0,100) 
  tau.block <- 1 / sigma.block^2

  # priors 
  for(i in 1:N.plots){     
    eps.plot[i]~dnorm(0,tau.plot)
  }

  for(j in 1:N.blocks){
    eps.block[j]~dnorm(0,tau.block)
  }

  # Likelihood
  for(k in 1:N){
    logit(p[k]) <- mu.alpha + eps.plot[plot[k]] + eps.block[block[k]]
    y[k] ~ dbin(p[k], n[k])

  }
}"
person mfidino    schedule 20.02.2020

arrow_upward
0
arrow_downward

Я почти уверен, что причина, по которой ваши апостериоры отличаются, заключается в том, что Jags не заботится о начальном наборе в R-коде.

Однако! Хотя set.seed() ничего не делает напрямую для Jags и бесполезен при вызове Jags напрямую через rjags, он распространяется, когда вы используете R2Jags.

Сравним:

  • rjags - это интерфейс нижнего уровня. Если вы не предоставите выбор случайного генератора и начальных чисел в ваших inits Jags, их инициализация будет основана на текущей временной метке.
  • R2Jags охватывает функции rjags. Функция jags() (R2Jags) вызывает jags.model() (rjags). Если вы проверите R-код функции jags(), вы увидите, что он генерирует начальное число для каждой цепочки, используя функцию runif() в R. Поскольку начальные числа Jags полагаются на вывод функции runif() в R, устанавливая начальное значение в R гарантирует, что вы будете получать одни и те же семена для Jags каждый раз, когда вы бежите.
person Philip De Vries    schedule 23.06.2020