Regresi Multinomial Bayesian menggunakan paket rjags

Aug 17 2020

Saya mencoba menyesuaikan model regresi logistik multinomial menggunakan rjagsuntuk hasil adalah variabel kategorikal (nominal) ( Hasil ) dengan 3 level, dan variabel penjelasnya adalah Usia (kontinu) dan Kelompok (kategorikal dengan 3 level). Dengan melakukan itu, saya ingin mendapatkan rata-rata Posterior dan 95% wilayah berbasis kuantil untuk Usia dan Grup .

Saya tidak terlalu hebat for loopyang menurut saya merupakan alasan mengapa kode tertulis saya untuk model tidak berfungsi dengan baik.

Priors beta saya mengikuti distribusi Normal, βj ∼ Normal (0,100) untuk j ∈ {0, 1, 2}.

Kode R yang dapat direproduksi

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
)

Pesan kesalahan :

Sumber yang saya cari bantuannya :

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

Model Dirichlet Multinomial di JAGS dengan kategori X

Referensi darihttp://www.stats.ox.ac.uk/~nicholls/MScMCMC15/jags_user_manual.pdf, halaman 31

Saya baru saja mulai mempelajari cara menggunakan rjagspaket ini sehingga petunjuk / penjelasan apa pun dan tautan ke sumber yang relevan akan sangat kami hargai!

Jawaban

3 Duck Aug 17 2020 at 23:20

Saya akan menyertakan pendekatan untuk masalah Anda. Saya telah mengambil prioritas yang sama dengan yang Anda tetapkan untuk koefisien. Saya hanya perlu menyebutkan bahwa karena Anda memiliki faktor, Groupsaya akan menggunakan salah satu levelnya sebagai referensi (dalam hal ini pink) sehingga pengaruhnya akan diperhitungkan oleh konstanta dalam model. Selanjutnya kodenya:

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)

Dengan hasil rantai untuk parameter yang disimpan Res, Anda dapat menghitung media posterior dan interval yang kredibel menggunakan kode berikut:

#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 

The bparameter milik masing-masing variabel dianggap ( agedan tingkat Group). Nilai akhir bisa berubah karena rantai campuran!