Прогнозирование результатов с использованием JAGS от R

Question

[Код обновлен и больше не соответствует сообщениям об ошибках]

Я пытаюсь понять, как JAGS предсказывает значения результатов (для смешанной марковской модели). Я обучил модель на наборе данных, который включает в себя результат m и ковариаты x1, x2 а также x3,

Прогнозирование результата без фиксации значений параметров работает в R, но вывод кажется совершенно случайным:

preds <- run.jags("model.txt",
                  data=list(x1=x1, x2=x2, x3=x3, m=m,
                            statealpha=rep(1,times=M), M=M, T=T, N=N), monitor=c("m_pred"),
                  n.chains=1, inits = NA, sample=1)

Компиляция модели rjags... Вызов симуляции с использованием метода rjags... Примечание: модель не требует адаптации Запись в модели на 4000 итераций... |**************************************************| 100% Выполнение модели за 1 итерацию... Моделирование завершено Выполнено моделирование

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

preds <- run.jags("model.txt",
                  data=list(x1=x1, x2=x2, x3=x3,
                            statealpha=rep(1,times=M), M=M, T=T, N=N, beta1=beta1), monitor=c("m"),
                  n.chains=1, inits = NA, sample=1)

Компиляция модели rjags... Ошибка: Произошла следующая ошибка при компиляции и адаптации модели с использованием rjags: Ошибка в rjags::jags.model(модель, data = dataenv, n.chains = длина (runjags.object$end.state),: RUNTIME ERROR: ошибка компиляции в строке 39. beta1[2,1] является логическим узлом и его нельзя наблюдать

beta1 в этом случае матрица оценок коэффициента 2x2.

1. Как JAGS предсказывает m в первом примере (без фиксированных параметров)? Это просто совершенно случайно выбирая m?
2. Как я могу включить ранее полученные оценки модели для моделирования новых значений результатов?

Модель является:

model{
 for (i in 1:N)
 {

for (t in 1:T)
  {
  m[t,i] ~ dcat(ps[i,t,])
  }

for (state in 1:M)
  {
  ps[i,1,state] <- probs1[state]

  for (t in 2:T)
    {
    ps[i,t,state] <- probs[m[(t-1),i], state, i,t]
    }

  for (prev in 1:M){
       for (t in 1:T) {
    probs[prev,state,i,t] <- odds[prev,state,i,t]/totalodds[prev,i,t]
    odds[prev,state,i,t] <- exp(alpha[prev,state,i] +
                                beta1[prev,state]*x1[t,i]
                                + beta2[prev,state]*x2[t,i]
                               + beta3[prev,state]*x3[t,i])
    }}

  alpha[state,state,i] <- 0

      for (t in 1:T) {
  totalodds[state,i,t] <- odds[state,1,i,t] + odds[state,2,i,t]
  }
}
alpha[1,2,i] <- raneffs[i,1]
alpha[2,1,i] <- raneffs[i,2]
raneffs[i,1:2] ~ dmnorm(alpha.means[1:2],alpha.prec[1:2, 1:2])
}

for (state in 1:M)
  {
  beta1[state,state] <- 0
  beta2[state,state] <- 0
  beta3[state,state] <- 0
  }

beta1[1,2] <- rcoeff[1]
beta1[2,1] <- rcoeff[2]
beta2[1,2] <- rcoeff[3]
beta2[2,1] <- rcoeff[4]
beta3[1,2] <- rcoeff[5]
beta3[2,1] <- rcoeff[6]

alpha.Sigma[1:2,1:2] <- inverse(alpha.prec[1:2,1:2])
probs1[1:M] ~ ddirich(statealpha[1:M])
for (par in 1:6)
{
alpha.means[par] ~ dt(T.constant.mu,T.constant.tau,T.constant.k)
rcoeff[par] ~ dt(T.mu, T.tau, T.k)
}

T.constant.mu <- 0
T.mu <- 0
T.constant.tau <- 1/T.constant.scale.squared
T.tau <- 1/T.scale.squared
T.constant.scale.squared <- T.constant.scale*T.constant.scale
T.scale.squared <- T.scale*T.scale
T.scale <- 2.5
T.constant.scale <- 10
T.constant.k <- 1
T.k <- 1
alpha.prec[1:2,1:2] ~ dwish(Om[1:2,1:2],2)
Om[1,1] <- 1
Om[1,2] <- 0
Om[2,1] <- 0
Om[2,2] <- 1

## Prediction
for (i in 1:N)
    {

   m_pred[1,i] <- m[1,i]

    for (t in 2:T)  
      {
      m_pred[t,i] ~ dcat(ps_pred[i,t,])
      }

    for (state in 1:M)
      {
      ps_pred[i,1,state] <- probs1[state]

      for (t in 2:T)
        {
        ps_pred[i,t,state] <- probs_pred[m_pred[(t-1),i], state, i,t]
        }

      for (prev in 1:M)
        {

       for (t in 1:T)
         {
        probs_pred[prev,state,i,t] <- odds_pred[prev,state,i,t]/totalodds_pred[prev,i,t]
        odds_pred[prev,state,i,t] <- exp(alpha[prev,state,i] +
                                    beta1[prev,state]*x1[t,i]
                                    + beta2[prev,state]*x2[t,i]
                                   + beta3[prev,state]*x3[t,i])
        }}

      for (t in 1:T) {
      totalodds_pred[state,i,t] <- odds_pred[state,1,i,t] + odds_pred[state,2,i,t]
       }
      }
  }

Answer 1

TL;DR: я думаю, что вы просто упускаете вероятность.

Ваша модель сложна, поэтому, возможно, я что-то упустил, но, насколько я могу судить, вероятности нет. Вы поставляете предикторы x1, x2, а также x3 как данные, но вы не даете никаких наблюдений m, Итак, в каком смысле JAGS может "соответствовать" модели?

Чтобы ответить на ваши вопросы:

1. Да, похоже что m выбирается случайным образом из категориального распределения, обусловленного остальной частью модели. Так как нет m Если данные представлены в виде данных, ни у одного из распределений параметров нет причин для обновления, поэтому ваш результат для m ничем не отличается от того, что вы получили бы, если бы вы просто делали случайные ничьи со всех априоров и распространяли их через модель в R или как угодно.

2. Хотя это все равно не будет соответствовать модели в любом смысле, вы могли бы предоставить значения для beta1 если они еще не были полностью определены в модели. JAGS жалуется, потому что в настоящее время beta1[i] = rcoeff[i] ~ dt(T.mu, T.tau, T.k)и параметры к распределению T все фиксированы. Если какой-либо из (T.mu, T.tau, T.k) вместо этого были даны приоры (определяя их как случайные), затем beta1 могут быть предоставлены в качестве данных, и JAGS будет рассматривать rcoeff[i] ~ dt(T.mu, T.tau, T.k) как вероятность. Но в текущем виде модели, насколько это касается JAGS, если вы поставите beta1 как данные, это находится в противоречии с фиксированным определением уже в модели.

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

model {
  b ~ dnorm(0, 1) # prior on b

  for(i in 1:N) {
    y[i] ~ dnorm(b * x[i], 1) # Likelihood of y | b (and fixed precision = 1 for the example)
  }

  for(i in 1:N_pred) {
    pred_y[i] ~ dnorm(b * pred_x[i], 1) # Prediction
  }
}

В этом примере модели, x, y, а также pred_x предоставляются в виде данных, неизвестный параметр b должен быть оценен, и мы желаем последующие предсказания pred_y при каждом значении pred_x, JAGS знает, что распределение в первом for петля является вероятностью, потому что y предоставляется в качестве данных. Задние образцы b будет ограничен этой вероятностью. Второй for петля выглядит аналогично, но так как pred_y не предоставляется в качестве данных, он ничего не может сделать, чтобы ограничить b, Вместо этого JAGS знает, как просто рисовать pred_y образцы кондиционированы на b и поставляется pred_x, Значения pred_x обычно определяются как наблюдаемые xдавая прогнозный интервал для каждой наблюдаемой точки данных или в виде регулярной последовательности значений вдоль оси x, чтобы сгенерировать плавный прогнозный интервал.