Berechnung der Wahrscheinlichkeit von x1> x2

Ich lerne selbst über die Wahrscheinlichkeit mithilfe von R, linearen Modellen und Wahrscheinlichkeitsberechnungen. Ich bin derzeit nicht sicher, wie ich zwei Vorhersagen aus einem Modell vergleichen kann. Die von mir verwendeten Daten werden (kostenlos) von hier heruntergeladen: 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

Das sind also die Setup-Daten. Nennen wir die Person mit dem höheren vorhergesagten CGPA (2,88) Rachel und die Person mit dem niedrigeren vorhergesagten CGPA (2,51) Tobias. Meine Frage ist, wie berechne ich die Wahrscheinlichkeit, dass Rachel einen höheren CGPA als Tobias hat? Ich habe in den Bereich unter der Kurve geschaut und bin mir nicht sicher, ob ich es richtig gemacht habe oder ob ich es richtig interpretiere. Flächenberechnung:

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

Der Unterschied zwischen den beiden Vorhersagen beträgt also ungefähr 12,5%. Wenn Rachel und Tobias jedoch dieselben Eingabevariablen hatten, um denselben CGPA zu erzeugen, beträgt die Wahrscheinlichkeit, dass einer von ihnen einen höheren CGPA aufweist, 50/50. Würde ich der Fläche 0,5 hinzufügen (62,5%), um die wahre Wahrscheinlichkeit zu erhalten? Bin ich weit weg und muss etwas anderes tun?

r regression predictive-models interpretation Kunio
quelle

Antworten:

Die Einstellung wird herkömmlicherweise in der Form ausgedrückt

y = X β + ε

$y = X\beta + \varepsilon$

für einen Vektor von Antworten, eine Modellmatrix und einen Vektor von Parametern unter der Annahme, dass die Zufallsfehler nicht mit gleichen Varianzen und Null bedeutet: das heißt $n$ $y$ $n\times k$ $X$ $k$ $\beta$ $\varepsilon = (\varepsilon_i)$ $\sigma^2$

E (ε) = 0; Var (ε) = σ^{2} I_{n} .

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

Wenn dies der Fall ist, ist die gewöhnliche Schätzung der kleinsten Quadrate

\hat{β} = (X^{'} X)^{-} X^{'} y .

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

Sei eine Matrix, deren Zeilen und die Werte der Regressoren für Rachel bzw. Thomas . Die vorhergesagten Antworten befinden sich im Vektor . Die tatsächlichen Antworten sind und wobei diese neuen Epsilons unkorrelierte Zufallsvariablen mit dem Mittelwert Null sind, unabhängig vom ursprünglichen und mit gemeinsamen Varianzen . $Z$ $2\times k$ $z_R$ $z_T$ $2$ $Z\hat\beta$ $z_R\beta+\varepsilon_R$ $z_T\beta+\varepsilon_T$ $\epsilon$ $\sigma^2$

Der Unterschied zwischen diesen Werten für Rachel minus Thomas, die ich nennen werde , ist einfach $\delta$

δ = (z_{R} β + ε_{R}) - (z_{T} β + ε_{T}) = (1, - 1) Z β + ε_{R} - ε_{T} .

$\delta=(z_R\beta+\varepsilon_R ) - (z_T\beta + \varepsilon_T) = (1,-1)Z\beta + \varepsilon_R - \varepsilon_T.$

Beide Seiten sind Matrizen - dh Zahlen - und offensichtlich sind sie aufgrund des Auftretens von auf der rechten Seite zufällig . (Die rechte Seite zeigt den geschätzten Unterschied zwischen den Antworten von Rachel und Thomas sowie die Abweichung zwischen den tatsächlichen und vorhergesagten Antworten von Rachel abzüglich der Abweichung zwischen den tatsächlichen und den vorhergesagten Antworten von Thomas.) Wir können den Erwartungsbegriff nach Begriff berechnen: $1\times 1$ $y$ $\varepsilon_R$ $\varepsilon_T$

\begin{aligned} E (δ) & = E ((1, - 1) Z β + ε_{R} - ε_{T}) \\ = (1, - 1) Z β + 0 - 0 \\ = z_{1} β - z_{2} β . \end{aligned}

$\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. }$

Dies ist genau das, was man annehmen würde: Der erwartete Unterschied ist der Unterschied in den vorhergesagten Werten. Sie kann geschätzt werden, indem die Parameter durch ihre Schätzungen ersetzt werden. Um dies anzuzeigen, setzen wir einen Hut über das " ": $E$

