jeudi 2 août 2018

Nested random effect in cumulative proportional odds model with JAGS

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