Random intercepts in hierarchical Dirichlet regression (jags)

Question

I have the following data structure:

• y: 3 columns that are observed proportions of deaths over the years.
• x1: GDP - continuous variable related to each year
• x2: Ages- related to deaths

Here the model specification:

Model

Here simulated data:

library(tidyverse)

Year <-2000:2010 
n.year <- length(Year)
Ages <- rep(seq(0,50,5),n.year)

# X2 AGES
n.ages <- length(Ages)

# Y Causes of Death 
Cause1 <- rpois(n.ages,60)
Cause2 <- rpois(n.ages,20)
Cause3 <- rpois(n.ages,90)

# Get proportion
C <- cbind(Cause1,Cause2,Cause3)
C_Porp <- C/rowSums(C)


Year.tot <- rep(Year,length(unique(Ages)))

DATA <- data.frame(Year=Year.tot,Ages=sort(Ages),
                   C1=C_Porp[,1],
                   C2=C_Porp[,2],
                   C3=C_Porp[,3]) 

# X1 GDP OVER YEARS
GDP <- rnorm(n.year,50,15)
GDP <- as.data.frame(cbind(GDP,Year))

# TOT
D <- as.matrix(merge(DATA,GDP))[,-1]
head(D)

Jags Model

library(runjags)
dirlichet.model = 
  
  
  "model {
  #setup priors for each species
  for(j in 1:N.c){
    m0[j] ~ dnorm(0, 1.0E-3) #intercept prior
    m1[j] ~ dnorm(0, 1.0E-3) #      mat prior
  }
  
  #implement dirlichet
  for(i in 1:N){
  
    y[i,1:N.c] ~ ddirch(a0[i,1:N.c]) 
     
    for(j in 1:N.c){
    
      a0[i,j] ~ dnorm(mu[i,j],0.001) 
      
      log(mu[i,j]) <- (m0[j]+alpha[Ages[i]] + m1[j] * GDP[i])
    
    }
    
     # Effect of site:
    for (s in 1:S){
      
      alpha[s] ~ dnorm(0,0.001)
      
    }}
    
  }
"

jags.data <- list(y = D[,c(3,4,5)],mat = D[,5], Ages=D[,1],N = nrow(D[,c(3,4,5)]), 
                  N.c = ncol(D[,c(3,4,5)]),S=(length(unique(D[,1]))))

jags.out <- run.jags(dirlichet.model,
                     data=jags.data,
                     adapt = 500,
                     burnin = 5000,
                     sample = 10000,
                     n.chains=5,
                     monitor=c('m0','m1',"mu"))
out <- summary(jags.out)

I got this error:

Errore: The following error occured when compiling and adapting the model using rjags:
 Error in rjags::jags.model(model, data = dataenv, n.chains = length(runjags.object$end.state),  : 
  RUNTIME ERROR:
Compilation error on line 12.
Index out of range taking subset of  alpha

Any suggestion for model specification?

Answer 1

There were several other small errors, too. This code should work now, though. The problems were:

1. Age was defined as a sequence from 0 to 50 by 5, but then you use Age to index alpha. I think what you really want is one different alpha for each value of Age. I solved this by doing the following in the data: Ages=match(D[,1], unique(D[,1])),
2. alpha was being re-defined because the loop over N enclosed both the loops over N.c and S. By closing the loop over N before the loop over S starts, that solves the problem.
3. GDP was defined in the data as mat, so I replaced mat = D[,5] with GDP = D[,5].

After those changes, the model runs.

library(tidyverse)

Year <-2000:2010 
n.year <- length(Year)
Ages <- rep(seq(0,50,5),n.year)

# X2 AGES
n.ages <- length(Ages)

# Y Causes of Death 
Cause1 <- rpois(n.ages,60)
Cause2 <- rpois(n.ages,20)
Cause3 <- rpois(n.ages,90)

# Get proportion
C <- cbind(Cause1,Cause2,Cause3)
C_Porp <- C/rowSums(C)


Year.tot <- rep(Year,length(unique(Ages)))

DATA <- data.frame(Year=Year.tot,Ages=sort(Ages),
                   C1=C_Porp[,1],
                   C2=C_Porp[,2],
                   C3=C_Porp[,3]) 

# X1 GDP OVER YEARS
GDP <- rnorm(n.year,50,15)
GDP <- as.data.frame(cbind(GDP,Year))

# TOT
D <- as.matrix(merge(DATA,GDP))[,-1]
head(D)



library(runjags)
dirlichet.model = 
  
  
  "model {
  #setup priors for each species
  for(j in 1:N.c){
    m0[j] ~ dnorm(0, 1.0E-3) #intercept prior
    m1[j] ~ dnorm(0, 1.0E-3) #      mat prior
  }
  
  #implement dirlichet
  for(i in 1:N){
  
    y[i,1:N.c] ~ ddirch(a0[i,1:N.c]) 
     
    for(j in 1:N.c){
    
      a0[i,j] ~ dnorm(mu[i,j],0.001) 
      
      log(mu[i,j]) <- (m0[j]+alpha[Ages[i]] + m1[j] * GDP[i])
    
    }
  }
     # Effect of site:
    for (s in 1:S){
      
      alpha[s] ~ dnorm(0,0.001)
      
    }
    
  }
"

jags.data <- list(y = D[,c(2,3,4)],
                  GDP = D[,5], 
                  Ages=match(D[,1], unique(D[,1])),
                  N = nrow(D[,c(3,4,5)]), 
                  N.c = ncol(D[,c(3,4,5)]),
                  S=(length(unique(D[,1]))))

jags.out <- run.jags(dirlichet.model,
                     data=jags.data,
                     adapt = 500,
                     burnin = 5000,
                     sample = 10000,
                     n.chains=5,
                     monitor=c('m0','m1',"mu"))
out <- summary(jags.out)