Following the general advice of starting with something simple before building the full model, I wanted to verify that I understood how to write a simple linear regression with a varying intercept. Here’s the Stan code:

data { int n_species; int n_sites; int n_obs; int species[n_obs]; int site[n_obs]; real intercept_scale; real beta_scale; real sigma_scale; vector[n_obs] x; vector[n_obs] y; } parameters { real beta_0; real beta_1; real sigma; real sigma_species; real sigma_species_site; vector[n_species] beta_species; matrix[n_species,n_sites] beta_species_site; } transformed parameters { vector[n_obs] mu_obs; vector[n_species] mu_species; matrix[n_species,n_sites] mu_species_site; for (j in 1:n_species) { mu_species[j] = beta_0 + beta_species[j]; for (k in 1:n_sites) { mu_species_site[j,k] = mu_species[j] + beta_species_site[j,k]; } } for (i in 1:n_obs) { mu_obs[i] = mu_species_site[species[i], site[i]] + beta_1*x[i]; } } model { for (i in 1:n_obs) { y[i] ~ normal(mu_obs[i], sigma); } beta_0 ~ normal(0.0, intercept_scale); for (j in 1:n_species) { beta_species[j] ~ normal(0.0, sigma_species); for (k in 1:n_sites) { beta_species_site[j,k] ~ normal(0.0, sigma_species_site); } } beta_1 ~ normal(0.0, beta_scale); sigma ~ cauchy(0.0, sigma_scale); sigma_species ~ gamma(1.0, 1.0); //cauchy(0.0, sigma_scale); sigma_species_site ~ gamma(1.0, 1.0); //cauchy(0.0, sigma_scale); }

To check my understanding, I compared the results I obtained from running this code with those I obtained from running stan_lmer() on the same data:

fit.lmer <- stan_lmer(y ~ x + (1|species/site), data=dat, prior=normal(0.0, beta_scale), prior_intercept=normal(0.0, intercept_scale), prior_aux=cauchy(0.0, sigma_scale), prior_covariance=decov(regularization=regularize, concentration=1.0, shape=1.0, scale=1.0), adapt_delta=0.99)

Estimates for the intercept, regression coefficient (slope), and residual standard deviation match pretty closely. The estimates for the random effect standard deviations and the corresponding random effects differ quite a bit.

In the simulations I’ve done (code in Github at https://kholsinger.github.io/mixed-models/) the number of species is set to 5 and the number of sites is set to 3. I’m not surprised that the estimates of the random effect standard deviations would be well off the mark. I was surprised at the differences between my implementation in STan and the results from stan_lmer(). The difference in results must be because of a difference in the prior on the random effects, but as you’ll see if you look at the code in Github, I’ve tried several different prior options, and I can’t get them to match.

For my purposes, I’m only interested in the regression coefficients, but I’d appreciate any insight on the difference between the priors as I’ve implemented them in Stan and as they’re implemented in stan_lmer(). As I build up the more complex model I’m working on, I need to get the random effects right since the number of units at each level is so small.