How is Group Level Term Estimated in Multilevel Model? I understand from reading Gelman and Hill that for a multilevel model such as this one
$$ y_i \sim N(\alpha_j + \beta X, \sigma^2) $$
$$ \alpha_j \sim N(\mu, \tau^2) $$
The $\alpha_j$ group-level intercept is a weighted average of the overall mean of the groups, $\mu$, and the deviation within the group (Gelman & Hill, p. 253). However, I am confused about how this $$ \alpha_j = \frac{(n_j/\sigma^2_y)\bar{y}_j + (1/\sigma^2_{\alpha})y_{all}}{(n_j/\sigma^2_y) + (1/\sigma^2_{\alpha})} $$
arises from the specification. Isn't $$ \alpha_j \sim N(\mu, \tau^2) $$ saying that you are taking a random draw from a normal distribution with mean $\mu$ and variance $\tau^2$, in which case, wouldn't each $\alpha_j$ be estimated randomly rather than from the data, even if $\mu$ is estimated as the average across all the responses for the different groups? Why does this specification make sense instead of simply expressing $\alpha_j$ as a weighted average?
 A: Mr. GoldIsFine,
It is not 100% clear to me if you're asking for a mathematical response or an intuitive response, so I'll do my best to give both.  I'll start with the parameter estimation, and the we'll dig into the intuition, and try to get at what is bothering you.
Let's start with the formulation of a MAP estimate.  The probability of our data ${(x_i, y_i)}$ given $(\alpha_j, \beta, \sigma, \tau, \mu)$, [note that I'm using $\sigma, \tau$ instead of $\sigma_y, \sigma_{alpha}$, both of which are used in the OP], is given by
$$
P(x,y | \alpha_j, \beta, \sigma, \tau, \mu) = P(x,y | \alpha_j, \beta, \sigma) = \prod_i \exp((y_i - \alpha_{j[i]} - \beta x_i)^2/\sigma^2),
$$
where $j[i]$ is the group with which data point $i$ is associated.  Note that this quantity doesn't depend on $\mu, \tau$ once $\alpha_j$ are fixed.  For the sake of brevity I will define $(x,y) = D$, $(\alpha_j, \beta, \mu) = \Theta$, and $(\tau, \sigma)$ = $\Pi$.  Then we have
$$
P(\Theta | D, \Pi) = P(D|\Theta, \Pi)\cdot P(\Theta | \Pi) / P(D | \Pi)
$$
This is the equation for the posterior distribution on $\Theta$ given the prior $P(\Theta | \Pi)$, and the observed data $D$.  We will attempt to maximize this with respect to $\Theta$.  Note that the denominator does not depend on $\Theta$ and is positive so we will ignore it.  Also note that we have already found the first factor above.  Now,
$$
P(\Theta | \Pi) = \prod_j \exp((\alpha_j - \mu)^2/\tau^2),
$$
so what we want to maximize finally (not really finally) is
$$
\prod_i \exp((y_i - \alpha_{j[i]} - \beta x_i)^2/\sigma^2) \cdot  \prod_j \exp((\alpha_j - \mu)^2/\tau^2).
$$
Since this function is clearly positive we can take logs and maximize this quantity instead.  Let $f(\Theta)$ be the logarithm of this quantity as a function of $\Theta$.  Since we are interested in $\alpha_j$ if we take the partial derivative of $f$ with respect to $\alpha_j$ and divide by 2 we arrive at:
$$
(\alpha_j - \mu)/\tau^2 + \sum_i^*(y_i - \alpha_j - \beta x)/\sigma^2,
$$
where the sum is taken over all $i$ corresponding to group $j$.  Setting this equal to zero and manipulating a bit we arrive at the equation:
$$
\alpha_j = \frac{\mu/\tau^2 + n_j\gamma_j/\sigma^2}{1/\tau^2 + n_j/\sigma^2},
$$
where $n_j$ is the number of data points in group $j$ and $\gamma_j$ is the mean of $y_i - \beta x_i$ taken over group $j$ (the group intercept).
Intuition:
Yes you are correct that if I tell you something is drawn from a normal distribution, you should guess that its true value is randomly somewhere in that distribution.  But now we have observed data!  And that data depended on the true value of $\alpha_j$.  So we can leverage this data to make a more informed estimation of $\alpha_j$.  The weighted average view of things is to give you sense of the effect that our priors have on our estimates.  Increasing $\tau$ forces the $\alpha_j$ closer to $\mu$, and decreasing $\tau$ allows the $\alpha_j$ to move more in the direction of the sample group intercept.
