Help with Varying Intercept, Varying Slope Model

Modeling
1

Hi everyone,

I’m trying to learn how to program multilevel models in Stan and I’m having difficulty programing a varying intercept and varying slope model with a group-level predictor. I’m using the Radon data and following the Gelman and Hill book. I’m trying to replicate the results of this model in Stan:

lmer(formula = log_radon ~ floor + log_uranium + (1 + floor | county) - 1, data = radon))

where log_uranium is the group-level predictor.

So far my Stan program looks like this:

data {
  int n; // number of observations
  int n_county; // number of counties
  int<lower = 0, upper = n_county> county[n]; 
  
  // level 1 variables
  vector[n] log_radon;
  vector[n] vfloor;
  
  // level 2 variables
  vector[n] log_uranium;
}

parameters {
  vector[n_county] a; // vector of county intercepts
  vector[n_county] b_floor; // floor slope parameter
  real b_uranium; // uranium slope parameter
  real mu_a; // mean of counties
  real<lower = 0> sigma_a; // variance of counties
  real<lower = 0> sigma_y; // model residual variance
}

transformed parameters {
  vector[n] mu;
  vector[n] alpha;
  
  // level 2 equation
  alpha = a[county] + b_uranium * log_uranium;
  
  // level 1 equation
  mu = alpha + b_floor[county] .* vfloor;
}

model {
  // priors
  mu_a ~ normal(0, 100);
  b_floor ~ normal(0, 100);
  b_uranium ~ normal(0, 100);
  sigma_a ~ uniform(0, 10);
  sigma_y ~ uniform(0, 10);
  
  // level 2 likelihood
  a ~ normal(mu_a, sigma_a);
  
  // level 1 likelihood
  log_radon ~ normal(mu, sigma_y);
}

This model runs, but the random effects estimates don’t match the results of the lmer function. I suspect that this is because I am not estimating the covariance matrix in the Stan model, but I’m not sure how to add that in.

Do you think that is the problem with the model, and if so, how can I add the covariance matrix to this code? Any help would be greatly appreciated!

1 Like
2

@bgoodri and/or @jonah might be able to help.

3

You are going to need a bivariate normal prior on a and b_floor. It is probably easiest to start with

brms::make_stancode(log_radon ~ floor + log_uranium + (1 + floor | county) - 1, data = radon)
1 Like