Фактор Байеса в R с Jaggs

Я пытаюсь выбрать между моделями линейной регрессии с использованием байесовского фактора в Jaggs через R. Для простоты я работаю над набором данных mtcars, а модели:

(1) милю на галлон только с перехватом (2) милю на галлон по весу (3) милю на галлон с дисплеем (4) милю на галлон по весу и дисплею

Я использую плоские априорные значения и пытаюсь реализовать сравнение моделей с помощью следующего кода, но сталкиваюсь с некоторыми проблемами и ошибками:

reg.bf.model <- "model{
    M ~ dcat(model.p[])
    model.p[1] <- 0.25
    model.p[2] <- 0.25
    model.p[3] <- 0.25
    model.p[4] <- 0.25

  # Likelihood
  for (i in 1:n){
      y[i] ~ dnorm(mu[i,M],tau[M])
        mu[i,1] <- alpha[1]
        mu[i,2] <- alpha[2] + beta1[1] * wt[i]
        mu[i,3] <- alpha[3] + beta1[2] * disp[i]
        mu[i,4] <- alpha[4] + beta1[3] * disp[i]+ beta2* wt[i]
  }

  for (m in 1:4){
      alpha[m] ~ dnorm(0,0.0001)
      sigma2[m] <- (1/tau[m])
  }

  for (j in 1:3) {
      beta1[j] ~ dnorm(0,0.0001)
  }

    beta2 ~ dnorm(0,0.0001)

    tau[1] ~ dgamma(prior.shape1[M],prior.rate1[M])
    tau[2] ~ dgamma(prior.shape2[M],prior.rate2[M])
    tau[3] ~ dgamma(prior.shape3[M],prior.rate3[M])
    tau[4] ~ dgamma(prior.shape4[M],prior.rate4[M])
}"

jags.data <- list(y=mtcars$mpg,n=nrow(mtcars), 
                  wt = mtcars$wt, disp = mtcars$disp,
                  prior.shape1=c(0.0001,1),prior.rate1=c(0.0001,1),
                  prior.shape2=c(1,0.0001),prior.rate2=c(1,0.0001),
                  prior.shape3=c(1,0.0001),prior.rate3=c(1,0.0001),
                  prior.shape4=c(1,0.0001),prior.rate4=c(1,0.0001))

jags.bf <- jags.model(file=textConnection(reg.bf.model),
                      #                      inits=jags.inits,
                      data=jags.data, n.chains=3)

update(jags.bf, 100)
jags.bf.out <- coda.samples(jags.bf,
                            variable.names=c("alpha","beta1","beta2"
                                             ,"sigma2","M"),
                            n.iter=50000, thin=400)

Но я получаю следующие сообщения об ошибках:

«Ошибка в update.jags (jags.bf, 100): ОШИБКА ЛОГИКИ: SimpleRange :: leftOffset. Индекс вне допустимого диапазона. Отправьте отчет об ошибке по адресу [email protected]enter code hereforge.net»

а также

"Ошибка в jags.samples (модель, имя_переменной, n.iter, thin, type =" trace ",: не установлены допустимые мониторы Дополнительно: Предупреждающие сообщения: 1: В FUN (X [[i]], .. .): Не могу установить монитор Нет модели!

2: В FUN (X [[i]], ...): невозможно установить монитор. Нет модели!

3: В FUN (X [[i]], ...): невозможно установить монитор. Нет модели!

4: В FUN (X [[i]], ...): невозможно установить монитор. Нет модели!

5: В FUN (X [[i]], ...): невозможно установить монитор. Нет модели! "

Мы очень ценим любые идеи или подталкивания в правильном направлении!


r bayesian jags modeling
person m1gnoc    schedule 18.02.2020    source источник
comment
Я думаю, что проблема связана с параметром модели M и его использованием. Должен признать, хотя у меня есть некоторый опыт работы с rjags, мне трудно понять, что делает строка M ~ dcat(model.p[]) и последующее использование M в качестве индекса. Во всяком случае, параметр модели prior.shape1 имеет только 2 элемента (в jags.data), и ошибка сообщает мне, что он пытался получить доступ к элементу [3] ...   -  person dario    schedule 18.02.2020

Ответы (1)


arrow_upward
0
arrow_downward

В качестве обновления, очевидно, фактор Байеса необходимо использовать в качестве альтернативы проверке гипотезы частотного анализа между двумя моделями. Поэтому следующая корректировка кода заставила все работать.

reg.bf.model <- "model{
  M ~ dcat(model.p[])
  model.p[1] <- prior1
  model.p[2] <- 1-prior1


  # Likelihood
  for (i in 1:n){
    y[i] ~ dnorm(mu[i,M],tau[M])
    mu[i,1] <- alpha[1]
    mu[i,2] <- alpha[2] + beta1[1] * wt[i]
  }

  for (m in 1:2){
      alpha[m] ~ dnorm(0,0.0001)
      sigma2[m] <- (1/tau[m])
  }

  for (j in 1:1) {
    beta1[j] ~ dnorm(0,0.0001)
  }

  beta2 ~ dnorm(0,0.0001)

  tau[1] ~ dgamma(prior.shape1[M],prior.rate1[M])
  tau[2] ~ dgamma(prior.shape2[M],prior.rate2[M])

}"
prior1 = 0.5
jags.data <- list(y=mtcars$mpg,n=nrow(mtcars), prior1 = prior1,
                  wt = mtcars$wt, disp = mtcars$disp,
                  prior.shape1=c(0.0001,1),prior.rate1=c(0.0001,1),
                  prior.shape2=c(1,0.0001),prior.rate2=c(1,0.0001))

jags.bf <- jags.model(file=textConnection(reg.bf.model),
                      #                      inits=jags.inits,
                      data=jags.data, n.chains=3)

update(jags.bf, 100)
jags.bf.out <- coda.samples(jags.bf,
                            variable.names=c("alpha","beta1","beta2"
                                             ,"sigma2","M"),
                            n.iter=50000, thin=400)

# calculating bayes factor
posterior.M <- unlist(jags.bf.out[,"M"])
posterior.odds <- mean(posterior.M==1)/mean(posterior.M==2)
prior.odds <- prior1/(1-prior1)
bayes.factor <- posterior.odds/prior.odds
bayes.factor
person m1gnoc    schedule 18.02.2020