Байесовский фактор в R с Джаггсом - PullRequest
Новая
       73

Байесовский фактор в R с Джаггсом

0 голосов
/ 18 февраля 2020

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

(1) миль на галлон только при перехвате (2) mpg на wt (3) mpg на disp (4) mpg на wt и disp

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

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): LOGI C ОШИБКА: SimpleRange :: leftOffset. Index вне допустимого диапазона Пожалуйста, отправьте отчет об ошибке в martyn_plummer@users.sourceenter code here forge. net "

и

" Ошибка в jags.samples (model, variable.names, 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]], ...): Не удается установить монитор. Никакая модель! "

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

1 Ответ

0 голосов
/ 18 февраля 2020

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

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
...