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!