ошибка rjags Недопустимый аргумент вектора для ilogit

Я хотел бы сравнить регрессию betareg с той же регрессией с использованием rjags.

library(betareg)
d = data.frame(p= sample(c(.1,.2,.3,.4),100, replace= TRUE),
               id = seq(1,100,1))

# I am looking to reproduce this regression with jags
b=betareg(p ~ id, data= d, 
          link = c("logit"), link.phi = NULL, type = c("ML"))
summary(b)

Ниже я пытаюсь сделать ту же регрессию с помощью rjags.

#install.packages("rjags")
library(rjags)
jags_str = "
model {
#model

y ~ dbeta(alpha, beta)
alpha <- mu * phi
beta  <- (1-mu) * phi
logit(mu) <- a + b*id

#priors
a  ~ dnorm(0, .5)
b  ~ dnorm(0, .5)
t0 ~ dnorm(0, .5)
phi <- exp(t0)
}" 
id = d$id
y = d$p
model <- jags.model(textConnection(jags_str), 
                    data = list(y=y,id=id)
)
update(model, 10000, progress.bar="none"); # Burnin for 10000 samples
samp <- coda.samples(model, 
                     variable.names=c("mu"), 
                     n.iter=20000, progress.bar="none")

summary(samp)
plot(samp)

Я получаю ошибку в этой строке

 model <- jags.model(textConnection(jags_str), 
                        data = list(y=y,id=id)
    )

Error in jags.model(textConnection(jags_str), data = list(y = y, id = id)) : 
  RUNTIME ERROR:
Invalid vector argument to ilogit

Вы можете посоветовать

(1) как исправить ошибку

(2) как установить априорные значения для бета-регрессии

Спасибо.


bayesian jags rjags hierarchical-bayesian
person user3022875    schedule 14.08.2017    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Эта ошибка возникает из-за того, что вы передали вектор id скалярной функции logit. В Jags функции обратной ссылки не могут быть векторизованы. Чтобы решить эту проблему, вам нужно использовать цикл for для прохождения каждого элемента id. Для этого я, вероятно, добавил бы дополнительный элемент в ваш список данных, который обозначает длину id.

d = data.frame(p= sample(c(.1,.2,.3,.4),100, replace= TRUE),
           id = seq(1,100,1), len_id = length(seq(1,100,1)))

Оттуда вам просто нужно внести небольшое редактирование в код jags.

for(i in 1:(len_id)){
y[i] ~ dbeta(alpha[i], beta[i])
alpha[i] <- mu[i] * phi
beta[i]  <- (1-mu[i]) * phi
logit(mu[i]) <- a + b*id[i]
}

Однако, если вы отследите mu, это будет матрица 20000 (количество итераций) на 100 (длина id). Вас, вероятно, больше интересуют фактические параметры (a, b и phi).

person mfidino    schedule 15.08.2017