\begin{matrix} (1) & \hat{E} (δ) = (1, - 1) Z \hat{β} = z_{1} \hat{β} - z_{2} \hat{β} . \end{matrix}

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

Das ist die die in der Frage erscheint. $2.88-2.51$

Wir können die Analyse des Unterschieds zwischen Rachel und Thomas fortsetzen, indem wir die beiden Komponenten der Unsicherheit über diese Verteilung ausdrücken: Eine ist, weil und aus zufälligen Daten geschätzt werden, und die andere ist das Auftreten dieser zufälligen Abweichungen und . $\beta$ $\sigma$ $\varepsilon_R$ $\varepsilon_T$

\begin{matrix} (2) & \begin{aligned} Var (Rachel - Thomas) & = Var ((1, - 1) Z \hat{β} + ε_{R} - ε_{T}) \\ = (1, - 1) Z Var (\hat{β}) Z^{'} (1, - 1)^{'} + Var (ε_{R}) + Var (ε_{T}) \\ = (1, - 1) Z Var (\hat{β}) Z^{'} (1, - 1)^{'} + 2 {\hat{σ}}^{2} . \end{aligned} \end{matrix}

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

Die Varianzen der Epsilons werden durch geschätzt . Wir kennen da es von abhängt . Es ist Routine, diese Varianz zu schätzen, indem durch die Schätzung der kleinsten Quadrate , wodurch eine Menge erzeugt wird, die manchmal geschrieben wird . $\hat\sigma^2$ $\operatorname{Var}(\hat\beta)$ $\sigma$ $\sigma^2$ $\hat\sigma^2$ $\widehat{\operatorname{Var}}(\hat\beta)$

Diese Schätzungen können nur in Wahrscheinlichkeiten umgewandelt werden, indem spezifischere Annahmen über die bedingten Verteilungen von auf . $y$ $X$ Bei weitem am einfachsten ist es anzunehmen, dass multivariate Normal ist, denn dann ist (eine lineare Transformation des Vektors ) selbst Normal und daher bestimmen sein Mittelwert und seine Varianz seine Verteilung vollständig. Die geschätzte Verteilung wird erhalten, indem die Hüte auf und . $y$ $\delta$ $y$ $E$ $\operatorname{Var}$

Schließlich haben wir alle Informationen zusammengestellt, die für eine Lösung erforderlich sind. Das OLS-Verfahren schätzt die Verteilung von Rachels Antwort minus Thomas 'Antwort auf Normal mit einem Mittelwert, der der Differenz der vorhergesagten Werte und einer Varianz, die durch geschätzt wird und die geschätzte Fehlervarianz und die Varianz-Kovarianz-Matrix der Koeffizientenschätzungen . $(1)$ $(2)$ $\hat\sigma^2$ $\operatorname{Var}(\hat\beta)$

