# Posterior predictive regression lines in Bayesian linear regression seem to exclude most data points

I have some data to which I'd like to fit a GLM using PyMC3. Here are the posterior predictive regression lines plotted along with the data points:

Domain knowledge indicates that a linear model with zero-intercept is appropriate and the data looks linearly related. However, it feels as though the fit is not OK since the post predictive lines seem to cover only a few data points, leaving most out. Is this intuition correct? In this case, 64% of the data points are not covered by the regression lines.

The PyMC3 model I used is really simple:

with pm.Model() as geom_glm:
pm.GLM.from_formula('delta_p ~ 0 + sot', dados)
step = pm.Metropolis()
trace = pm.sample(draws = 10**4, step = step)

pm.plot_posterior_predictive_glm(trace,
lm = lambda x, sample: sample['sot'] * x,
label = 'posterior predictive check',
samples = 1000)


If I allow the intercept to be inferred from the data as well (even though domain knowledge indicates this is inappropriate), the situation doesn't seem to improve much:

It's nothing to worry about. Bayesian linear regression is (to simplify) just linear regression with random variables for the slope and intercept.

Now, the posterior distributions for the slope and intercept random variables permit many credible values for those parameters. And those combinations of credible parameters permit many credible regression lines to be drawn:

(Source for stylised image: http://www.indiana.edu/~jkkteach/WorkshopUWM2012.html)

The shaded area on your chart is simply a similar representation for the above.

I would argue the opposite to you: I would argue that you want the shaded region to be as tight as possible, to have more "confidence" (in the non-technical sense of the word) in the fitted regression line.

• I’d like the shaded region to be as tight as the data (and the prior) allow, but no further. My question is whether the data allow it to be as tight as it currently is. In fact, the chart you showed also strikes me as “too tight” (lots of data points outside the shaded area). I’d like to know if my intuition about this is off. – LmnICE Feb 3 '18 at 3:34
• The regression line represents the conditional mean of $Y$ given $X$. There's no reason that the credible area for such means should cover any data points at all. Again, it's not a problem. – A. G. Feb 3 '18 at 7:06
• No, you're already predicting $Y$ using the conditional mean: $\hat{Y} \approx E[Y|X]$. This is the very core of the linear regression model - it's a model for the mean, and this is how it is defined in statistics and machine learning textbooks. But there's nothing wrong with this because the conditional mean is the optimal predictor for $Y$ under a squared loss. This is why it is so widely used. Here's a proof: stats.stackexchange.com/a/71869/190932 – A. G. Feb 5 '18 at 10:34
• Yes, I understand it can be confusing. Sure, you are indeed trying to model the relationship between the dependent variable $Y$ and covariate(s) $X$. But linear regression can only predict $E[Y|X]$. Have a look at this stylised representation: fao.org/docrep/w5449e/w5449egf.gif If you want a model between $Y$ and $X$ that works in a different way, then you should pick a different model. – A. G. Feb 5 '18 at 11:36
• It should theoretically be possible to do without changing the model, like in this rJAGs example: doingbayesiandataanalysis.blogspot.co.uk/2015/10/… But I'm not sure how you'd do it in PyMC3, as I'm no expert in this particular library. Perhaps you can try asking a follow-up question about this, quoting this link, and see if anyone has any better ideas? – A. G. Feb 5 '18 at 12:38