Regressione multinomiale bayesiana utilizzando il pacchetto rjags

Question

Sto cercando di adattare un modello di regressione logistica multinomiale utilizzando rjagsper il risultato una variabile categorica (nominale) ( Risultato ) con 3 livelli e le variabili esplicative sono Età (continua) e Gruppo (categorica con 3 livelli). In tal modo, vorrei ottenere le medie posteriori e le regioni basate sui quantili al 95% per età e gruppo .

Non sono molto bravo, il for loopche penso sia il motivo per cui il mio codice scritto per il modello non funziona correttamente.

I miei beta priori seguono una distribuzione Normale, βj ∼ Normal(0,100) per j ∈ {0, 1, 2}.

Codice R riproducibile

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
)

Messaggio di errore :

Fonti che ho cercato aiuto :

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

Modello multinomiale di Dirichlet in JAGS con X categorico

Riferimento dahttp://www.stats.ox.ac.uk/~nicholls/MScMCMC15/jags_user_manual.pdf, pagina 31

Ho appena iniziato a imparare come utilizzare il rjagspacchetto, quindi qualsiasi suggerimento/spiegazione e collegamento a fonti pertinenti sarebbe molto apprezzato!

Answer 1

Includerò un approccio al tuo problema. Ho preso gli stessi precedenti che hai definito per i coefficienti. Devo solo menzionare che, poiché hai un fattore Group, userò uno dei suoi livelli come riferimento (in questo caso pink), quindi il suo effetto sarà preso in considerazione dalla costante nel modello. Avanti il ​​codice:

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)

Con i risultati delle catene per i parametri salvati in Res, puoi calcolare i media posteriori e gli intervalli credibili utilizzando il codice successivo:

#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 

I bparametri appartengono a ciascuna delle variabili considerate ( ageei livelli di Group). I valori finali potrebbero cambiare a causa delle catene miste!