Search code examples
rlinear-regressionjagsrjags

Jags: Attempt to redefine node error, mixed effect regression

I want to perform a mixed effect regression in rjags, with a random slope and intercept. I define the following toy dataset:

library(ggplot2)
library(data.table)

global_slope <- 1
global_int <- 1

Npoints_per_group <- 50
N_groups <- 10
pentes <- rnorm(N_groups,-1,.5)

centers_x <- seq(0,10,length = N_groups)
center_y <- global_slope*centers_x + global_int

group_spread <- 2
group_names <- sample(LETTERS,N_groups)

df <- lapply(1:N_groups,function(i){
  x <- seq(centers_x[i]-group_spread/2,centers_x[i]+group_spread/2,length = Npoints_per_group)
  y <- pentes[i]*(x- centers_x[i])+center_y[i]+rnorm(Npoints_per_group)
  data.table(x = x,y = y,ID = group_names[i])
}) %>% rbindlist()

ggplot(df,aes(x,y,color = as.factor(ID)))+
  geom_point()

enter image description here

This is a typical situation of Simpson paradox: an overall increasing trend when you have a decreasing trend within each group (given by the ID variable).

I define the following model:

library(rjags)

model_code_simpson <- 
" model
{

 # first level
  for (i in 1:n) {
    y[i] ~ dnorm(alpha[i] + beta[i] * x[i], tau)
    alpha[i] = alpha[group[i]] # random intercept
    beta[i] = beta[group[i]] # random slope
  }

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

# Priors
mu.alpha ~ dnorm(0,0.001)
mu.beta ~ dnorm(0,0.001)
sigma ~ dunif(0,10)
sigma.alpha ~ dunif(0,10)
sigma.beta ~ dunif(0,10)

# Derived quantities
tau <- pow(sigma,-2)
tau.alpha <- pow(sigma.alpha,-2)
tau.beta <- pow(sigma.beta,-2)

}
"

# Choose the parameters to watch
model_parameters <- c("mu.alpha","tau.alpha","tau.beta","tau")

# define numeric grouping variable
df[,ID2 := .GRP,by = ID]

model_data <- list(n = nrow(df), 
                   y = df$y, 
                   x = df$x,
                   group = df$ID2,
                   J = df[,uniqueN(ID)])

model <- jags.model(textConnection(model_code_simpson),
                    data = model_data,
                    n.chains = 2)

I get the following error:

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Deleting model

Error in jags.model(textConnection(model_code_simpson), data = model_data,  : 
  RUNTIME ERROR:
Compilation error on line 8.
Attempt to redefine node beta[1]

I do not understand what is happening, and related questions did not help me much.

Solution

  • You defined beta twice. First, beta is a vector of length n when you are looping through the data. Second, beta is a vector of length J when you are creating the random effects. This "redefining" is causing this issue, but it is an easy fix. You just need to remove that first instance of beta in your model and it will compile (i.e., just move your nested indexing inside of dnorm() and you are good to go).

    model_code_simpson <- 
  " model
{

 # first level
  for (i in 1:n) {
    y[i] ~ dnorm(
      alpha[group[i]] + beta[group[i]] * x[i],
      tau
    )
  }

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

# Priors
mu.alpha ~ dnorm(0,0.001)
mu.beta ~ dnorm(0,0.001)
sigma ~ dunif(0,10)
sigma.alpha ~ dunif(0,10)
sigma.beta ~ dunif(0,10)

# Derived quantities
tau <- pow(sigma,-2)
tau.alpha <- pow(sigma.alpha,-2)
tau.beta <- pow(sigma.beta,-2)

}
"