How to deal with treatment effect at state level but observations at school level? I have a panel of schools and I want to estimate the effect of a state level policy on school level outcomes (all schools in a treated state are treated). I can't use a fixed effects model since the treatment is at the state level and so the fixed effects term absorb all the treatment effect. I can't use a difference in difference model model since I do not have many pre-treatment periods. I am also worried about heterogeneity between the schools.
Lately, I have been thinking about multi level models (schools nested in states). What do you think about that? Any other suggestions on how to take on this issue?
 A: I take it you're measuring at the school level and that


*

*Schools are nested within states

*States are either intervened upon or not


In what follows, I'll just assume the likelihoods are normal for simplicity.  Under the assumption that states are exchangeable, then I think we can model this as follows...
The outcome for school $i$ in state $j$ is 
$$y_i \sim \mathcal{N}( \mu_j, \sigma^2_j)$$
That is, the outcomes in each state can be considered to be normally distirbuted with some state level average outcome $\mu_j$.  I've allowed for the outcomes to be heteroskedastic by also allowing the states to have their own variance.  The appropriateness of this assumption is up to you and can be validated via a posterior predictive check.
We can assume that the state level effects are drawn from some country level distribution.  Namely
$$ \mu_j \sim \mathcal{N}(M_0 + x_jM_1, S^2) $$
Here, $M_0$ is the country level average outcome, $x$ is a binary indicator for if the state was intervened upon or not,  $M_1$ is the effect of the intervention, and $S^2$ the country level variance.  It might be the case that the effect of the intervention may depend on the state, in which case we could include a random effect of the intervention.  Again, this will depend on your expert knowledge of the problem and can be evaluated via posterior predictive checks.
We can then place priors on $M_0$, $M_1$, $S^2$ and $\sigma^2_j$.
Fitting this model should be straightforward in something like Stan or brms.  I can't give you an entire tutorial on either in a CV answer, but I can point you to appropriate resources should you find this approach interesting.  I would warn you that approaching the problem from this way might be very involved.
EDIT: So Because I'm bored, I went ahead and simulated some data to fit this model.  Here is some Stan code to simulate the data.
data{
   // Number of observations.  Could be number of schools
   // Depends on your problem
   int N_obs;
   // Number of states; 51
   int N_states;
   // x is a binomial variable; 1 -- intervene, 0 -- control
   real x[N_states];
}
parameters{}
model{}
generated quantities{
  real M0 = normal_rng(80,5);
  real M1 = fabs(normal_rng(0,2));
  real S = fabs(cauchy_rng(0,1));

  vector[N_states] sigma;
  vector[N_states] mu;
  matrix[N_states,N_obs] y;

  for (j in 1:N_states){
    mu[j] = normal_rng(M0 + x[j]*M1, S );
    sigma[j] = lognormal_rng(log(5.0),0.5);
    for(i in 1:N_obs){
      y[j,i] = normal_rng(mu[j], sigma[j]);
    }
  }
}

If you run this Stan code, you might get something that looks like this

I've cooked the books to ensure that the intervention (i.e. $M_1$) is positive, so if you squint you can see that the light blue dots are, on average, higher than the dark blue dots.
Now, for the modelling piece, I am going to do something different than what I've written.  I'm actually going to average the observations across states and instead posit that I've observed a draw from that state's sampling distribution.  That means that my likelihood is
$$ \bar{y}_j \sim \mathcal{N}(\mu_j, s^2) $$
where $s$ is the standard error for the mean of the observations within each state.
The Stan code for the model is
data{
  int N_states;
  int N_obs;

  vector[N_states] x;

  real sey[N_states];
  real ybar[N_states];
}
parameters{
  real M0;
  real M1;
  real<lower=0> S;
  // Random effects
  vector[N_states] z;
}
transformed parameters{
  vector[N_states] mu;
  //Vectorized means for each state
  mu = M0 + x*M1 + z*S;

}
model{
  M0 ~ normal(80,5);
  M1 ~ normal(0,1);
  S ~ cauchy(0,1);
  z ~ normal(0,1);


  ybar ~ normal(mu, sey);

  }

Once we run the code on the data generated above, we can plot the marginal posterior of the intervention

The black bar is the true value of the intervention from the simulation.  As you can see, the model does pretty well!
