4
votes

I am attempting to call the following jags model in R:

model{
  # Main model level 1
  for (i in 1:N){
    ficon[i] ~ dnorm(mu[i], tau)
    mu[i] <- alpha[country[i]]
  }

  # Priors level 1 
  tau ~ dgamma(.1,.1)

  # Main model level 2
  for (j in 1:J){
    alpha[j] ~ dnorm(mu.alpha, tau.alpha)
  }

  # Priors level 2
  mu.alpha ~ dnorm(0,.01)
  tau.alpha ~ dgamma(.1,.1)

  sigma.1 <- 1/(tau)
  sigma.2 <- 1/(tau.alpha)

  ICC <- sigma.2 / (sigma.1+sigma.2)  
}

This is a hierarchical model, where ficon is a continuous variable 0-60, that may have a different mean or distribution by country. N = number of total observations (2244) and J = number of countries (34). When I run this model, I keep getting the following error message:

Compilation error on line 5.
Subset out of range: alpha[35]

This code worked earlier, but it's not working now. I assume the problem is that there are only 34 countries, and that's why it's getting stuck at i=35, but I'm not sure how to solve the problem. Any advice you have is welcome!

The R code that I use to call the model:

### input files JAGS ###
data <- list(ficon = X$ficon, country = X$country, J = 34, N = 2244)

inits1 <- list(alpha = rep(0, 34), mu.alpha = 0, tau = 1, tau.alpha = 1)
inits2 <- list(alpha = rep(1, 34), mu.alpha = 1, tau = .5, tau.alpha = .5)
inits <- list(inits1, inits2)

# call empty model 
eqlsempty <- jags(data, inits, model.file = "eqls_emptymodel.R",
                  parameters  = c("mu.alpha", "sigma.1", "sigma.2", "ICC"), 
                  n.chains = 2, n.iter = itt, n.burnin = bi, n.thin = 10)
1
There is a typographical error: on line 5 you attempt to index country with a value i ranging up to N=2244. That's obviously not right. - whuber
I disagree with @whuber. country comes from X$country which is a factor with N rows. The problem presumably is that within X$country some countries have been coded with numbers > 34 (in particular, it seems that there is a country coded as 35. - Jacob Socolar
@user2868853 Thank you; I think you're correct--and your diagnosis is consistent with joethorley's answer. - whuber

1 Answers

8
votes

To solve the problem you need to renumber your countries so they only have the values 1 to 34. If you only have 34 countries and yet you are getting the error message you state then one of the countries must have the value 35. To solve this one could call the following R code before bundling the data:

x$country <- factor(x$country)
x$country <- droplevels(x$country)
x$country <- as.integer(x$country)

Hope this helps