Search code examples
rbayesianjagsr2jags

How to fit a model with and without an interaction in a JAGS regression model


I'm using this tutorial to wrap my head around JAGS code. In the section 'Same model with an additional categorical predictor' it states that "This model includes an interaction between sex and body length". How can I remove this so that there's no interaction?

Here's the full setup and model in R and JAGS.

First the data:

set.seed(42)

samplesize <- 50 # Larger sample size because we're fitting a more complex model
b_length <- sort(rnorm(samplesize)) # Body length
sex <- sample(c(0, 1), size = samplesize, replace = T) # Sex (0: female, 1: male)

int_true_f <- 30 # Intercept of females
int_true_m_diff <- 5 # Difference between intercepts of males and females
slope_true_f <- 10 # Slope of females
slope_true_m_diff <- -3 # Difference between slopes of males and females

mu <- int_true_f + sex * int_true_m_diff + (slope_true_f + sex * slope_true_m_diff) * b_length # True means
sigma <- 5 # True standard deviation of normal distributions

b_mass <- rnorm(samplesize, mean = mu, sd = sigma) # Body mass (response variable)

# Combine into a data frame:
snakes2 <- data.frame(b_length = b_length, b_mass = b_mass, sex = sex)
head(snakes2)

jagsdata_s2 <- with(snakes2, list(b_mass = b_mass, b_length = b_length, sex = sex, N = length(b_mass)))

JAGS code:

lm2_jags <- function(){
    # Likelihood:
    for (i in 1:N){
        b_mass[i] ~ dnorm(mu[i], tau) # tau is precision (1 / variance)
        mu[i] <- alpha[1] + sex[i] * alpha[2] + (beta[1] + beta[2] * sex[i]) * b_length[i]
    }
    # Priors:
    for (i in 1:2){
        alpha[i] ~ dnorm(0, 0.01)
        beta[i] ~ dnorm(0, 0.01)
    }
    sigma ~ dunif(0, 100)
    tau <- 1 / (sigma * sigma)
}

Initial values and run:

init_values <- function(){
    list(alpha = rnorm(2), beta = rnorm(2), sigma = runif(1))
}

params <- c("alpha", "beta", "sigma")

fit_lm2 <- jags(data = jagsdata_s2, inits = init_values, parameters.to.save = params, model.file = lm2_jags,
               n.chains = 3, n.iter = 12000, n.burnin = 2000, n.thin = 10, DIC = F)

Solution

  • The interaction term is contained in your calculation of mu. The sex changes how the formula between body length and body mass is defined, via the slope terms. To build a model where sex and body length are treated as independent with respect to how they affect body mass, you could do something like this:

    mu <- int_true_f + (sex * int_true_m_diff) + b_length 
    

    The JAGS code would then become

    lm2_jags <- function(){
      # Likelihood:
      for (i in 1:N){
        b_mass[i] ~ dnorm(mu[i], tau) # tau is precision (1 / variance)
        mu[i] <- alpha[1] + (sex[i] * alpha[2]) + (b_length[i] * alpha[3])
      }
      # Priors:
      for (i in 1:3){
        alpha[i] ~ dnorm(0, 0.01)
      }
      sigma ~ dunif(0, 100)
      tau <- 1 / (sigma * sigma)
    }