CODEWITHSUNDEEP
rstan
goldisfine
goldisfine

Reputation: 4850

Nested Model in STAN?

Say I want to model a random effect at two levels, i.e. I have two levels of nesting: individuals within a parent group and parent groups within a grandparent group. I know how to write a basic model for a single random effect (below) from examples like these but I don't know how to write the equivalent to

lmer(resp ~ (1|a/b), data = DAT)

in lmer.

STAN code for single RE. Question is, how to nest a within a higher level b?

data{
  int<lower=0> N;
  int<lower=0> K;
  matrix[N,K] X;
  vector[N] price;
  int J;
  int<lower=1,upper=J> re[N];
}
parameters{
  vector[J] a;
  real mu_a;
  real tau;
  real<lower=0> sigma_a;
  real<lower=0> sigma;
  vector[K] beta;
}
transformed parameters{
  vector[N] mu_hat;
  for(i in 1:N)
    mu_hat[i] <- a[re[i]];
}
model {
  mu_a ~ normal(0,10);
  tau ~ cauchy(0,5);
  a ~ normal(mu_a,sigma_a);
  for(i in 1:N)
    price[i] ~ normal(X[i]*beta + mu_hat[i], sigma);
}
"

Upvotes: 6

Views: 1253

Answers (2)

Sid
Sid

Reputation: 440

If your model is simple, the brms package is worth looking at. It compiles your formulae down to stan and runs the model. It also has expressive syntax borrowed from lmer. What I like is that you can compile the model as a stan file and then build on it if you already have the lmer formulae

And of course it has the added benefit(from stan) that the confusing difference around estimation of "fixed effect" vs "random effects" is gone and both are estimated essentially as parameters with posterior distributions.

Upvotes: 0

Bob Carpenter
Bob Carpenter

Reputation: 3753

I'm not sure what the a/b notation is in lmer, but if you want nested levels multiple layers deep, then it's easy with a predictor. Say you have an IRT model with students (j in 1:J) nested in schools (school[j] in 1:S) and schools nested in cities (city[s] in 1:C).

[Update 14 April 2017]

You can now vectorize everything. So rather than this:

for (j in 1:J) 
  theta[j] ~ normal(alpha[school[j]], sigma_theta);
for (s in 1:S)
  alpha[s] ~ normal(beta[city[s]], sigma_alpha);
beta ~ normal(0, 5);

you can have

theta ~ normal(alpha[school], sigma_theta);
alpha ~ normal(beta[city], sigma_alpha);
beta ~ normal(0, 5);

Upvotes: 4

Related Questions