Calcular la Probabilidad de que x1 > x2

Question

Soy auto-aprendizaje acerca de la probabilidad mediante R, modelos lineales, y cálculos de probabilidad. Actualmente estoy atascado en la forma de comparar 2 predicciones a partir de un modelo. Los datos que estoy usando es descargar(gratis) desde aquí: 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

Así que esa es la configuración de datos. Deje que el nombre de la persona con el mayor predijo CGPA (2.88) Rachel y la persona con el menor predijo CGPA (2.51) Tobias. Mi pregunta es, ¿cómo puedo calcular la probabilidad de que Rachel tiene una mayor CGPA de Tobias? He mirado en el área bajo la curva, y no se si hice lo correcto, o si estoy interpretando correctamente. Cálculo del área:

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

Así que la diferencia entre los 2 predicciones es de 12,5% aproximadamente. Sin embargo, si Rachel y Tobias tenido las mismas variables de entrada para producir la misma CGPA, la probabilidad de que 1 de ellos tiene una mayor CGPA es de 50/50. Podría yo añadir 0,5 para el área (62.5%) para obtener el verdadero probabilidad? Estoy lejos y necesita hacer algo más?

Answer 1

El valor convencionalmente se expresa en la forma

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

for an $n$-vector $s$ 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\veces 1$ matrices--that is, numbers--and evidently they are random by virtue of the appearance of $s$ 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) y= 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.\la etiqueta{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$ y por lo tanto podrían ser descuidado. Que no siempre será el caso.)

Este es el mecanismo que subyace en el cálculo de los intervalos de predicción: podemos calcular intervalos de predicción para la diferencia entre Rachel y Thomas CGPA el uso de esta distribución.

Answer 2

El problema puede parecer fácil, pero es sorprendentemente complicado.

Con el fin de evaluar la probabilidad de que Rachel CPGA ( $y_1$ ) es mayor que Tobias' ($y_2$), mientras que conocer cuál es su hgpa, sat y ltrs-de los puntajes, es la misma que la escritura $P(y_2 - y_1 > 0 | X)$ donde $X$ son sus puntuaciones. Porque se puede escribir $y_i = \hat{y_i} + \epsilon_i$, también podemos decir

\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*}

Aquí es donde te quedas atascado, porque no sabemos $\sigma_y^2$ seguro. Lo mejor que podemos hacer aquí, es estimado mediante el cálculo de la varianza de la regresión de los residuos. Si la muestra es lo suficientemente grande ($\rightarrow \infty$), esto va a converger a $\sigma_y^2$.

Si desea omitir la estimación de error en $\hat{\sigma_y^2}$, se puede implementar esta en 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