I have to fit a cumulative proportional odds model with JAGS. I want to explain an ordinal variable (deb, on the scale from 0 to 5) with a categorical variable (PRO, 25 levels) and taking into account the variability of the family (DESC, 362 levels) as a random effect which is nested in PRO. Each PRO doesn't have the same number of DESC and each DESC doesn't have the same number of individuals. I want to see if PRO affects deb and I want to estimate the family variability.
The frequentist model would be , with {lme4} :
mod <- lmer(deb ~ PRO + (1|PRO/DESC), data=data)
I wrote this code, with the help of this discussion : Nested Random effect in JAGS/ WinBUGS
prov <- as.numeric(data$PRO)
famille <- as.numeric(data$DESC)
model <- "model{
for(i in 1:N){
mu[i] <- betapro[prov[i]]
logit(Q[i,1]) <- alpha[1]-mu[i]
p[i,1] <- Q[i,1]
for(j in 2:4){
logit(Q[i,j]) <- alpha[j]-mu[i]
p[i,j] <- Q[i,j] - Q[i,j-1]
}
p[i,5] <- 1 - Q[i,4]
y[i] ~ dcat(p[i,])
}
## priors over thresholds
for(j in 1:4){
alpha0[j] ~ dnorm(0,1.0E-3)
}
alpha <- sort(alpha0)
for(i in 1:N){
betapro[prov[i]] ~ dnorm(0, taupro)
betadesc[famille[i]] ~ dnorm(betapro[prov[i]], taudesc)
}
taupro ~dgamma(0.01, 0.01)
taudesc ~dgamma(0.01, 0.01)
}
"
writeLines(model,con="model_PRO_DESC.txt")
datalist <- list(N=nrow(data), y=data$deb,
prov=prov, famille=famille)
initiale <- list(alpha0 = c(-0.8, -0.3, 0.2, 0.7), betapro=rep(0,25), betadesc=rep(0, 362))
nchain <- 3
init <- list(list(alpha0=initiale$alpha0-2,
betapro=initiale$betapro-2,
betadesc=initiale$betadesc-2),
initiale,
list(alpha0=initiale$alpha0+2,
betapro=initiale$betapro+2,
betadesc=initiale$betadesc+2))
Nechantillons <- 10000
burnin <- 1000
eclair<-1
parameters <- c("alpha", "betapro")
DEB13_PRO_DESC <-jags(model.file="model_PRO_DESC.txt", data=datalist, n.chain = nchain, inits=init, parameters.to.save=parameters,
n.iter=Nechantillons, n.burnin=burnin, n.thin=eclair)
And I got this error message but I don't understand where is the problem.
Error in jags.model(model.file, data = data, inits = init.values, n.chains = n.chains, :
RUNTIME ERROR:
Compilation error on line 25.
Attempt to redefine node betapro[11]
Is anybody can help me ?
I hope hearing from you soon !
Aucun commentaire:
Enregistrer un commentaire