Rjags paketini kullanan Bayes Çok Terimli Regresyon
Bir sığdırıyorum multinomial lojistik regresyon modeli kullanılarak rjags
sonuç için bir kategorik (nominal) değişken (olduğu Sonuç 3 düzeyleri ile) ve açıklayıcı değişkenler Yaş (sürekli) ve Grup (3 seviyeleri ile kategorik). Bunu yaparken, Yaş ve Grup için Posterior ortalamaları ve% 95 kuantil bazlı bölgeleri elde etmek istiyorum .
Ben de gerçekten büyük değilim for loop
maket için benim yazılı kodu düzgün çalışmadığını gösterir nedeni olduğunu düşünüyorum hangi.
Beta önceliklerim j ∈ {0, 1, 2} için Normal dağılım, βj ∼ Normal (0,100) izliyor.
Tekrarlanabilir R kodu
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
)
Hata mesajı :

Yardım için aradığım kaynaklar :
http://people.bu.edu/dietze/Bayes2018/Lesson21_GLM.pdf
Kategorik X ile JAGS'de Dirichlet Multinomial model
Referans danhttp://www.stats.ox.ac.uk/~nicholls/MScMCMC15/jags_user_manual.pdf, sayfa 31

rjags
Paketi nasıl kullanacağımı yeni öğrenmeye başladım, bu nedenle herhangi bir ipucu / açıklama ve ilgili kaynaklara bağlantı büyük ölçüde takdir edilecektir!
Yanıtlar
Sorununuza bir yaklaşım ekleyeceğim. Katsayılar için tanımladığınız aynı önceleri aldım. Sadece bir faktörünüz Group
olduğu için onun seviyelerinden birini referans olarak kullanacağım (bu durumda pink
), böylece etkisinin modeldeki sabit tarafından dikkate alınacağını belirtmek istiyorum. Ardından kod:
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)
Kaydedilen parametreler için zincirlerin sonuçlarıyla, Res
sonraki kodu kullanarak arka ortamı ve güvenilir aralıkları hesaplayabilirsiniz:
#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
Parametreler olarak değişken (her birine ait age
ve seviyeleri Group
). Karma zincirler nedeniyle nihai değerler değişebilir!