Régression multinomiale bayésienne utilisant le package rjags

Aug 17 2020

J'essaie d'adapter un modèle de régression logistique multinomial en utilisant rjagspour le résultat une variable catégorielle (nominale) ( Résultat ) à 3 niveaux, et les variables explicatives sont l' âge (continu) et le groupe (catégoriel à 3 niveaux). Ce faisant, j'aimerais obtenir les moyennes postérieures et les régions basées sur les quantiles à 95 % pour l'âge et le groupe .

Je ne suis pas vraiment doué for loopet je pense que c'est la raison pour laquelle mon code écrit pour le modèle ne fonctionne pas correctement.

Mes a priori bêta suivent une distribution normale, βj ∼ Normal(0,100) pour j ∈ {0, 1, 2}.

Code R reproductible

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
)

Message d'erreur :

Sources que j'ai consultées pour obtenir de l'aide :

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

Modèle multinomial de Dirichlet dans JAGS avec X catégoriel

Référence dehttp://www.stats.ox.ac.uk/~nicholls/MScMCMC15/jags_user_manual.pdf, page 31

Je viens juste de commencer à apprendre à utiliser le rjagspackage, donc tout conseil/explication et lien vers des sources pertinentes seraient grandement appréciés !

Réponses

3 Duck Aug 17 2020 at 23:20

Je vais inclure une approche à votre problème. J'ai pris les mêmes priors que vous avez définis pour les coefficients. J'ai seulement besoin de mentionner que, comme vous avez un facteur Group, j'utiliserai l'un de ses niveaux comme référence (dans ce cas pink) afin que son effet soit pris en compte par la constante du modèle. Ensuite le 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)

Avec les résultats des chaînes pour les paramètres enregistrés dans Res, vous pouvez calculer les médias postérieurs et les intervalles crédibles à l'aide du code suivant :

#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 

Les bparamètres appartiennent à chacune des variables considérées ( ageet les niveaux de Group). Les valeurs finales pourraient changer à cause des chaînes mixtes !