Linear Regression Confidence Interval - MathOverflow most recent 30 from http://mathoverflow.net 2013-05-24T23:15:53Z http://mathoverflow.net/feeds/question/11907 http://www.creativecommons.org/licenses/by-nc/2.5/rdf http://mathoverflow.net/questions/11907/linear-regression-confidence-interval Linear Regression Confidence Interval Fernanda 2010-01-15T19:43:50Z 2010-07-09T00:38:41Z <p>I am fitting a linear regression line to my data and computing the confidence interval for some predicted $y$ (independent variable) (<a href="http://people.stfx.ca/bliengme/ExcelTips/RegressionAnalysisConfidence2.htm" rel="nofollow">http://people.stfx.ca/bliengme/ExcelTips/RegressionAnalysisConfidence2.htm</a>).</p> <p>Now I want to do the inverse. I need a way to measure the confidence for some predicted $x$ (dependent variable) using the same regression. How can I achieve this?</p> <p>Thank you!</p> http://mathoverflow.net/questions/11907/linear-regression-confidence-interval/11921#11921 Answer by kweinert for Linear Regression Confidence Interval kweinert 2010-01-15T21:21:42Z 2010-01-15T21:21:42Z <p>Usually you predict a dependent variable y and calculate a confidence interval, i.e. given x0, you calculate [y-, y+] where y will probably lie in.</p> <p>For the reverse, if you have a y0 and want to find [x-, x+] for whatever reasons, regression will not help.</p> <p>The appropriate tool for this kind of analysis could be structural equations modeling [1]</p> <p>[1] <a href="http://en.wikipedia.org/wiki/Structural_equation_modeling" rel="nofollow">http://en.wikipedia.org/wiki/Structural_equation_modeling</a></p> http://mathoverflow.net/questions/11907/linear-regression-confidence-interval/12192#12192 Answer by Simon BarthelmÃ© for Linear Regression Confidence Interval Simon BarthelmÃ© 2010-01-18T14:43:03Z 2010-01-18T14:43:03Z <p>It's an interesting question. There's a way to come up with a reasonable Bayesian answer, I think, but it requires a bit more work than linear regression. First, assume that you know the functional relationship that links y to x. That is, you know f in: y = f(x) + noise You make the usual assumption that the noise is normally distributed with mean 0 and variance s. That gives a probability distribution for y given f, x and s: p(y|x,f,s) So you can compute how likely it is to observer, say, y = 3.5, when <code>$f(x)=2x$</code>, and s=1, by computing: <code>$\mathcal{N}(3.5;2*x,1)$</code> where $\mathcal{N}(x;m,s)$ is the normal distribution function (mean m, std.dev s). Again, assuming both f and s are known exactly, the probability that $x$ was $a$ given the output $y$ is given by Bayes' theorem: $p(x=a|y,f,s)=p(y|x=a,f,s)p(x=a)/\int{p(y|x,f,s)p(x)dx}$</p> <p>(Note: $p(y|x=a,f,s)=\mathcal{N}(y;f(x),s)$)</p> <p>$p(x)$ is a prior distribution over the inputs, representing what you know about how inputs are generated, e.g. if you know that they are between 0 and 1 you can set $p(x)$ to the uniform distribution over [0,1]. Now $p(x|y,f,s)$ gives you a posterior distribution over possible inputs, so that you can find your confidence interval simply by plotting the distribution. </p> <p>The problem is that neither f nor s are known exactly. If you have tons of data, that assumption is not a problem: estimate f by non-parametric regression, and estimate s through the standard deviation of the residuals. If you don't have enough data, then your problem gets a lot more complicated - you have to "integrate out" f and s: $p(x|y)=\int{p(x|y,f,s)p(f,s)dfds}$. For that you need to compute a distribution over functions and noise levels, and integrate over it (tricky). Gaussian processes can help you there, read Rasmussen and Williams' book: <a href="http://www.gaussianprocess.org/gpml/" rel="nofollow">http://www.gaussianprocess.org/gpml/</a></p> http://mathoverflow.net/questions/11907/linear-regression-confidence-interval/31129#31129 Answer by Michael Hardy for Linear Regression Confidence Interval Michael Hardy 2010-07-08T22:38:42Z 2010-07-08T22:38:42Z <p>The question is ambiguous on an important point. The linked web page deals with a confidence interval for the regression line, and thus for the estimated expected value of the new observation $y$. But "predicted $y$" suggests something else: an interval for a new $y$ independent of those that went into the least-squares estimates. Usually the term "prediction interval" is used for that, to distinguish it from a <em>confidence interval</em> such as what the linked page treats. (OK, I see that prediction intervals are also treated on that page, at the bottom.) The prediction interval would be much longer than the confidence interval, because of the variability in the new observation (indeed, with a very large data set, the length of the confidence interval might be for all practical purposes zero, but that of the prediction interval would have a huge irreducible variance component: the variance of the error in each individual observation). (On the linked page, notice the $n$ in the denominators. Those terms go to zero as the sample-size grows. In the prediction interval there's one term with no $n$ in the denominator; the <em>confidence</em> interval there's not.)</p> <p>Let $\rho$ be the sample correlation; let $\overline{x}$ and $\overline{y}$ be the mean observed $x$- and $y$-values; let $s_x$ and $s_y$ be the sample standard deviations. Then the least-squares line for estimating the average $y$-value given the $x$-value is $$y - \overline{y} = \rho \frac{s_y}{s_x}(x - \overline{x}).$$ The line used for predicting $x$ based on $y$ is found similarly. If you know the slope of the line above and the sample correlation then you can find it.</p> <p>To finish the job, you need to know the sample SDs too.</p> http://mathoverflow.net/questions/11907/linear-regression-confidence-interval/31132#31132 Answer by Michael Hardy for Linear Regression Confidence Interval Michael Hardy 2010-07-08T23:36:45Z 2010-07-08T23:42:32Z <p>Just to recall some basic stuff: Suppose $Y_i = \alpha + \beta x_i + \varepsilon_i$ for $i=1,\dots,n$ and $\varepsilon_i \sim N(0,\sigma^2)$ and these are independent. We will observe the $x$'s and $y$'s but not $\alpha$, $\beta$, or $\varepsilon_i$. The $Y$'s are random only because the $\varepsilon$'s are random. These are weaker assumptions than that we have $(x,y)$ pairs that are <em>jointly</em> normally distributed.</p> <p>Let <code>$$X = \begin{bmatrix} 1 &amp; x_1 \\ \vdots &amp; \vdots \\ 1 &amp; x_n \end{bmatrix}$$</code> be the "design matrix" (so called because if the experimenter can choose the $x$ values, then this is how the experiment is designed).</p> <p>Then the least-squares estimates of $\alpha$ and $\beta$ are given by <code>$$\begin{bmatrix} \hat\alpha \\ \hat\beta \end{bmatrix} = (X^T X)^{-1}X^T Y$$</code> and therefore the probability distribution of the least-squares estimators is given by <code>$$\begin{bmatrix} \hat\alpha \\ \hat\beta \end{bmatrix} \sim N\left( \begin{bmatrix} \alpha \\ \beta \end{bmatrix}, \sigma^2 (X^T X)^{-1} \right)$$</code> (the variance is thus of course a $2\times 2$ positive-definite matrix). The predicted $y$-value for a given $x$ value is therefore $\hat\alpha + \hat\beta x$, and this therefore has a probability distribution given by <code>$$\hat y = \hat\alpha + \hat\beta x = \begin{bmatrix}1, &amp; x\end{bmatrix} \begin{bmatrix} \hat\alpha \\ \hat\beta \end{bmatrix} \sim N\left( \begin{bmatrix}1, &amp; x\end{bmatrix} \begin{bmatrix} \alpha \\ \beta \end{bmatrix}, \sigma^2 \begin{bmatrix} 1, &amp; x \end{bmatrix} (X^T X)^{-1} \begin{bmatrix} 1 \\ x \end{bmatrix} \right)$$</code> <code>$$= N\left( \alpha + \beta x, \frac{\sigma^2}{n}\frac{\sum_i (x_i - x)^2}{\sum_i (x_i - \overline{x})^2} \right).$$</code> (OK, check my algebra here; it's trivial but laborious.)</p> <p>There's a simple geometric intuition behind the dependence of the variance on $x$ and in particular the fact that the variance is smallest when $x=\overline{x}$, so think about that too.</p> <p>Now $\sigma^2$ must be estimated based on the data. The <em>errors</em> $y_i - (\alpha + \beta x_i)$ are unobservable but but the <em>residuals</em> $y_i - (\hat\alpha + \hat\beta x_i)$ (i.e. the <em>estimated</em> errors) are the components of the random vector $$\hat\varepsilon = (I - H)Y = (I - X(X^T X)^{-1} X^T)Y.$$ ("H" stands for "hat", for reasons that should be apparent.) It is easy to see that the $n\times n$ hat matirx $H = X(X^T X)^{-1} X^T$ is the matrix of the orthogonal projection of rank $2$ onto the 2-dimensional column space of $X$. And $I - H$ is the rank-$(n-2)$ projection onto the orthogonal complement of that space. Diagonalized, this matrix just has $n-2$ instances of 1 on the diagonal and 0 in the other two positions. Therefore $$\frac{\hat\sigma^2}{\sigma^2} = \frac{\| \hat\varepsilon \|^2}{\sigma^2}$$ is distributed like a sum of squares of $n - 2$ independent $N(0,1)$ random variables. It therefore has a chi-square distribution with $n-2$ degrees of freedom.</p> <p>Finally, we need this: $\hat\varepsilon$ and <code>$\begin{bmatrix}\hat\alpha \\ \hat\beta \end{bmatrix}$</code> are probabilistically independent. This is true because both are linear transformations of the same vector of independent identically distributed normal random variables and their covariance vanishes: <code>$$\operatorname{cov}\left(\hat\varepsilon, \begin{bmatrix}\hat\alpha \\ \hat\beta \end{bmatrix} \right) = \operatorname{cov}\left( (I - H)Y , (X^T X)^{-1}X^ T \right)$$</code> <code>$$= (I - H) \operatorname{cov}(Y,Y) X(X^T X)^{-1} = \sigma^2 (I - H) X(X^T X)^{-1}$$</code> and this is the $n\times 2$ zero matrix, by definition of $H$.</p> <p>Now all our lemmas are in place and we can draw some conclusions:</p> <p>Firstly <code>$$\frac{\hat y - (\alpha + \beta x)}{\sqrt{ \frac{\sigma^2}{n}\frac{\sum_i (x_i - x)^2}{\sum_i (x_i - \overline{x})^2} }} \sim N(0,1).$$</code></p> <p>Hence if $\sigma$ were miraculously known, we could say that <code>$$\hat y \pm A \sqrt{ \frac{\sigma^2}{n}\frac{\sum_i (x_i - x)^2}{\sum_i (x_i - \overline{x})^2} }$$</code> are the endpoints of a 90% confidence interval for $\alpha + \beta x$ if $\pm A$ are the endpoints of the interval above which lies 90% of the area under the bell-curve.</p> <p>But $\sigma$ is not known. Since $\hat\sigma^2$ is indpendent of the random variable in the numerator and has a chi-square distribution with $n-2$ degrees of freedom, we can put $\hat\sigma$ in place of $\sigma$ and instead of the normal distribution use the Student's t-distribution with $n-2$ degrees of freedom.</p> <p>That's the conventional frequentist confidence interval.</p> <p>For the prediction interval, just remember that the new value of $Y$ is independent of those we used above, so the variance of the difference between that and the predicted value is $\sigma^2$ plus the variance of the predicted value, found above.</p> http://mathoverflow.net/questions/11907/linear-regression-confidence-interval/31136#31136 Answer by Michael Hardy for Linear Regression Confidence Interval Michael Hardy 2010-07-09T00:38:41Z 2010-07-09T00:38:41Z <p>One little fact from algebra is worth mentioning here. Let $$\overline{x} = \frac{1}{n}\sum_{i=1}^n x_i.$$ Then $$\sum_{i=1}^n (x_i - x)^2 = n(x-\overline{x})^2 + \sum_{i=1}^n (x_i - \overline{x})^2.$$</p> <p>As a consequence, the confidence interval as I wrote it is the same as the one on the page that the original poster linked to.</p> <p>Another consequence is that the nature of the dependence of the variance of the predicted $y$-value on $x$ is seen to be quite simple.</p> <p>Another is that the sample mean is the maximum-likelihood estimate of the population mean when the population is assumed to be normally distributed.</p>