rjags 패키지를 사용한 베이지안 다항 회귀
결과 를 위해 사용 하는 다항 로지스틱 회귀 모델 을 맞추려고 노력하고 있습니다 .3 수준 rjags의 범주 형 (명목) 변수 ( 결과 )이고 설명 변수는 연령 (연속) 및 그룹 (3 단계 범주)입니다. 이를 통해 Age 및 Group에 대한 사후 평균 및 95 % 분위수 기반 영역을 얻고 싶습니다 .
나는 for loop모델에 대해 작성한 코드가 제대로 작동하지 않는 이유라고 생각하는 것이 정말 훌륭 하지 않습니다.
내 베타 사전은 j ∈ {0, 1, 2}에 대한 정규 분포 βj ∼ Normal (0,100)을 따릅니다.
재현 가능한 R 코드
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
)
오류 메시지 :
내가 도움을 구한 출처 :
http://people.bu.edu/dietze/Bayes2018/Lesson21_GLM.pdf
범주 형 X를 사용하는 JAGS의 Dirichlet 다항 모델
참조 에서http://www.stats.ox.ac.uk/~nicholls/MScMCMC15/jags_user_manual.pdf, 페이지 31
rjags패키지 사용 방법을 배우기 시작 했으므로 힌트 / 설명 및 관련 소스에 대한 링크는 크게 감사하겠습니다!
답변
귀하의 문제에 대한 접근 방식을 포함하겠습니다. 나는 계수에 대해 정의한 것과 동일한 사전을 취했습니다. 요소 Group가 있으므로 해당 수준 중 하나를 참조 (이 경우 pink)로 사용하므로 그 효과가 모델의 상수에 의해 고려 된다는 점만 언급하면됩니다 . 다음 코드 :
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)
에 저장된 매개 변수에 대한 체인 결과로 Res다음 코드를 사용하여 사후 미디어 및 신뢰할 수있는 간격을 계산할 수 있습니다.
#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 b parameters belong to the each of the variables considered (age and the levels of Group). Final values could change because of the mixed chains!