Bayesianische multinomiale Regression unter Verwendung des rjags-Pakets

Question

Ich versuche, ein multinomiales logistisches Regressionsmodell anzupassen , das rjagsfür das Ergebnis eine kategoriale (nominale) Variable ( Ergebnis ) mit 3 Ebenen verwendet, und die erklärenden Variablen sind Alter (kontinuierlich) und Gruppe (kategorisch mit 3 Ebenen). Dabei möchte ich die Posterior-Mittelwerte und 95 %-Quantil-basierten Regionen für Alter und Gruppe erhalten .

Ich bin nicht wirklich großartig darin, for loopwas meiner Meinung nach der Grund ist, warum mein geschriebener Code für das Modell nicht richtig funktioniert.

Meine Beta-Prioritäten folgen einer Normalverteilung, βj ∼ Normal(0,100) für j ∈ {0, 1, 2}.

Reproduzierbarer R-Code

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
)

Fehlermeldung :

Quellen, die ich zur Hilfe gesucht habe :

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

Dirichlet-Multinomialmodell in JAGS mit kategorischem X

Referenz vonhttp://www.stats.ox.ac.uk/~nicholls/MScMCMC15/jags_user_manual.pdf, Seite 31

Ich habe gerade angefangen zu lernen, wie man das rjagsPaket benutzt, also wäre jeder Hinweis/Erklärung und Link zu relevanten Quellen sehr willkommen!

Answer 1

Ich werde eine Herangehensweise an Ihr Problem enthalten. Ich habe die gleichen Prioritäten genommen, die Sie für Koeffizienten definiert haben. Ich muss nur erwähnen, dass ich, da Sie einen Faktor in Grouphaben, eine seiner Ebenen als Referenz verwenden werde (in diesem Fall pink), sodass seine Wirkung von der Konstante im Modell berücksichtigt wird. Als nächstes der Code:

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)

Mit den Ergebnissen von Ketten für Parameter, die in gespeichert Ressind, können Sie spätere Medien und glaubwürdige Intervalle mit dem nächsten Code berechnen:

#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 

Die bParameter gehören zu jeder der betrachteten Variablen ( ageund den Ebenen von Group). Endwerte können sich aufgrund der gemischten Ketten ändern!