Hello everybody,

This is my first Stan question—I’m still learning, so please excuse me if the question is basic: this is about non-centered parameterization in a hierarchical setting. I want to know if my modelling was done correctly.

I was wondering about an extension of the 8-school example. Suppose we have three schools instead with classes nested in schools, and we want to estimate treatment effects per class given this knowledge. So for example, let the following be our data:

```
schools_df <- data.frame(
school = c("A", "A", "A", "B", "B", "C", "C", "C"),
classes = c(1, 2, 3, 1, 2, 1, 2, 3),
y = c(28, 8, -3, 7, -1, 1, 18, 12),
sigma = c(15, 10, 16, 11, 9, 11, 10, 18))
```

I set up the following Stan program and would like to know if I’ve approached it correctly:

```
data {
int N; // Number of observations
int<lower=0> J; // Number of schools (now 3)
int<lower=0> g[N]; // Group assignment
real y[N]; // Estimated treatment effects
real<lower=0> sigma[N]; // s.e. of effect estimates
}
parameters {
vector[J] mu; // By-school mean
real<lower=0> tau[J]; // By-school s.d.
real eta[N]; // Individual class errors
}
transformed parameters {
real theta[N]; // School effects
for (n in 1:N)
theta[g[n]] = mu[g[n]] + tau[g[n]] * eta[n];
}
model {
target += normal_lpdf(eta | 0, 1);
target += normal_lpdf(y | theta, sigma);
}
```

Also, in the case that we would want to estimate school effects pooled within school across classes, the following should be implemented—right?.

```
data {
int N; // Number of observations
int<lower=0> J; // Number of schools (now 3)
int<lower=0> g[N]; // Group assignment
real y[N]; // Estimated treatment effects
real<lower=0> sigma[N]; // s.e. of effect estimates
}
parameters {
real mu; // Population mean
real<lower=0> tau; // Population s.d.
real eta[J]; // school-level errors
}
transformed parameters {
real theta[J]; // School effects
for (j in 1:J)
theta[j] = mu + tau * eta[j];
}
model {
target += normal_lpdf(eta | 0, 1);
for(n in 1:N){
y[n] ~ normal(theta[g[n]], sigma[n]);
}
}
```

(I’m not too sure how to do the `target+=`

expression for the `y`

here.)

The original Stan model that I used as a baseline, the non-centered parameterization when we had 8 schools (so one per observation) was the following:

```
data {
int<lower=0> J; // Number of schools
real y[J]; // Estimated treatment effects
real<lower=0> sigma[J]; // s.e. of effect estimates
}
parameters {
real mu; // Population mean
real<lower=0> tau; // Population s.d.
real eta[J]; // school-level errors
}
transformed parameters {
real theta[J]; // School effects
for (j in 1:J)
theta[j] = mu + tau * eta[j];
}
model {
target += normal_lpdf(eta | 0, 1);
target += normal_lpdf(y | theta, sigma);
}
```

I think I’m doing it right, but I’m still confused by multiple things, so I’d like to clear up my confusions as I go.

Thank you so much.