Байесовская полиномиальная регрессия с использованием пакета rjags

Я пытаюсь подогнать под модель полиномиальной логистической регрессии, используя rjags, поскольку результат представляет собой категориальную (номинальную) переменную (Результат) с 3 уровнями, а независимые переменные <сильные > Возраст (непрерывный) и Группа (трехуровневая категория). При этом я хотел бы получить апостериорные средние и 95% квантильных регионов для возраста и группы.

Я не очень хорошо разбираюсь в for loop, что, я думаю, является причиной того, что мой написанный код для модели не работает должным образом.

Мои априорные бета-версии подчиняются нормальному распределению, βj ∼ Normal (0,100) для j ∈ {0, 1, 2}.

Воспроизводимый код R

library(rjags)

set.seed(1)
data <- data.frame(Age = round(runif(119, min = 1, max = 18)),
                   Group = c(rep("pink", 20), rep("blue", 18), rep("yellow", 81)), 
                   Outcome = c(rep("A", 45), rep("B", 19), rep("C", 55)))

X <- as.matrix(data[,c("Age", "Group")]) 
J <- ncol(X)
N <- nrow(X)

## Step 1: Specify model
cat("
model {
for (i in 1:N){

    ##Sampling model
    yvec[i] ~ dmulti(p[i,1:J], 1)
    #yvec[i] ~ dcat(p[i, 1:J])  # alternative
    for (j in 1:J){
      log(q[i,j]) <- beta0 + beta1*X[i,1] + beta2*X[i,2] 
      p[i,j] <- q[i,j]/sum(q[i,1:J])  
    } 
    
    ##Priors
    beta0 ~ dnorm(0, 0.001)
    beta1 ~ dnorm(0, 0.001)
    beta2 ~ dnorm(0, 0.001)
}
}",
file="model.txt")

##Step 2: Specify data list 
dat.list <- list(yvec = data$Outcome, X=X, J=J, N=N) 

## Step 3: Compile and adapt model in JAGS 
jagsModel<-jags.model(file = "model.txt",
                      data = dat.list,
                      n.chains = 3,
                      n.adapt = 3000
)

Сообщение об ошибке:

введите описание изображения здесь

Источники, к которым я обращаюсь за помощью:

http://people.bu.edu/dietze/Bayes2018/Lesson21_GLM.pdf

Полиномиальная модель Дирихле в JAGS с категориальным X

Ссылка с сайта http://www.stats.ox.ac.uk/%7Enicholls/MScMCMC15/jags_user_manual.pdf, стр. 31

введите описание изображения здесь

Я только начал изучать, как использовать пакет rjags, поэтому я буду очень признателен за любые подсказки / объяснения и ссылки на соответствующие источники!


r logistic-regression bayesian rjags multinomial
person Minh Chau    schedule 17.08.2020    source источник
comment
Я думаю, вам просто нужно переместить априорные элементы за пределы for циклов. Так что просто переместите один из } из после приоров до них.   -  person MrFlick    schedule 17.08.2020
comment
Спасибо, хорошее место! Я попробовал ваше решение, но теперь получаю другое сообщение об ошибке. Невозможно разрешить следующие параметры: X [1,2] (строка 7) X [2,2] (строка 7) X [3,2] (строка 7) ....   -  person Minh Chau    schedule 17.08.2020

Ответы (1)


arrow_upward
3
arrow_downward

Я включу подход к вашей проблеме. Я взял те же априорные решения, которые вы определили для коэффициентов. Мне нужно только упомянуть, что, поскольку у вас есть коэффициент Group, я буду использовать один из его уровней в качестве эталона (в данном случае pink), поэтому его влияние будет учтено константой в модели. Далее код:

library(rjags)
#Data
set.seed(1)
data <- data.frame(Age = round(runif(119, min = 1, max = 18)),
                   Group = c(rep("pink", 20), rep("blue", 18), rep("yellow", 81)), 
                   Outcome = c(rep("A", 45), rep("B", 19), rep("C", 55)))

#Input Values we will avoid pink because it is used as reference level
#so constant absorbs the effect of that level
r1 <- as.numeric(data$Group=='pink')
r2 <- as.numeric(data$Group=='blue')
r3 <- as.numeric(data$Group=='yellow')
age <- data$Age
#Output 2 and 3
o1 <- as.numeric(data$Outcome=='A')
o2 <- as.numeric(data$Outcome=='B')
o3 <- as.numeric(data$Outcome=='C')
#Dim, all have the same length
N <- length(r2)

## Step 1: Specify model

model.string <- "
model{
for (i in 1:N){ 

## outcome levels B, C
o1[i] ~ dbern(pi1[i])
o2[i] ~ dbern(pi2[i]) 
o3[i] ~ dbern(pi3[i]) 

## predictors
logit(pi1[i]) <- b1+b2*age[i]+b3*r2[i]+b4*r3[i]
logit(pi2[i]) <- b1+b2*age[i]+b3*r2[i]+b4*r3[i]
logit(pi3[i]) <- b1+b2*age[i]+b3*r2[i]+b4*r3[i]

} 
## priors
b1 ~ dnorm(0, 0.001)
b2 ~ dnorm(0, 0.001)
b3 ~ dnorm(0, 0.001)
b4 ~ dnorm(0, 0.001)
}
"
#Model
model.spec<-textConnection(model.string)

## fit model w JAGS
jags <- jags.model(model.spec,
                   data = list('r2'=r2,'r3'=r3,
                               'o1'=o1,'o2'=o2,'o3'=o3,
                               'age'=age,'N'=N),
                   n.chains=3,
                   n.adapt=3000)

#Update the model
#Update
update(jags, n.iter=1000,progress.bar = 'none')
#Sampling
results <- coda.samples(jags,variable.names=c("b1","b2","b3","b4"),n.iter=1000,
                        progress.bar = 'none')
#Results
Res <- do.call(rbind.data.frame, results)

С результатами цепочек для параметров, сохраненными в Res, вы можете вычислить апостериорные медиа и достоверные интервалы, используя следующий код:

#Posterior means
apply(Res,2,mean)

         b1          b2          b3          b4 
-0.79447801  0.00168827  0.07240954  0.08650250

#Lower CI limit
apply(Res,2,quantile,prob=0.05)

         b1          b2          b3          b4 
-1.45918662 -0.03960765 -0.61027923 -0.42674155

#Upper CI limit
apply(Res,2,quantile,prob=0.95)

         b1          b2          b3          b4 
-0.13005617  0.04013478  0.72852243  0.61216838

Параметры b относятся к каждой из рассматриваемых переменных (age и уровни Group). Конечные значения могут измениться из-за смешанных цепочек!

person Duck    schedule 17.08.2020
comment
Большое спасибо за четкое объяснение каждой строчки кода и пояснение категориальной объясняющей переменной. Это очень полезно, и я понимаю все детали, поэтому еще раз спасибо! - person Minh Chau; 18.08.2020