Interaction terms

Modeling
1

Hi,

  1. Can you advise how to add an interaction term between the categorical variable “groups2” which has 3 categories (N_GROUPS2=3) and the continuous variable “X”? In this case there should be 3 beta2 coefficients. There is some incomplete code below for this.

  2. To allow for variation between group1 categories, can you advise how allow each category in group1 to have their own interaction terms for the “group2” by “X” interaction. N_GROUP1=4 so in this case there should be 12 coefficients. To reduce code I am keeping assuming all group coefficients come from a fixed prior distribution (not centered or non-centered).

Thank you

data{
  int<lower=1> N;
  int<lower=1> N_GROUPS1;
  int<lower=1> N_GROUPS2;
  int<lower=1, upper=N_GROUPS1> groups1[N];
  int<lower=1, upper=N_GROUPS2> groups2[N];
  real X[N];
  real Y[N];
}

parameters{
  vector[N_GROUPS1] alpha;
  vector[N_GROUPS1] beta1;
  vector[N_GROUPS2] beta2;
  real sigma;
}


model{
  alpha~normal(0,5);
  beta1~normal(0,5); 
  beta2~normal(0,5); 
  sigma~normal(0,5);
  for(i in 1:N){
    Y[i]~normal(alpha[groups1[i]]+beta1[groups1[i]]*X[i]  + beta2[groups2[i]]* ??,sigma);
  }
}
2

  1. The usual approach here is varying slopes. So rather than an interaction predictor, you use a different “slope” parameter for each of the groups. That means beta1[groups1[i]], so the slope already varies by group 1. To add variation by group 2, just follow the same plan as for (2) below, but with the coefficient.

  2. Sticking to your notation you introduce an interaction like this:

matrix[N_GROUPS1, N_GROUPS2] beta12;

Then the extra term is just beta12[groups1[i], groups2[i]].

You might want to consider hierarchical priors for these.

P.S. The loop for the likelihood can be vectorized for much more efficient code:

 for(i in 1:N){
    Y[i] ~ normal(alpha[groups1[I]]
                  + beta1[groups1[i]] * X[i]
                  + beta2[groups2[I]]
                  + beta12[groups1[i], groups2[i]],
                  sigma);
  }

to

Y ~ normal(alpha[groups1] + beta1[groups1] .* X + beta2[groups2], sigma)