# Calculating Probability of x1 > x2

I'm self-learning about probability using R, linear models, and probability calculations. I'm currently stuck on how to compare 2 predictions from a model. Data I'm using is downloaded(free) from here: wmbriggs.com/public/sat.csv

df <- read.csv("sat.csv")              # Load data
lm <- lm(cgpa~hgpa+sat+ltrs,data=df)   # model to predict College GPA
new.df <- data.frame(hgpa=c(4,3),sat=c(1168,1168),ltrs=c(6,6))  # 2 scenario data. Same SAT and LTRS, differing Highschool GPA
predict(lm,new.df)                     # plug our scenario data into the model to predict cgpa based on input
1        2
2.881214 2.508154


So that's the setup data. Let's name the person with the higher predicted CGPA (2.88) Rachel and the person with the lower predicted CGPA (2.51) Tobias. My question is, how do I calculate the probability that Rachel has a higher CGPA than Tobias? I've looked into area under the curve, and not sure if I did it's correct, or if I'm interpreting it correctly. Area calculation:

area <- pnorm(2.881214,1.9805,0.7492264)-pnorm(2.508154,1.9805,0.7492264) # area under the curve between the 2 predicted CGPAs
[1] 0.1259893


So the difference between the 2 predictions is 12.5% roughly. However, if Rachel and Tobias had the same input variables to produce the same CGPA, the probability that either 1 of them has a higher CGPA is 50/50. Would I add 0.5 to the area (62.5%) to get the true probability? Am I way off and need to be doing something else?

The setting is conventionally expressed in the form

$$y = X\beta + \varepsilon$$

for an $n$-vector $y$ of responses, an $n\times k$ model matrix $X$, and a $k$-vector of parameters $\beta$, under the assumptions that the random errors $\varepsilon = (\varepsilon_i)$ are uncorrelated with equal variances $\sigma^2$ and zero means: that is,

$$E(\varepsilon)=0; \ \operatorname{Var}(\varepsilon) = \sigma^2 I_{n}.$$

When this is the case, the ordinary least squares estimate is

$$\hat\beta = (X^\prime X)^{-} X^\prime y.$$

Let $Z$ be a $2\times k$ matrix whose rows $z_R$ and $z_T$ give the values of the regressors for Rachel and Thomas, respectively. The predicted responses are in the $2$-vector $Z\hat\beta$. The actual responses are $z_R\beta+\varepsilon_R$ and $z_T\beta+\varepsilon_T$ where these new epsilons are zero-mean uncorrelated random variables, independent of the original $\epsilon$, and with common variances $\sigma^2$.

The difference between those values for Rachel minus Thomas, which I will call $\delta$, is simply $$\delta=(z_R\beta+\varepsilon_R ) - (z_T\beta + \varepsilon_T) = (1,-1)Z\beta + \varepsilon_R - \varepsilon_T.$$

Both sides are $1\times 1$ matrices--that is, numbers--and evidently they are random by virtue of the appearance of $y$ on the right hand side. (The right hand side is the estimated difference between Rachel's and Thomas's responses, plus the deviation $\varepsilon_R$ between Rachel's actual and predicted responses, minus the deviation $\varepsilon_T$ between Thomas's actual and predicted responses.) We may compute its expectation term by term:

\eqalign{ E(\delta) &= E\left((1,-1)Z\beta + \varepsilon_R - \varepsilon_T\right)\\ &= (1,-1)Z\beta +0 - 0\\ &= z_1\beta - z_2\beta. }

This is exactly what one would suppose: the expected difference is the difference in predicted values. It can be estimated by replacing the parameters by their estimates. To indicate this, let's place a hat over the "$E$":

$$\hat{E}(\delta) = (1,-1)Z\hat\beta = z_1\hat\beta - z_2\hat\beta.\tag{1}$$

That's the $2.88-2.51$ appearing in the question.

We may continue the analysis of the difference between Rachel and Thomas by expressing the two components of uncertainty about that distribution: one is because $\beta$ and $\sigma$ are estimated from random data and the other is the appearance of those random deviations $\varepsilon_R$ and $\varepsilon_T$.

\eqalign{ \operatorname{Var}(\text{Rachel}-\text{Thomas}) &= \operatorname{Var}\left((1,-1)Z\hat\beta + \varepsilon_R - \varepsilon_T\right) \\ &= (1,-1)Z \operatorname{Var}(\hat\beta) Z^\prime (1,-1)^\prime + \operatorname{Var}(\varepsilon_R) + \operatorname{Var}(\varepsilon_T) \\ &=(1,-1)Z \operatorname{Var}(\hat\beta) Z^\prime (1,-1)^\prime + 2\hat\sigma^2.\tag{2} }

The variances of the epsilons are estimated by $\hat\sigma^2$. We don't know $\operatorname{Var}(\hat\beta)$ because it depends on $\sigma$. It is routine to estimate this variance by replacing $\sigma^2$ by its least-squares estimate $\hat\sigma^2$, producing a quantity sometimes written $\widehat{\operatorname{Var}}(\hat\beta)$.

These estimates can be converted into probabilities only by making more specific assumptions about the conditional distributions of $y$ on $X$. By far the simplest is to assume $y$ is multivariate Normal, for then $\delta$ (being a linear transform of the vector $y$) itself is Normal and therefore its mean and variance completely determine its distribution. Its estimated distribution is obtained by placing the hats on $E$ and $\operatorname{Var}$.

