After having performed a regression, how can I estimate the probability of a specific outcome?

Let's assume that I have a dataset which gives me the two variables: height ($x_1$), daily calorie intake ($x_2$) and weight ($y$) of a person. In this dataset, we assume I have a large enough number ($n$) of examples in my dataset to make it statistically accurate.

Now, I can do a regression with multiple variables (in this case: 2), to obtain a prediction of a person's weight based on their height and calorie intake. It's rather easy: If I want to find out the predicted weight of a person, I need simply insert the value $x_1$ and $x_2$ into my fitted equation to obtain $y$.

However, I am wondering how I can do something a bit different:

Given a height $x_1$ and a daily calorie intake $x_2$, what is the probability of a person having a weight greater than or less than a certain value, or having a weight in a certain range.

For example: Assume I have $2000$ calories intake per day and a body height of $180 \text{ cm}$, what is the probability that the person weighs between $80\text{ kg}$ and $90\text{ kg}$.

Is such an estimation possible? And if so, how can I do it (analytically or numerically)? Is there any way for me to bound the error that this estimate will have?

• How about estimating the distribution of residuals and reading the answer from that? Nov 4 '14 at 9:57

If you're taking a frequentist approach, it sounds you're trying to compute the distribution that's used in a prediction interval.

You can do that with a regression, but you do rely on the distributional assumption more than for other forms of inference.

The usual way is to construct a version of a pivotal quantity, but there are some differences from the usual definition of pivotal quantity because instead of trying to get a pivot for a parameter, you need something that will work for a future observation.

In regression, if $y_f$ is that future observation and $\hat{y}_f=x_f\hat{\beta}$ its predicted value, then $Q=\frac{y_f-\hat{y}_f}{s_f}$ has a t-distribution with $n-p-1$ d.f. where $p$ is the number of predictors and $s_{f}$ is the standard error of the prediction.

The formula for $\text{Var}(y_f-\hat{y}_f)$ is often given in derivations of the prediction interval; $s_f$ is the square root of its sample estimate. If $x_f$ is the row-vector of predictors (including the constant) for the future observation, $y_f$, then

$\text{Var}(y_f-\hat{y}_f) = \sigma^2(1+x_f (X'X)^{-1} x_f')$

To obtain $s_f$, replace $\sigma^2$ by $s^2$ in the above formula and take the square root.

You can make probability statements about $Q$, and from those, derive statements about $y_f$ in terms of $\hat{y}_f$ and $s_f$ that should have the appropriate coverage if the assumptions hold.

Packages that will give one-sided prediction intervals can be used to help with the sort of computation you need, but it's easy enough to do directly.

• Dear Glen_b, could you maybe define what y_f^hat is? I admit I am somewhat of a beginner when it comes to statistics (regressions are not exactly high end), and I am not familiar with the notation. How can I calculate it? Nov 5 '14 at 10:49
• I've made a brief addition to the answer. Hopefully that's clearer now. Nov 5 '14 at 10:57

You can do it easily if you know the distribution of the error terms (or make an assumption). In case you describe, normality could be reasonable. Then $Y \sim N(\beta_1 x_1 + \beta_2 x_2, \sigma^2)$. Lets say that $\beta_1=0.4$, $\beta_2=0.004$ and $\sigma^2 = 100$. Then for the height and intake you describe, the weight is distributed as $Y \sim N(80, 10^2)$, and hence probability that his weight is between 80-90 kilos is roughly $\Phi(1) - \Phi(0)\approx 0.34$). But be advised that there are multiple assumptions which should be checked so that this inference makes sense, such as constant variance and the normality or errors.

• Thank you for your comment! There is just one thing I don't quite understand: By your definition, the Normal distribution is centered around the output of the function I compute. This, so far, I understand. Edit: Sorry, hit enter by accident. Now, am I correct in assuming that I estimate sigma^2 by summing over y_i - y'_i / n-p-1? (with n = training examples, p = parameters, in this case 2, y_i known (correct) outputs and y'_i the output of my function) Second: Is there any way for me to give an error bound to this? i.e.: A person is 80-90kilos with 0.34 +/- c? Now, Nov 4 '14 at 11:32
• I think error bounds on probability estimates makes little sense: the lack of certainty is already encoded in the statement. Nov 5 '14 at 14:25