In lecture notes on introductory graduate course on Bayesian statistics, there is a short discussion of how Multiple linear regression may be treated in the paradigm of "borrowing strength" aka "hierarchical models".

The outline of the treatment is given in the lecture note, and the rest is left as an easy exercise for the reader, which I am not able to complete.

We are given the following multiple regression:

$$Y_{[1 \times n]} = \beta_{[1 \times p]} X_{[p \times n]}+ \varepsilon_{[1 \times n]}\,,$$

where the subscript indicates the dimensionality of the quantities, $Y$ is known (response) $X$ is known (predictors) (both are centered so we are not concerned with the intersect), we want to estimate $\beta$ vector.

The errors $$\{\varepsilon_i\}_{i=1}^n\,~\text{are i.i.d drawn from}~~N \left(0, \frac{1}{\tau}\right)\,,$$ with $N(\mu, \sigma^2)$ being the normal distribution (that is the variance of $\varepsilon_i$ is $\frac{1}{\tau}$).

The system is assumed to be under-determined, meaning $p > n$. Using "hierarchical approach" in this situation is claimed to be reasonable. Specifically, we assume that the elements of vector $\beta$ are distributed as

$$ \beta_i \sim N\left(0, \frac{\sigma^2}{\tau}\right)\quad \left(\text{i.e. the variance is } \frac{\sigma^2}{\tau}\right)\,. $$

Here things start to become a bit vague. The claim is that one can use "empirical Bayes approach", that is to estimate "hyper parameters" $\tau$ and $\sigma^2$ using MLE, which is claimed to be easy, that would give estimates $\hat{\tau}, \hat{\sigma^2}$. And then the posterior distribution of $\beta | \hat{\sigma}, \hat{\tau}, X, Y$ is available in closed form. I was not able to carry this scheme out, my (long) attempts are detailed below. Any help or reference to a solution would be appreciated.

What I try to do is the following

  1. write likelihood function $f_{Y_i}(y_i| \tau, \sigma; \{x_{ij}\})$

  2. maximize with respect to $\tau, \sigma^2$, to find MLEs $\hat{\tau}$, $\hat{\sigma^2}$. (I get stuck here)

  3. plug in the MLE estimates from the above and write posterior distribution.

(Painful) Details:

So I want to write the likelihood function $L(\tau, \sigma | \{x_{ij}\}, y_i)$, for this I need the probability density $f_{Y_i}(y_i| \tau, \sigma)$. The given model implies that

$$ Y_i = \sum_{j=1}^p\beta_j x_{i,j} + \varepsilon_i\,, $$

thus, $$ Y_i| \tau, \sigma \sim N\left(0, \text{var}\left[{\sum_{j=1}^p\beta_j x_{i,j}}\right] + \text{var}[\varepsilon_i]\right) = N\left(0, \frac{\sigma^2}{\tau}\left(\sum_{i=1}^p x_{i,j}^2\right) + \frac{1}{\tau}\right) = N\left(0, \frac{\sigma^2}{\tau}\sigma^2_{i} + \frac{1}{\tau}\right)\,, $$ where I used a shorthand $\sigma^2_{i} :=\sum_{i=1}^p x_{i,j}^2$.

The probability density of $Y_i$ is

$$ f(y_i | \tau, \sigma, \{x_{i,j}\}) = \frac{\sqrt{\tau}}{\sqrt{2 \pi \left( \sigma^2\sigma^2_{i} + 1 \right)}} \exp\left[\frac{\tau y_i^2}{2 \left( \sigma^2\sigma^2_{i} + 1 \right)}\right]\,. $$

This yields the log likelihood function

$$ L(\tau, \sigma^2: \{x\}, Y) = \sum_{i=1}^n \log\left(\frac{\sqrt{\tau}}{\sqrt{2 \pi \left( \sigma^2\sigma^2_{i} + 1 \right)}} \exp\left[-\frac{\tau y_i^2}{2 \left( \sigma^2\sigma^2_{i} + 1 \right)}\right]\right) $$ $$ = 0.5 n \cdot \log(\tau) - 0.5 \sum_{i=1}^n \log\left( \sigma^2\sigma^2_{i} + 1 \right) - 0.5 n \cdot \log(\pi) - \sum_{i=1}^n \frac{\tau y_i^2} { 2 \left( \sigma^2\sigma^2_{i} + 1 \right) }\,.\tag{1} $$

The MLE for $\tau$ is easy to get (derivative $\frac{\partial}{\partial \tau}$ and equating to zero).

$$ 0.5 n \cdot \frac{1}{\tau} - \sum_{i=1}^n \frac{y_i^2} { 2 \left( \sigma^2\sigma^2_{i} + 1 \right)} = 0\,, $$

$$ \boxed{ \hat{\tau} = \frac{2}{n}\left( \sum_{i=1}^n \frac{y_i^2} { 2 \left( \sigma^2\sigma^2_{i} + 1 \right)} \right)^{-1}}\,. $$

For the MLE of $\sigma^2$, derivative of $(1)$, $\frac{\partial}{\partial \sigma^2}$ yields

$$ - \sum_{i=1}^n \frac{\sigma^2_{i}}{\left( \sigma^2\sigma^2_{i} + 1 \right)} + \sum_{i=1}^n \frac{\tau y_i^2 \sigma^2_{i} } { \left( \sigma^2\sigma^2_{i} + 1 \right)^2\,. } $$ Rearranging and equating to zero $$ \left( \sum_{i=1}^n \frac{\sigma^2_{i}}{\left( \sigma^2\sigma^2_{i} + 1 \right)} - \sum_{i=1}^n \frac{\tau y_i^2 \sigma^2_{i} } { \left( \sigma^2\sigma^2_{i} + 1 \right)^2 } \right) = 0\,. $$

Here I get stuck I don't see how I can solve for $\sigma^2$ (in the denominators). All the development seems to be too cumbersome, contrary to what the claim in the note is - that it is easy.

  • $\begingroup$ Numerical optimization? $\endgroup$ – jaradniemi Aug 23 '16 at 16:44
  • $\begingroup$ I do strongly believe there must be an analytical solution to the original problem. I have probably overlooked some straight forward simple book solution. I thought that may be one can re-introduce $\beta$ and then find the $f(y|\tau , \sigma^2)$ thorough integrating this out $\int_\beta f(y|\tau , \sigma^2, \beta) f(\beta|\sigma^2, \tau ) d\beta$ and make use of that, but so far I am not able to solve this way either... $\endgroup$ – them Aug 23 '16 at 16:54
  • $\begingroup$ You have to be careful doing this. Sometimes the empirical Bayes solution says, "I found you a nice, big optimum for the likelihood: just don't penalize the coefficients at all!" $\endgroup$ – eric_kernfeld Sep 8 '17 at 14:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.