Finally we have assembled all the information needed for a solution. The OLS procedure estimates the distribution of Rachel's response minus Thomas's response to be Normal with a mean equal to the difference in predicted values $(1)$ and with a variance estimated by $(2)$, which involves the estimated error variance $\hat\sigma^2$ and the variance-covariance matrix of the coefficient estimates, $\operatorname{Var}(\hat\beta)$.

This R code directly carries out the calculations exhibited in formulas $(1)$ and $(2)$:

fit <- lm(cgpa ~ hgpa + sat + ltrs, data=df)         # model to predict College GPA
Z <- as.matrix(data.frame(intercept=1, hgpa=c(4,3), sat=c(1168,1168),ltrs=c(6,6)))

cont <- matrix(c(1,-1), 1, 2)             # Rachel - Thomas "contrast".
beta.hat <- coef(fit)                     # Estimated coefficients for prediction
delta.hat <- cont %*% Z %*% beta.hat      # Predicted mean difference
sigma.hat <- sigma(fit)                   # Estimated error SD
var.delta.hat <- cont %*% Z %*% vcov(fit) %*% t(Z) %*% t(cont) + 2 * sigma.hat^2
pnorm(0, -delta.hat, sqrt(var.delta.hat)) # Chance Rachel > Thomas


The output for these data is $0.67$: OLS estimates that there is a $67\%$ chance that Rachel's CGPA exceeds that of Thomas. (It turns out in this case, because Rachel and Thomas are so similar, the model fits so well, and the amount of data is so large, that $\widehat{\operatorname{Var}}(\hat\delta)$ is tiny compared to $2\hat\sigma^2$ and so could be neglected. That will not always be the case.)

This is the mechanism that underlies the computation of prediction intervals: we can compute prediction intervals for the difference between Rachel's and Thomas's CGPA using this distribution.

• @Taylor the model asserts that any individual response is in the form $z\beta+\epsilon$. The hats appear only when working with model estimates. I see that I have written it confusingly--it's a vestige of making a transition between two formulations of the model. Let me fix it and we'll see whether it looks consistent then.
– whuber
Aug 15, 2017 at 21:41
• @whuber: question: why '-delta.hat' (negative)? And can we replace pnorm with own estimated cdf via ecdf {stats}? Any implication for the lm estimation? (lm does not assume normality). Aug 18, 2017 at 8:14
• @Max (1) pnorm computes the chance that a variable will be less than its argument whereas we want the chance of being greater. Technically, then, I should have invoked pnorm(0, delta.hat, sqrt(var.delta.hat), lower.tail=FALSE), but I exploited its symmetry to shorten the statement. (2) It's unclear what values you propose for your ecdf. (3) For non-Normal response distributions you likely would need a generalized linear model or some other generalization.
– whuber
Aug 18, 2017 at 13:09

Your problem may sound easy, but it's surprisingly complicated.

In order to assess the probability that Rachel's CPGA (call it $y_1$) is larger than Tobias' ($y_2$), while knowing what their hgpa, sat and ltrs-scores are, is the same as writing $P(y_2 - y_1 > 0 | X)$, where $X$ are their scores. Because we can write $y_i = \hat{y_i} + \epsilon_i$, we can also say

\begin{align*} P(y_2 - y_1 > 0 | X) = & P( \underbrace{\epsilon_2 - \epsilon_1}_{\sim N(0, 2\sigma_y^2)} + \underbrace{\hat{y_2} - \hat{y_1}}_{ = 2.8812 - 2.5082} > 0 | X) \\ = &P(\epsilon_2 - \epsilon_1 < 0,373) \end{align*}

This is where you get stuck, because we don't know $\sigma_y^2$ for sure. The best we can do here, is estimate it by calculating the variance of your regression residuals. If your sample is large enough ($\rightarrow \infty$), this will converge to $\sigma_y^2$.

If you want to ignore the estimation error in $\hat{\sigma_y^2}$, you can implement this in R:

sigma_hat <- summary(lm)$sigma e2_min_e1 <- diff(predict(lm, new.df)) * -1 pnorm(e2_min_e1, 0, 2*sigma_hat) # 0.6255  • it is not true that$y_i = \hat{y_i} + \epsilon_i$. Aug 15, 2017 at 20:02 • why not?$\hat{y_i} \equiv E(y_i | X_i)$(actually just the linear projection, but under the normal linear regression assumptions, this is also the cond. exp) and it always holds then that$y_i = E(y_i | X_i) + \epsilon_i$and epsilon has zero mean Aug 15, 2017 at 20:05 •$\hat{y}_i = \widehat{E(y_i|x_i)}$Aug 15, 2017 at 20:09 • @KenS. Thank you Ken. I know I can get Standard Error in the 'predict()' by simply adding 'se.fit=TRUE'. I tried it with your code though and it gave me an error message: 'Error in r[i1] - r[-length(r):-(length(r) - lag + 1L)] : non-numeric argument to binary operator' Aug 15, 2017 at 20:20 • One of the standard assumption of OLS is that the linear functional form is correctly specified. If that assumption holds, then$y_i = E(y_i | X_i) + \epsilon_i\$. I'm not sure I am getting your point. Could it just be a notational difference? Aug 15, 2017 at 20:23