# Calculate variance of a prediction for TLS

## Definitions

I have a set of data. Let's assume the true underlying model to be $$\eta_i=\beta_0+\beta_1\zeta_i$$. Where $$\beta_i$$ are the true coefficients of the model, $$\eta_i$$ are the true independent variables and $$\zeta_i$$ are the dependent variables. These cannot be observed.

The measurements (the data I have) are assumed to both have random errors ($$x_i = \zeta_i + \delta_i$$ and $$y_i = \eta_i + \gamma_i$$), which we assume to be normally distributed: $$\delta_i \sim N(0,\sigma_{x_i}^2)$$ and $$\gamma_i \sim N(0,\sigma_{y_i}^2)$$.

Through values from datasheets and error propagation, I have estimates of this distrubution: $$s_{x_i} \approx \sigma_{x_i}$$ and $$s_{y_i} \approx \sigma_{y_i}$$.

## Problem

Because my measurements have errors on both axis, I use Total Least Squares (TLS) to fit a linear model to my data: $$\hat y_i = b_0 + b_1 x_0$$.

Confusingly, I found people calling this also Orthogonal Distance Regression (ODR) or Bivariate Least Squares (BLS). I use the scipy.odr (which is using ODRPACK) implementation to be exact.

Now, let's assume I have a new independent data point which I assume follows the same model ($$\eta_0=\beta_0+\beta_1\zeta_0$$). Again, this true new data point cannot be observed, so I actually have $$x_0$$. Again, through looking at data sheets and propagating the error, I estimated $$s_{x_0} \approx \sigma_{x_0}$$.

I use the fitted model to predict $$\hat y_0 = b_0 + b_1 x_0$$. Now I want to calculate the variance of the prediction $$s_{\hat y_0}^2$$. How do I do that?

I don't think this is relevant to the problem, but I include it to avoid an XY Problem.

I want to use $$s_{y_0}$$ to construct a prediction interval for $$x_0$$: $$\hat y_0 \pm t_{\alpha,n-2}s_{\hat y_0}$$.

Then I check if the measured $$y_0$$ falls outside of the prediction interval to detect an anomaly in the data. I need the probability, as one needs to know the probability that this new data point is an anomaly. I asked about this method in this question.

## My journey

I already found a bunch of implementations, all of which give different results!

#### Approach 1

This question leads to these lecture notes and this website. The website references Wolberg, J., Data Analysis Using the Method of Least Squares, 2006, Springer for their approach.

I don't fully understand how they end up at their solution, I think they propagate the error of the prediction $$\hat y_0 = b_0 + b_1 x_0$$ using the Taylor expansion method, but I am unsure. These sources don't use different $$s_{y_i}$$. They assume one $$s_y$$ for all data points, including the data point to be predicted. The code calculates: $$s_{\hat y_0}^2 = s_y^2 + s_{b_0}^2 + x_0^2 s_{b_1}^2 + 2x_0\text{cov}(b_0, b_1)$$ If I read Wolberg, Section 2.6 correctly, this method does not take into account the errors in $$x$$, but I am not sure.

#### Approach 2

I also came across this question, which found this paper.

Here, they somehow expand the implementation of the normal Least Squares method to arrive at: $$s_{\hat y_0}^2 = [1 + \pmb X_0^T(\pmb X^T \pmb W^{-1} \pmb X)^{-1} \pmb X_0 + s_{x_0}^2b_1^2]s^2$$ where $$\pmb W$$ is calculated from the weights for the TLS, and $$s$$ is the experimental error.

#### Approach 3

The same paper from Approach 2 claims that the same value can be also derived by another method, and references Meloun M, Militky J, Forina M. Chemometrics for Analytical Chemistry. Volume 1: PC-aided Statistical Data Analysis. Ellis Horwood: Chichester, 1992. To my knowledge they also propagate the error of the prediction using the Taylor expansion method. But they arrive at: $$s_{\hat y_0}^2 = s^2 + s_{b_0}^2 + x_0^2s_{b_1}^2 + b_1^2s_{x_0}^2+ 2x_0\text{cov}(b_0, b_1)$$ Interestingly, this approach matches the approach 2, if I set $$s=1$$.

#### Side notes

The paper from Approach 2 and 3 actually calculate the variance of the response mean value (which in my eyes can be used for the confidence interval) with $$s_{\hat y_0}^2 = [\pmb X_0^T(\pmb X^T \pmb W^{-1} \pmb X)^{-1} \pmb X_0 + s_{x_0}^2b_1^2]s^2$$ or $$s_{\hat y_0}^2 = s_{b_0}^2 + x_0^2s_{b_1}^2 + b_1^2s_{x_0}^2+ 2x_0\text{cov}(b_0, b_1)$$. They then calculate the variance of the response variable (used for prediction interval) as the mean of $$q$$ observations with $$s_{\hat y_0}^2 = [\frac1q+\pmb X_0^T(\pmb X^T \pmb W^{-1} \pmb X)^{-1} \pmb X_0 + s_{x_0}^2b_1^2]s^2$$. As I only have one observation, $$q=1$$. And as approach 2 should match approach 3, I used the same logic for approach 3 to get the prediction interval.

Additionally, when researching what could have gone wrong, I looked at Melound et. al., Equation (1.40). And if I am not mistaken, the second term has a second partial derivative, which might be neglected in the paper? But I am very sure, that all of these sources have much more knowledge in statistics, so I am pretty sure I simply don't understand what's going on.

• Perhaps if you would clearly characterize what you mean by "prediction interval" and the precise meaning of "$x_0$" then your question would answer itself. After all, if you wish to predict a response from values of explanatory variables, it does you no good to use the true values of the explanatory variables, because you don't know them. If, then, you are using the measured values, then their errors don't appear to be relevant for this purpose: you want to run a regression conditional on the measured values.
– whuber
Commented Jul 24 at 20:07
• @whuber My bad, I thought a prediction interval is a well defined term. I mean the interval, in which a new measurement of the same distribution will fall, with a specified probability. This should take the variance of the process as well as the the variance of the model into account. Commented Jul 25 at 12:46
• Also, all values I defined are the measured values. As you said, I don't know the real values. But I don't see how the errors of the measured values should not play a role in the prediction interval? If the true value is actually a bit shifted, then the response will be different as well. I don't understand how a regression conditional would help me. Commented Jul 25 at 12:48
• "Prediction interval" indeed has a standard clear definition. But it's not evident you are using it in a standard way. Indeed, I still can't figure out what you really want. If you are making the prediction based on a measured value, then you need to regress the response against the measured values, period, with no consideration of measurement errors in the explanand. If you want an unconditional prediction interval, then the value of the explanand is irrelevant.
– whuber
Commented Jul 25 at 12:49
• Sorry, I don't understand where the confusion is coming from. In short, I want to calculate what the paper is calculating. Commented Jul 25 at 13:14