I am building a hierarchical linear model with varying intercepts. It takes the form for each unit $i$ in group $j$:

$$y_{ij} = \alpha_j + \beta_1 x_{ij,1} + \beta_2 x_{ij,2} \quad (1) $$

I am developing this hierarchical linear model using a complete Bayesian Analysis using stan. In stan, I am using fairly non-informative priors. All variables -- for each $j$, $\alpha_j$, $\beta_1$, and $\beta_2$ -- are given normal distributions with hyper-parameters. Additionally the output $y$ follows $N(\hat{y}, \sigma_y)$ where $\hat{y} = y_{ij}$ from above.

These hyper-parameters are all given uniform distributions over $[0,\infty)$ with the exception of

$$\sigma_y \sim U(0, 100)$$ $$\sigma_a \sim U(0, 100)$$ $$\sigma_b \sim U(0, 100)$$

By contrast, if I build a linear regression model with complete pooling, then the model takes the form

$$y_{ij} = \alpha + \beta_1 x_{ij,1} + \beta_2 x_{ij,2} \quad(2) $$

Question: The non-informative decisions for priors has the effect that the parameter space for $\alpha_j$, for each $j$, $\beta_1$, and $\beta_2$ should consider the coefficients obtained from regression in (2). Then it should follow that hierarchical linear regression performs as well as ordinary regression. Running this in py-stan, I find that the answer to this assertion is "no." Why?

Edit: By "perform as well as", I mean than when comparing the RMSE between predictions using (1) and (2), (2) always outperforms (1). In fact, the opposite is seen in the scatter plot below which displays for, years $y \in [2000,2015]$, the RMSE of (1) and (2) trained on data from years $< y$ and then tested on data from year $y$.

enter image description here


Here is a related example that shares the same spirit of the phenomenon I expect to observe in this example:

Suppose that a regression model of $k$ predictors trained on an output $y$ has $R^2 = G_k$. If we add one predictor $p$ to this model, making the model have a total of $k+1$ predictors, then let $R^2 = G_{k+1}$. Since the set of solution(s) $S_2$ to linear regression in the second instance (with predictor $p$) is a superset of the first regression model's solution space $S_1$, it follows that $G_{k+1} \geq G_{k}$.

That is, if $\beta_1, \beta_2$ are the predictors for the first and second regression problems, then

$$ G_{k} = \frac{\sum_{j=1}^{n} (\beta_1^Tx_j - y_j)^2}{\sum_{j=1}^{n} (\bar{y} - y_j)^2} \leq \frac{\sum_{j=1}^{n} (\beta_2^Tx_j - y_j)^2}{\sum_{j=1}^{n} (\bar{y} - y_j)^2} = G_{k+1} $$

  • $\begingroup$ what do you mean by "performs as well as"? $\endgroup$ – Ben Bolker Mar 28 '18 at 13:29
  • $\begingroup$ I was unclear when I wrote that. I have made that bit more clear. Please see my edit. $\endgroup$ – Muno Mar 28 '18 at 13:35
  • $\begingroup$ when you compute RMSE, what are you using as your point estimate? The mean? I think the phenomenon you're seeing could have something to do with the dispersion of the posterior, which would be larger in the more complex model (cf. the "penalty" terms in BIC and DIC) $\endgroup$ – Ben Bolker Mar 28 '18 at 14:19
  • $\begingroup$ Yes, each $y_{ij}$ has a posterior mean. I use that as my point estimate. $\endgroup$ – Muno Mar 28 '18 at 15:59

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.