Multilevel model with two predictors and their interaction

Modeling
,
1

Hey folks!

I’m a Stan novice, so I wanted to check whether I had coded a multilevel model with an interaction term at level 1 correctly. I’m using cmdstanr. If any advice or resources to refer to would be very very welcome!

I’m running an experiment where I’ve recruited participants from three sources. The basic elements of the experiment are participants are randomly assigned to one of two conditions, complete a self-report measure of sexual attractiveness, and respond about how likely they are to engage in a certain behavior. I want to test the effects of experimental condition, sexual attractiveness, and their interaction on their behavioral intentions. Because of the three recruitment sources, I want to run a multilevel model (assuming some dependency in these variables attributable to the recruitment source).

I’ve prepared the data as a matrix (using this handy tutorial) where columns represent the intercept, binary experimental condition (either 0 or 1), a continuous self-report measure, and their interaction (so condition*the measure), and rows represent individual participants.

The model I’ve put together is this:

data {
  int n; //number of observations
  int s; //number of samples
  int k; // number of predictors + intercept, so: 1 intercept and 3 slopes 
  array[n] int id; //participant id
  matrix[n, k] x; // matrix of predictors 
  array[n] real y; //outcome vector
} 

parameters {
matrix[k, s] zp; //matrix of intercepts and slopes
vector<lower = 0>[k] sigma_s; //variance of slopes and intercepts
vector[k] beta; //multivariate normal priors for sample intercepts and slopes
cholesky_factor_corr[k] lp; //correlation matrix of slopes and intercepts
real<lower = 0> sigma; //population sigma
}

transformed parameters{
  matrix[k, s] z; //non-centered sample parameters. column 1 = sample intercepts, column 2-4 = sample slopes
  z = diag_pre_multiply(sigma_s, lp)*zp;
}

model {
vector[n] mu; // mean

// priors
beta ~ normal(0, 1); //priors for intercepts and slopes
lp ~ lkj_corr_cholesky(2); //prior for their correlation
sigma_s ~ exponential(1); //sample variance
sigma ~ exponential(1); // population variance
to_vector(zp) ~ normal(0, 1); //sample parameters are drawn from a population distribution of intercepts and slopes


// likelihood 

for(i in 1:n) {
  mu[i] = beta[1] + z[1, id[i]] + //intercept
  (beta[2] + z[2, id[i]])*x[i,2] + //effect of condition
  (beta[3] + z[3, id[i]])*x[i,3] + //effect of desirability
  (beta[4] + z[4, id[i]])*x[i,4]; //their interaction
}

y ~ normal(mu, sigma);

}

generated quantities {
  matrix[k,k] omega; //the correlation matrix for intercepts and slopes
  omega = multiply_lower_tri_self_transpose(lp);
 }

I’ve run the model on simulated data, and it seems to mix fine. My biggest concern is how I’ve model the interaction, which is the substantive relationship of interest in this study.

Thanks heaps in advance!

2

Sorry this didn’t get answered sooner—this is a relatively easy one, because your model looks good as is.

  1. I would suggest including bounds on your data for error checking, like int<lower=0> n.

  2. You can use compound declare define:

matrix[k, s] z = diag_pre_multiply(sigma_s, lp) * zp;

I would strongly recommend space around operators for readability.

We really need to find a way to make this whole operation easier.

  1. You can simplify your mu calculation to
vector[n] mu = beta[1] + z[1, id]
             + (beta[2] + z[2, id]) .* x[ , 2] 
             + (beta[3] + z[3, id]) .* x[ , 3]
             + (beta[4] + z[4, id]) .* x[ , 4];

Then you can make it more efficient by transposing x and defining it as an array of vectors (that eliminates all the copies otherwise involved in pulling out a column like x[ , 2]. Same for z.

Otherwise, this all looks fine.

3

Thanks Bob! I’ve put in your suggestions. :) Definitely agree the spaces between operators saves some eye strain.

I had another question and sorry if this is a super basic follow-up. Would you be able to recommend how best to transpose x to an array of vectors? I’ve looked up the Stan documentation for transposing matrices, but can’t quite wrap my head around how to do that in my code.

4

Sure. Your original was

matrix[N, K] x;

I converted your n and k to caps to make indexing easier below. Here’s a transposed version implemented as an array of vectors.

  array[K] vector[N] x;

You can also make these row_vector types if that’s more convenient. The benefit is that x[k] here does not need to do any copying—it just returns the k-th array element. On the other hand, when you have the original definition and look at x[ , k], it has to allocate a vector of size K, then walk over the k-th column of x and copy each element into the new vector. It’s not so much the speed, which is fast, but the memory pressure you want to avoid, since Stan programs are typically more memory than CPU bound.

1 Like