Блог

Я установил следующую simple linear regression Bayesian модель, используя rjags . Я смог запустить модель, указав все предикторы отдельно (например, для lm объекта). Теперь я хочу узнать, как указать предсказатели, введя их в виде матрицы вместо того, чтобы указывать их по отдельности.

Итак, я запустил следующий код, но он выдал несколько ошибок.

Я использовал tobbaco набор данных в rrr пакете, чтобы предоставить воспроизводимый пример.

 library(rrr)
require(dplyr)
library(rjags)
tobacco <- as_data_frame(tobacco)
N1 = length(tobacco$Y1.BurnRate)
x1 = model.matrix(Y1.BurnRate~X2.PercentChlorine X3.PercentPotassium ,data = tobacco)
        
bayes_model_mul1=
 "model {
      for(i in 1:N1){
        Y1.BurnRate[i]~dnorm(mu1[i],tau1)
        
        for(j in 1:3){
          mu1[i]=beta1[j]*x1[i,j]  
        }
      }
      
      for (l in 1:3) { beta1[l] ~dnorm(0, 0.001) }        
      tau1 ~ dgamma(.01,.01)
      sigma_tau1 = 1/tau1 
        
    }"
        
        
model3 <- jags.model(textConnection(bayes_model_mul1), 
                     data = list(Y1.BurnRate=tobacco$Y1.BurnRate, x1=x1, N1=N1),
                     n.chains=1)
  

После запуска model3 я получил следующую ошибку.

Ошибка в jags.model(textConnection(bayes_model_mul1), data = list(Y1.BurnRate = tobacco $Y1.BurnRate, :
ОШИБКА ВРЕМЕНИ ВЫПОЛНЕНИЯ:
ошибка компиляции в строке 6.
Попытка переопределить узел mu1[1]

Кто-нибудь может помочь мне разобраться в этом? Связано ли это с введением предикторов в виде матрицы?

Есть несколько способов сделать это, вот два:

1. Используйте умножение матрицы вне цикла правдоподобия

 m1 =
 "model {
     mu1 = x1 %*% beta1 # ---> this
     for(i in 1:N1){
       Y1.BurnRate[i] ~ dnorm(mu1[i], tau1)
    }

    for (l in 1:3) { beta1[l] ~ dnorm(0, 0.001) }
    tau1 ~ dgamma(.01,.01)
    sigma_tau1 = 1/tau1 
 }"
  

2. Используйте inprod для умножения параметров на матрицу проектирования

 m2 = 
  "model {
    for(i in 1:N1){
      mu1[i] = inprod(beta1, x1[i,]) #----> this
      Y1.BurnRate[i] ~ dnorm(mu1[i], tau1)
    }

    for (l in 1:3) { beta1[l] ~ dnorm(0, 0.001) }
    tau1 ~ dgamma(.01,.01)
    sigma_tau1 = 1/tau1 
 }"
  

Вы получали ошибку с for(j in 1:3){ mu1[i] = beta1[j]* x1[i,j] } тем, что каждый раз, когда вы выполняете цикл, хотя индекс параметра j вы перезаписываете mu1[i] . Это также не суммирует отдельные термины. Возможно, вы также сможете индексировать mu1 с помощью j , а затем sum но не протестированы…