# How to incorporate the uncertainty of the model coefficients in the prediction interval of a multiple linear regression?

I'm dealing with the modeling of small experimental physics data sets (specifically the stickiness of glue-compounds). As most experimental work does not generate thousands of samples, but rather a handful, I need to be inventive in how to deal with this small number of data sets (say 10-20). At this point I have a model-framework (regression see below at PSS) which can deal with this rather well.

However, to have a better picture of the accuracy of my predictions, I want to have an error-bar on my predicted values, this to check how well my predictions predict new experiments. As this work is numerical in nature, the error-bar will be originating from the underlying theoretical model, how do these errors propagate (i.e., error-analysis as one is used to in experimental physics)

For the sake of simplicity, assume that I am dealing with a multiple linear regression model, say (in reality there will be many many more terms): $$y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 \tag{1}$$ What I am looking for is an algebraic way of calculating (numerically) error-bars (in actuality its the prediction interval (PI) or confidence interval(CI), as both are related). In statistics literature, there are references to such a problem, and examples of how the PI and CI can be calculated. However, these only consider the variability of the $$x$$'s. The PI and CI are then related to (cf. question question 147242): $$\hat{V}_f=s^2\cdot\mathbf{x_0}\cdot\mathbf{(X^TX)^{-1}}\cdot\mathbf{x_0^T} + s^2 \tag{2}$$

In contrast to these, each of my model coefficients[see: PSS below] ($$\beta_0, \beta_1$$ and $$\beta_2$$) in this case have an error-bar (extracted via bootstrapping from a distribution, with the distributions being numerical in nature not analytic, and the distributions are specific for each of the three coefficients). Is there a way to incorporate the uncertainty of the $$\beta_i$$'s (c.q. the "error-bars") in the calculation of the PI (and CI).

To put it very simple, how can the equation $$\hat{V}_f=s^2\cdot\mathbf{x_0}\cdot\mathbf{(X^TX)^{-1}}\cdot\mathbf{x_0^T} + s^2 \tag{3}$$ be modified to also incorporate the fact that the coefficients themselves are a mean of a distribution.

(PS: One could create an ensemble of various model instances with the $$\beta_i$$ drawn from their respective distributions, and based on the distribution of obtained $$y_0$$ calculate the CI of the $$y_0$$, but this is not really computationally efficient and brings a lot of other issues which I would like to avoid.)

(PPS: The regression model presented is not the result of a direct regression toward a single data set, instead it is constructed as follows:

1. Create an ensemble of N data sets.
2. On each data set a regression gives rise to a linear model as indicated in the post above. This gives rise to N values for each of the coefficients $$\beta$$.
3. The mean of each of the three sets is calculated.
4. These three mean coefficients are the coefficients of the model presented above.
5. The goal here: find the prediction interval for the averaged model above taking into account the fact that the coefficients $$\beta$$ are calculated from numerical distributions.)
• This question should probably be posted here: stats.stackexchange.com Apr 1, 2020 at 19:28
• Reposting a question in the wrong forum, because an expert (such as whuber) pointed out that the question does not have enough focus is probably the wrong strategy. Apr 1, 2020 at 20:08
• I'm voting to close this question as off-topic because it doesn't appear to be about physics. Apr 2, 2020 at 8:03
• Meta post about the closing of this question. Apr 2, 2020 at 14:54
• I've removed a number of comments that were attempting to answer the question and/or responses to them. Please keep in mind that comments should be used for suggesting improvements and requesting clarification on the question, not for answering. Apr 10, 2020 at 6:00

I don't totally understand the post you linked, it seems they are implicitely assuming they have a model for how $$\vec{x}_0$$ is generated, which is not true in the generic case... However, if I understand your question, the most generic and simplest solution to achieve what you want is bootstrapping your prediction intervals. The basic idea is to use each of your $$N$$ sets of data to produce a vector $$\vec{\beta}$$, then stack your $$\vec{\beta}$$ into a matrix