Dieser RCode führt direkt die in den Formeln und gezeigten Berechnungen durch : $(1)$ $(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

Die Ausgabe für diese Daten beträgt : OLS schätzt, dass die Wahrscheinlichkeit, dass Rachels CGPA die von Thomas übersteigt, bei . (In diesem Fall stellt sich heraus, dass Rachel und Thomas so ähnlich sind, das Modell so gut passt und die Datenmenge so groß ist, dass im Vergleich winzig ist bis und könnte so vernachlässigt werden. Das wird nicht immer der Fall sein.) $0.67$ $67\%$ $\widehat{\operatorname{Var}}(\hat\delta)$ $2\hat\sigma^2$

Dies ist der Mechanismus, der der Berechnung von Vorhersageintervallen zugrunde liegt : Mit dieser Verteilung können wir Vorhersageintervalle für die Differenz zwischen Rachels und Thomas 'CGPA berechnen.

whuber
quelle

@ Taylor Das Modell behauptet, dass jede einzelne Antwort in der Form vorliegt . Die Hüte werden nur angezeigt, wenn mit Modellschätzungen gearbeitet wird. Ich sehe, dass ich es verwirrend geschrieben habe - es ist ein Überbleibsel des Übergangs zwischen zwei Formulierungen des Modells. Lassen Sie mich das beheben und wir werden sehen, ob es dann konsistent aussieht.

z β + ϵ

$z\beta+\epsilon$

whuber

@whuber: frage: warum '-delta.hat' (negativ)? Und können wir pnorm durch ecdf {stats} durch eigenes geschätztes cdf ersetzen? Irgendwelche Implikationen für die lm-Schätzung? (lm nimmt keine Normalität an).

Maximilian

@Max (1) pnormberechnet die Wahrscheinlichkeit, dass eine Variable kleiner als ihr Argument ist, während wir die Chance haben möchten, größer zu sein . Technisch hätte ich mich also aufrufen sollen pnorm(0, delta.hat, sqrt(var.delta.hat), lower.tail=FALSE), aber ich habe seine Symmetrie ausgenutzt, um die Aussage zu verkürzen. (2) Es ist unklar, welche Werte Sie für Ihr ecdf vorschlagen. (3) Für nicht normale Antwortverteilungen benötigen Sie wahrscheinlich ein verallgemeinertes lineares Modell oder eine andere Verallgemeinerung.

whuber

Ihr Problem mag einfach klingen, ist aber überraschend kompliziert.

Um die Wahrscheinlichkeit , dass Rachels CPGA (nennen wir es zu beurteilen ) größer ist als Tobias' ( ), während zu wissen , was ihr , und -scores sind, ist das gleiche wie das Schreiben , wo sind ihre Punktzahlen. Da wir schreiben können, können wir auch sagen $y_1$ $y_2$ hgpasatltrs $P(y_2 - y_1 > 0 | X)$ $X$ $y_i = \hat{y_i} + \epsilon_i$

\begin{aligned} P (y_{2} - y_{1} > 0 | X) = & P (\underset{\sim N (0, 2 σ_{y}^{2})}{\underset{⏟}{ϵ_{2} - ϵ_{1}}} + \underset{= 2.8812 - 2.5082}{\underset{⏟}{\hat{y_{2}} - \hat{y_{1}}}} > 0 | X) \\ = & P (ϵ_{2} - ϵ_{1} < 0, 373) \end{aligned}

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

Hier stecken Sie fest, weil wir nicht genau kennen . Das Beste, was wir hier tun können, ist, es zu schätzen, indem wir die Varianz Ihrer Regressionsreste berechnen. Wenn Ihre Stichprobe groß genug ist ( ), konvergiert dies zu . $\sigma_y^2$ $\rightarrow \infty$ $\sigma_y^2$

Wenn Sie den Schätzfehler in ignorieren möchten , können Sie dies in R implementieren: $\hat{\sigma_y^2}$

sigma_hat <- summary(lm)$sigma
e2_min_e1 <- diff(predict(lm, new.df)) * -1

pnorm(e2_min_e1, 0, 2*sigma_hat)
# 0.6255

KenHBS
quelle

es ist nicht wahr, dass .

y_{i} = \hat{y_{i}} + ϵ_{i}

$y_i = \hat{y_i} + \epsilon_i$

Taylor

warum nicht? (eigentlich nur die lineare Projektion, aber unter den normalen linearen Regressionsannahmen ist dies auch die Bedingung exp) und es gilt dann immer, dass und epsilon hat den Mittelwert Null

\hat{y_{i}} \equiv E (y_{i} | X_{i})

$\hat{y_i} \equiv E(y_i | X_i)$

y_{i} = E (y_{i} | X_{i}) + ϵ_{i}

$y_i = E(y_i | X_i) + \epsilon_i$

KenHBS

{\hat{y}}_{i} = \hat{E (y_{i} | x_{i})}

$\hat{y}_i = \widehat{E(y_i|x_i)}$

Taylor

@ KenS. Danke Ken. Ich weiß, dass ich Standardfehler in 'Predict ()' erhalten kann, indem ich einfach 'se.fit = TRUE' hinzufüge. Ich habe es jedoch mit Ihrem Code versucht und es gab mir eine Fehlermeldung: 'Fehler in r [i1] - r [-Länge (r) :-( Länge (r) - Verzögerung + 1L)]: nicht numerisches Argument für Binär Betreiber '

Kunio

Eine der Standardannahmen von OLS ist, dass die lineare Funktionsform korrekt angegeben ist. Wenn diese Annahme zutrifft, ist . Ich bin mir nicht sicher, ob ich deinen Standpunkt verstehe. Könnte es nur ein Unterschied in der Notation sein?

y_{i} = E (y_{i} | X_{i}) + ϵ_{i}

$y_i = E(y_i | X_i) + \epsilon_i$

KenHBS