Cálculo da probabilidade de x1> x2

Sou autodidata sobre probabilidade usando R, modelos lineares e cálculos de probabilidade. Atualmente, estou preso em como comparar duas previsões de um modelo. Os dados que estou usando são baixados (grátis) a partir daqui: 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

Então esses são os dados de configuração. Vamos nomear a pessoa com o CGPA previsto mais alto (2.88) Rachel e a pessoa com o CGPA previsto mais baixo (2.51) Tobias. Minha pergunta é: como faço para calcular a probabilidade de Rachel ter um CGPA maior que Tobias? Examinei a área sob a curva e não tenho certeza se fiz corretamente ou se estou interpretando corretamente. Cálculo de á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

Portanto, a diferença entre as duas previsões é de 12,5% aproximadamente. No entanto, se Rachel e Tobias tivessem as mesmas variáveis de entrada para produzir o mesmo CGPA, a probabilidade de um deles ter um CGPA maior é 50/50. Eu adicionaria 0,5 à área (62,5%) para obter a verdadeira probabilidade? Estou longe e preciso fazer outra coisa?

r regression predictive-models interpretation Kunio
fonte

Respostas:

A configuração é expressa convencionalmente na forma

y = X β + ε

$y = X\beta + \varepsilon$

para um vetor de respostas , uma matriz modelo e um vetor de parâmetros , sob as premissas de que os erros aleatórios não estão correlacionados com variâncias iguais e zero significa: isto é, $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}.$

Nesse caso, a estimativa dos mínimos quadrados ordinários é

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

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

Seja uma matriz cujas linhas e fornecem os valores dos regressores para Rachel e Thomas, respectivamente. As respostas previstas estão no vetor . As respostas reais são e onde esses novos epsilons são variáveis aleatórias não correlacionadas com média zero, independentes do original e com variações comuns . $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$

A diferença entre esses valores para Rachel menos Thomas, que chamarei de , é simplesmente $\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.$

Ambos os lados são matrizes - isto é, números - e, evidentemente, são aleatórios em virtude da aparência de no lado direito. (O lado direito é a diferença estimada entre as respostas de Rachel e Thomas, mais o desvio entre as respostas reais e previstas de Rachel, menos o desvio entre as respostas reais e previstas de Thomas.) Podemos calcular seu termo de expectativa por termo: $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. }$

Isso é exatamente o que se poderia supor: a diferença esperada é a diferença nos valores previstos. Pode ser estimado substituindo os parâmetros por suas estimativas. Para indicar isso, vamos colocar um chapéu sobre o " ": $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}$

Esse é o aparece na pergunta. $2.88-2.51$

Podemos continuar a análise da diferença entre Rachel e Thomas, expressando os dois componentes da incerteza sobre essa distribuição: um é porque e são estimados a partir de dados aleatórios e o outro é a aparência desses desvios aleatórios e . $\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} }$

As variações dos epsilons são estimadas por . Não conhecemos porque depende de . É rotineiro estimar essa variação substituindo por sua estimativa de mínimos quadrados , produzindo uma quantidade às vezes escrita . $\hat\sigma^2$ $\operatorname{Var}(\hat\beta)$ $\sigma$ $\sigma^2$ $\hat\sigma^2$ $\widehat{\operatorname{Var}}(\hat\beta)$

Estas estimativas podem ser convertidos em probabilidades apenas por fazer suposições mais específicos sobre as distribuições condicionais de sobre . $y$ $X$ De longe, o mais simples é assumir que é Normal multivariada, pois então (sendo uma transformação linear do vetor ) em si é Normal e, portanto, sua média e variância determinam completamente sua distribuição. Sua distribuição estimada é obtida colocando os chapéus em e . $y$ $\delta$ $y$ $E$ $\operatorname{Var}$

Finalmente, reunimos todas as informações necessárias para uma solução. O procedimento OLS estima que a distribuição da resposta de Rachel menos a resposta de Thomas seja Normal com uma média igual à diferença nos valores previstos e com uma variação estimada por , que envolve a variação de erro estimada e a matriz de variância-covariância das estimativas do coeficiente, . $(1)$ $(2)$ $\hat\sigma^2$ $\operatorname{Var}(\hat\beta)$

Este Rcódigo realiza diretamente os cálculos exibidos nas fórmulas e : $(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

A saída para esses dados é : a OLS estima que há uma probabilidade de que o CGPA de Rachel exceda o de Thomas. (Neste caso, porque Rachel e Thomas são muito parecidos, o modelo se encaixa muito bem e a quantidade de dados é tão grande que é minúsculo em comparação para e, portanto, pode ser negligenciado. Isso nem sempre será o caso.) $0.67$ $67\%$ $\widehat{\operatorname{Var}}(\hat\delta)$ $2\hat\sigma^2$

Este é o mecanismo subjacente ao cálculo dos intervalos de previsão : podemos calcular intervalos de previsão para a diferença entre o CGPA de Rachel e Thomas usando essa distribuição.

whuber
fonte

@ Taylor, o modelo afirma que qualquer resposta individual está no formato . Os chapéus aparecem apenas ao trabalhar com estimativas de modelo . Vejo que escrevi de maneira confusa - é um vestígio de fazer uma transição entre duas formulações do modelo. Deixe-me corrigi-lo e veremos se parece consistente.

z β + ϵ

$z\beta+\epsilon$

whuber

@whuber: question: why '-delta.hat' (negativo)? E podemos substituir o pnorm pelo próprio cdf estimado via ecdf {stats}? Alguma implicação para a estimativa lm? (lm não assume normalidade).

Maximilian

O @Max (1) pnormcalcula a chance de uma variável ser menor que seu argumento, enquanto queremos a chance de ser maior . Tecnicamente, então, eu deveria ter invocado pnorm(0, delta.hat, sqrt(var.delta.hat), lower.tail=FALSE), mas explorei sua simetria para encurtar a afirmação. (2) Não está claro quais valores você propõe para seu ecdf. (3) Para distribuições de resposta não normais, você provavelmente precisaria de um modelo linear generalizado ou de alguma outra generalização.

whuber

Seu problema pode parecer fácil, mas é surpreendentemente complicado.

A fim de avaliar a probabilidade de que CPGA de Rachel (chamemos-lhe ) é maior do que Tobias' ( ), sabendo que os seus , e -scores são, é o mesmo que escrever , onde são suas pontuações. Como podemos escrever , também podemos dizer $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*}$

É aqui que você fica preso, porque não sabemos ao . O melhor que podemos fazer aqui é estimar calculando a variação de seus resíduos de regressão. Se sua amostra for grande o suficiente ( ), isso convergirá para . $\sigma_y^2$ $\rightarrow \infty$ $\sigma_y^2$

Se você deseja ignorar o erro de estimativa em , pode implementá-lo em R: $\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
fonte

não é verdade que .

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

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

Taylor

Por que não? (na verdade apenas a projeção linear, mas sob as premissas normais de regressão linear, esta também é a condição exp) e sempre é válido que e epsilon tem média zero

\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. Obrigado Ken. Eu sei que posso obter o erro padrão no 'predict ()' simplesmente adicionando 'se.fit = TRUE'. Eu tentei com o seu código e ele me deu uma mensagem de erro: 'Erro em r [i1] - r [-length (r) :-( length (r) - lag + 1L)]: argumento não numérico para binário operador '

Kunio 15/08

Uma das suposições padrão do OLS é que a forma funcional linear está especificada corretamente. Se essa suposição , então . Não tenho certeza se estou entendendo o seu ponto. Poderia ser apenas uma diferença notacional?

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

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

KenHBS