$$B = \begin{bmatrix}\vec{\beta}_1 \\ \vec{\beta}_2 \\ \vdots \\ \vec{\beta}_N \end{bmatrix}.$$

Now your distribution of outputs is $$B\cdot\vec{x}_0$$, and you can do statistics on the elements of that vector present confidence intervals.

• The assumption is not implicit as it is a linear regression, which is performed according to a given algorithm. The fact that I use a Linear Regression on an ensemble to produce an equation of the same explicit form makes this a bit tricky. On the one hand I think you need to incorporate the error of the Linear regression model, and at the same time I need the error bars on the coefficients obtained from the ensemble. Apr 13, 2020 at 11:12
• The suggested bootstrapping is a valid idea (cf. the PS of my post) However, as I explained already, this is rather inefficient as you need to perform a large set of predictions, and need to keep track of the ensemble of models. Apr 13, 2020 at 11:13
• If you just have the distribution of $\vec{\beta}$, then since this is a multiple linear regression, the variance of the output over all $\vec{\beta}$ is just $$\sum \sigma_i^2 x_i^2,$$ where $x_i$ are the input in $\sigma_i^2$ is the variance of the $i$th parameter of $\vec{\beta}$ over the $N$ regressors. So you can compute the variance of each column of $B$ and use that to give you the variance of your output in only two steps instead of $N$. This only works for linear regression though, and it doesn't give you percentiles, which is what you wanted. Apr 13, 2020 at 19:34
• that looks interesting thanks. Sep 26, 2020 at 15:20

This is a problem that is essentially tailor-made for Bayesian analysis. The output of a Bayesian analysis is the joint distribution of all of your model coefficients. So, you can simulate samples from the predicted data by first drawing a sample from the model coefficients and then using those model coefficients to draw a sample from the data. This is called the "posterior predictive distribution". It is commonly used in Bayesian analysis to evaluate the validity of the model. If your model reasonably approximates your data generation process then your actual data should be reasonably similar to your posterior predicted data.

I recommend using the rstanarm package in R. IMO, even if you don't know R it is worth learning it just to use this package.

https://mc-stan.org/rstanarm/

• That probably works on a reasonable data set, but there are a few limitations: (1) it needs to be implemented within an existing python framework so R-libs are out of the question and (2) Indeed there are many ways were you can generate ensembles of predictions and continue from that, but that is the inefficient route I wish to avoid. It should be a single point calculation instead. Sep 26, 2020 at 15:23
• There is also PyStan which is closely related. Anyway, suit yourself regarding the computations. I believe that Bayesian statistics are the right approach for this type of problem, but it is your project
– Dale
Sep 26, 2020 at 16:17
• Thanks for the reference. I'll look into it. Sep 27, 2020 at 13:35

You should not mess your brain with statistics. There is Lies, Big Lies, and there is Statistics.

You should work on your direct task, what is causality of the effects which you obtain in your work.

We all know about that "Spurious Correlation" facts. Correlation is not causation. Stanley Cup is correlated with Staples sales[1]. So what? Nothing.

I don't understand why you need multiple linear regression, which is incredibly faulty because of internal theoretical inconsistencies. Mainly, there is no way you can use any result of any "regression" as proof of strong causality. But multiple mixed variable regression does not even let you find weak causality. You know what heteroscedasticity is? [2]

It is what 2003 Nobel Prize was awarded for. Work on physics, not statistics. You have Robert Engle for the second thing.

About Error bars which you need. Draw Error bars on paper with whatever size you feel is right. You are scientist. These are your bars, not somebody else's. Insert some noise inside your experimental signal line and conclude the error sizes you get.

• We are dealing with an ensemble model within the context of ML, and applied to a real world small data-set (i.e. no nice curve/distribution/... only numerical data). And we are not looking for a proof of strong or weak causality, just for a model that fits the data. Jan 25 at 10:39