7.2 Regressão Linear Simples
A palavra “regressão” não constitui escolha feliz sob o ponto de vista etimológico, mas já está tão arraigada na literatura estatística que não tentamos substituí-la por outra que exprimisse mais convenientemente suas propriedades essenciais. Foi introduzida por Galton, em conexão com a hereditariedade da estatura. Galton verificou que os filhos de pais cujas estaturas se afastam \(x\) polegadas em relação à estatura média de todos os pais, apresentam um desvio em relação à estatura média de todos os filhos inferior a \(x\) polegadas, isto é, verifica-se o que Galton chamou de “regressão à mediocridade”. Em geral, a idéia ligada à palavra “regressão” nada tem a ver com êste sentido, devendo ser considerada apenas como expressão convencional. (Yule and Kendall 1948, 246)
Para maiores detalhes veja (Galton 1889, 95–99) e (Stigler 1997).
7.2.1 Reta
A equação da reta é uma relação matemática utilizada para descrever uma reta no plano cartesiano. Pode ser apresentada de formas distintas, sendo que na Seção 7.2.2 é utilizada a notação da reta reduzida, fazendo \(a' = \beta_{1}\) e \(b' = \beta_{0}\).
| Tipo | Equação |
|---|---|
| Geral | \(ax + by + c = 0\) |
| Segmentária | \(\frac{x}{-c/a} + \frac{y}{-c/b} = 1\) |
| Reduzida | \(y = -\frac{a}{b} x - \frac{c}{b} \; \therefore \; y = a'x + b'\) |
Exemplo 7.7 (Equação reduzida) Considere a reta que passa pelos pontos \(A = (0,-3)\) e \(B = (1.5,0)\).
Uma maneira de descobrir a equação reduzida é substituir os pontos \(A\) e \(B\) em \(y = a'x + b'\):
| Tipo | Equação |
|---|---|
| Ponto \(A\) | \(-3 = a' \times 0 + b' \; \therefore \; b' = -3\) |
| Ponto \(B + b'\) | \(0 = a' \times 1.5 + (-3) \; \therefore \; a' = \frac{3}{1.5} = 2\) |
Assim, a equação reduzida da reta é \(y = 2x - 3\), onde a inclinação é \(a' = 2\) e o intercepto (constante ou coeficiente linear) é \(b' = -3\). Para cada aumento de 1 unidade em \(x\), \(y\) aumenta 2 unidades.
Exemplo 7.8 (Equação segmentária) Do Exemplo 7.7, pode-se obter a equação segmentária da reta a partir da forma reduzida. \[y = 2x - 3 \; \therefore \; 2x - y = 3 \; \therefore \; \frac{2}{3} x - \frac{1}{3} y = \frac{3}{3} \; \therefore \; \frac{x}{3/2} + \frac{y}{3/-1} = 1 \; \therefore \; \frac{x}{1.5} + \frac{y}{-3} = 1.\] Assim, \(-c/a = 1.5\) e \(-c/b = -3\). Note que \(x_{B} = 1.5\) e \(y_{A} = -3\).
Exemplo 7.9 (Equação geral) Do Exemplo 7.7, pode-se obter a equação geral da reta a partir da forma reduzida. \[y = 2x - 3 \; \therefore \; 2x - y - 3 = 0.\] Assim, \(a = 2\), \(b = -1\) e \(c = -3\). Note que \(a' = - \frac{2}{-1} = 2\) e \(b' = - \frac{-3}{-1} = -3\).
7.2.2 Modelo
O modelo de Regressão Linear Simples (RLS) populacional/universal é construído, na abordagem clássica, com todos os \(N\) pares ordenados da população ou universo, i.e., \((x_i,y_i), \; i \in \{1,\ldots,N\}\). Pode ser descrito pela Eq. (7.13).
\[\begin{equation} Y = \beta_0 + \beta_1 x + \varepsilon, \tag{7.13} \end{equation}\] onde \(\varepsilon \sim \mathcal{N}(0,\sigma_{\varepsilon})\).
Note que por simplicidade considera-se \(Y \equiv Y_i\), \(x \equiv x_i\) e \(\varepsilon \equiv \varepsilon_i\). É importante ainda ter em mente que a variável aleatória \(Y\) representa a esperança condicional aos valores observados de \(x\), denotada por \(E(Y|x)\). A estimativa de \(E(Y|x) = Y\) será anotada por \(\widehat{E(Y|x)} = y\).
Exemplo 7.10 No gráfico a seguir pode-se notar que a linha azul refere-se à média \(\bar{y}=180\). Este valor é constante – ou incondicional – para qualquer valor de \(x\) (no caso a temperatura máxima do dia na escala Celsius). Já a linha vermelha pontilhada indica que espera-se um valor baixo de \(y\) (no caso o número de garrafas vendidas no dia) quanto menor o \(x\), e alto \(y\) quanto maior o \(x\).
plot(dr$temp, dr$gar, las = 1)
abline(h = mean(dr$gar), col = 'blue')
abline(a = -19.11, b = 5.91, col = 'red', lty = 2)
legend(26, 235, legend=c('Incondicional','Condicional'),
col = c('blue', 'red'), lty = 1:2, cex = 0.8, box.lty = 0)
7.2.2.1 Estimativa dos parâmetros
Na maioria dos casos práticos trabalha-se com amostras, sendo necessário estimar os valores de \(\beta_0\) e \(\beta_1\). O método dos mínimos quadrados pode ser utilizado para calcular estas estimativas. O princípio do método é minimizar a soma de quadrado dos erros, i.e., \[\begin{equation} \text{minimizar} \sum_{i=1}^n \varepsilon_{i}^{2}. \tag{7.14} \end{equation}\]
Basicamente utiliza-se \(\varepsilon = Y - \beta_0 - \beta_1 x\) da Eq. (7.13) e deriva-se (parcialmente) em relação a \(\beta_0\) e \(\beta_1\), fazendo cada uma das derivadas parciais igual a zero. Para maiores detalhes recomenda-se (DeGroot and Schervish 2012, 689). As estimativas por mínimos quadrados são enfim dadas por \[\begin{equation} \hat{\beta}_1 = \frac{n \sum{xy} - \sum{x} \sum{y}}{n \sum{x^2} - (\sum{x})^2} \tag{7.15} \end{equation}\] e \[\begin{equation} \hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}. \tag{7.16} \end{equation}\]
Exemplo 7.11 Considere novamente os dados do Exemplo 7.2. Pode-se ajustar um modelo linear simples conforme Equações (7.13), (7.15) e (7.16). Note a estimativa na forma \(y = -19.11 + 5.91 x\).
##
## Call:
## lm(formula = gar ~ temp, data = dr)
##
## Coefficients:
## (Intercept) temp
## -19.110 5.915Modelo RPO
Um modelo na forma da Eq. (7.17) é chamado Regressão Pela Origem (RPO) pelo fato de a reta ajustada passar pelo ponto \((0,0)\), a origem do plano cartesiano. Basta fazer \(\beta_0=0\) na Equação (7.13), e para mais detalhes veja (Eisenhauer 2003). \[\begin{equation} Y = \beta_1^{RPO} x + \varepsilon, \tag{7.17} \end{equation}\] onde \(\varepsilon \sim \mathcal{N}(0,\sigma_{\varepsilon})\). A estimativa do parâmetro \(\beta_1^{RPO}\) é dada por \[\begin{equation} \hat{\beta}_1^{RPO} = \frac{\sum{xy} }{\sum{x^2} }. \tag{7.18} \end{equation}\]
Exemplo 7.12 Considere novamente os dados do Exemplo 7.2. Pode-se ajustar um modelo RPO conforme Equações (7.17) e (7.18), resultando em \(y = 5.36 x\). Na função stats::lm (linear model) do R é possível fazer o ajuste utilizando + 0 no lugar de - 1.
##
## Call:
## lm(formula = gar ~ temp - 1, data = dr)
##
## Coefficients:
## temp
## 5.359Exercício 7.7 Considere novamente os dados dos Exemplos 7.2 e 7.11.
a. Calcule à mão a estimativa dos parâmetros do modelo descrito na Equação (7.13) considerando as Equações (7.15) e (7.16).
b. Calcule à mão a estimativa dos parâmetros do modelo descrito na Equação (7.17) considerando a Equação (7.18).
\(\\\)
7.2.3 Diagnóstico inferencial
A análise inferencial consiste em avaliar a qualidade de um modelo em relação à estimação de seus parâmetros (\(\theta\)). Desta forma pode-se pensar no prisma do entendimento do modelo. Este entendimento usualmente passa pela avaliação da aderência do modelo aos dados, atribuindo graus de evidência a hipóteses associadas aos parâmetros modelados ou minimizando algum critério teórico de informação.
Teste para \(\beta_0\)
As hipóteses usuais do teste para \(\beta_0\) são \(H_0: \beta_0 = 0\) vs \(H_1: \beta_0 \ne 0\). Sob \(H_0\) \[\begin{equation} T_0 = \frac{\hat{\beta}_0}{ep(\hat{\beta}_0)} \sim t_{n-2}, \tag{7.19} \end{equation}\] onde \(\hat{\beta}_0\) é dado pela Eq. (7.16) e \[\begin{equation} ep(\hat{\beta}_0) = \sqrt{\hat{\sigma}^2 \left[ \frac{1}{n} + \frac{\bar{x}^2}{S_{xx}} \right] } = \sqrt{\frac{\sum_{i=1}^n (y_i - \hat{y}_i)^2}{n-2} \left[ \frac{1}{n} + \frac{\bar{x}^2}{\sum_{i=1}^n (x_i - \bar{x})^2} \right] }. \tag{7.20} \end{equation}\]
Teste para \(\beta_1\)
O teste para \(\beta_1\) é fundamental na análise de diagnóstico. É com ele que decide-se a respeito da presença ou ausência de relação linear entre \(X\) e \(Y\). As hipóteses usuais do teste para \(\beta_1\) são \(H_0: \beta_1 = 0\) vs \(H_1: \beta_1 \ne 0\). Sob \(H_0\) \[\begin{equation} T_1 = \frac{\hat{\beta}_1}{ep(\hat{\beta}_1)} \sim t_{n-2}, \tag{7.21} \end{equation}\] onde \(\hat{\beta}_1\) é dado pela Eq. (7.15) e \[\begin{equation} ep(\hat{\beta}_1) = \sqrt{\frac{\hat{\sigma}^2}{S_{xx}}} = \sqrt{\frac{\sum_{i=1}^n (y_i - \hat{y}_i)^2 / (n-2)}{\sum_{i=1}^n (x_i - \bar{x})^2} }. \tag{7.22} \end{equation}\]
Teste para \(\beta_1^{RPO}\)
\[\begin{equation} T_1^{RPO} = \frac{\hat{\beta}_1^{RPO}}{ep(\hat{\beta}_1^{RPO})} \sim t_{n-1}, \tag{7.23} \end{equation}\] onde \(\hat{\beta}_1^{RPO}\) é dado pela Eq. (7.18) e \[\begin{equation} ep(\hat{\beta}_1^{RPO}) = \sqrt{\frac{\hat{\sigma}^2}{S_{xx}}} = \sqrt{\frac{\sum_{i=1}^n (y_i - \hat{y}_i)^2 / (n-1)}{\sum_{i=1}^n x_{i}^{2} } }. \tag{7.24} \end{equation}\]
Critérios de informação
AIC - Critério de Informação de Akaike
(Akaike 1974) propõe um critério teórico de informação conforme Eq. (7.25). Quanto menor o valor da medida, melhor (no sentido de mais parcimonioso) é o modelo.
\[\begin{equation} \text{AIC} = -2 l(\hat{\beta}) + 2n_{par}, \tag{7.25} \end{equation}\]
onde \(n_{par}\) indica o número de parâmetros ajustados independentemente no modelo e \(l(\hat{\beta})\) é o logaritmo do valor de máxima verossimilhança do modelo. Para mais detalhes veja ?stats::AIC e ?stats::logLik.
BIC - Critério de Informação Bayesiano
(Schwarz 1978) propõe um critério teórico de informação conforme Eq. (7.26). O autor argumenta31 que tanto o AIC quanto o BIC fornecem uma formulação matemática do princípio da parcimônia na construção de modelos. Porém, para um grande número de observações, os dois critérios podem diferir bastante entre si, sendo mais indicado dar preferência ao BIC.
\[\begin{equation} \text{BIC} = -2 l(\hat{\beta}) + \log(n) n_{par}, \tag{7.26} \end{equation}\]
Exemplo 7.13 Considere o Exemplo 7.11. Pode-se utilizar a função genérica base::summary do R para realizar o diagnóstico dos modelos considerados. O modelo RLS na forma \(y=-19.110+5.915x\) foi atribuído ao objeto fit.
##
## Call:
## lm(formula = gar ~ temp, data = dr)
##
## Residuals:
## Min 1Q Median 3Q Max
## -44.204 -8.261 3.518 10.796 33.540
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -19.1096 22.9195 -0.834 0.411
## temp 5.9147 0.6736 8.780 1.57e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.24 on 28 degrees of freedom
## Multiple R-squared: 0.7336, Adjusted R-squared: 0.7241
## F-statistic: 77.1 on 1 and 28 DF, p-value: 1.565e-09O modelo RPO na forma \(y=5.3590x\) foi atribuído ao objeto fit0.
##
## Call:
## lm(formula = gar ~ temp - 1, data = dr)
##
## Residuals:
## Min 1Q Median 3Q Max
## -45.810 -10.537 3.503 11.979 30.267
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## temp 5.35904 0.09734 55.05 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.14 on 29 degrees of freedom
## Multiple R-squared: 0.9905, Adjusted R-squared: 0.9902
## F-statistic: 3031 on 1 and 29 DF, p-value: < 2.2e-16Os critérios de informação de Akaike (AIC) e Schwartz (BIC) podem ser calculados respectivamente via stats::AIC e stats::BIC. Os dois critérios sugerem que o modelo RPO é mais parcimonioso do que o modelo RLS.
## [1] 263.2736## [1] 267.4772## [1] 262.0094## [1] 264.8118Pode-se ainda realizar uma inspeção visual como complemento da análise inferencial. Para um modelo bem ajustado espera-se que os gráficos de resíduos por valores ajustados (primeira coluna) não indiquem qualquer tipo de tendência ou comportamento modelável. É também esperado para um bom ajuste que o gráfico ‘Normal Q-Q’ apresente os pontos o mais próximo possível da linha pontilhada, que indica normalidade perfeita dos resíduos. A distância de Cook do último gráfico indica possíveis pontos influentes, que podem estar contribuindo de maneira desproporcional ao modelo. Para mais detalhes veja a documentação de ?stats::lm.influence, (R. D. Cook 1977) e (Paula 2013).
Exercício 7.9 Considerando o Exemplo 7.13, verifique:
a. \(\hat{\beta}_0 = -19.110\).
b. \(\hat{\beta}_1 = 5.915\).
c. \(ep(\hat{\beta}_0) = 22.919\).
d. \(ep(\hat{\beta}_1) = 0.674\).
e. \(T_0=-0.83\).
f. \(T_1=8.78\).
g. O p-value associado à hipótese \(H_0:\beta_0=0\) é \(0.41\).
h. O p-value associado à hipótese \(H_0:\beta_1=0\) é \(1.6e-09\).
i. \(\hat{\beta}_1^{RPO} = 5.3590\).
j. \(ep(\hat{\beta}_1^{RPO}) = 0.0973\).
k. \(T_1^{RPO} = 55\).
l. O p-value associado à hipótese \(H_0:\beta_1^{RPO}=0\) é \(<2e-16\).
Exemplo 7.14 Considere novamente os dados do Exemplo 7.1.
a. Faça a regressão de M\(f0 por M\)f0_mhs e vice-versa.
b. O que você percebe após avaliar o summary dos modelos?
Exercício 7.10 Considere novamente os dados do Exemplo 7.2.
a. Realize à mão o teste dos parâmetros do modelo descrito na Equação (7.13) considerando as Equações (7.19), (7.20), (7.21) e (7.22).
b. Realize à mão o teste dos parâmetros do modelo descrito na Equação (7.17) considerando as Equações (7.23) e (7.24).
7.2.3.1 Coeficiente de determinação
Original \[\begin{equation} r^2 = 1 - \frac{SQ_{res}}{SQ_{tot}} = 1 - \frac{\sum_{i=1}^n (y_i-\hat{y}_i)^2}{\sum_{i=1}^n (y_i-\bar{y})^2} \tag{7.27} \end{equation}\]
Ajustado \[\begin{equation} r_{a}^2 = 1 - \frac{SQ_{res}/gl_{res}}{SQ_{tot}/gl_{tot}} = 1 - \frac{(\sum_{i=1}^n (y_i-\hat{y}_i)^2)/(n-p)}{(\sum_{i=1}^n (y_i-\bar{y})^2)/(n-1)} = 1 - (1-r^2) \frac{n-1}{n-p-1} \tag{7.28} \end{equation}\]
7.2.4 Diagnóstico preditivo
A análise preditiva consiste em avaliar a qualidade do modelo em relação à capacidade de predição de novos valores (\(X\)). Desta forma pode-se pensar no prisma da capacidade preditiva do modelo. Tal capacidade pode ser avaliada por diversas medidas, a depender do tipo de variável envolvida.
Variáveis categóricas tipicamente requerem ser agrupadas ou classificadas, e para mensurar a capacidade preditiva utilizam-se medidas tais como acurácia, sensibilidade, sensitividade, etc. Para mais detalhes veja a documentação de ?caret::confusionMatrix e este artigo da Wikipedia.
Modelos para variáveis numéricas tipicamente fazem uso do erro quadrático médio (de previsão), ou EQM, dado pela Eq. (7.29).
\[\begin{equation} \text{EQM} = \frac{\sum_i^{n} (y_i - \hat{y}_i)^2}{n} \tag{7.29} \end{equation}\]
Pode-se escrever uma função para o cálculo do EQM.
Exemplo 7.15 Considere novamente os dados do Exemplo 7.2. Utilizando a ideia de treino-teste, pode-se comparar os modelos quanto aos seus EQMs. Para isso são realizadas \(M=100\) replicações ajustando os modelos RLS e RPO com \(ntrain=20\) observações de treino, confrontadas com as \(ntest=10\) observações de teste restantes. Note que o sorteio da amostra é realizado através da função base::sample associada a sed.seed indexada por i, permitindo a replicação exata dos resultados.
M <- 100 # Número de repetições
ntrain <- 20 # Tamanho da amostra de treino
fit_eqm <- vector(length = M) # Vetor para guardar os M EQMs de fit
fit0_eqm <- vector(length = M) # Vetor para guardar os M EQMs de fit0
for(i in 1:M){
set.seed(i); itrain <- sort(sample(1:nrow(dr), ntrain))
train <- dr[itrain,]
test <- dr[-itrain,]
xtest <- data.frame(temp = dr[-itrain,'temp'])
fit_temp <- lm(gar ~ temp, data = train)
fit0_temp <- lm(gar ~ temp-1, data = train)
fit_eqm[i] <- eqm(test$gar, predict(fit_temp, xtest))
fit0_eqm[i] <- eqm(test$gar, predict(fit0_temp, xtest))
}Pode-se comparar numericamente os EQMs dos modelos através da função genérica base::summary e visualmente por graphics::boxplot.
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 145.9 301.4 362.1 386.7 455.0 815.2## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 137.8 271.6 349.6 365.4 453.6 779.3
Nota-se que os valores de EQM associados ao modelo RPO fit0 são ligeriamente inferiores aos do modelo RLS fit. Desta maneira sugere-se o uso da RPO em detrimento da RLS, ainda que os testes de Student e Mann-Whitney (com seus pressupostos básicos devidamente admitidos) sugiram que não haja diferença significativa entre as médias e medianas dos EQMs dos modelos avaliados.
##
## Welch Two Sample t-test
##
## data: fit_eqm and fit0_eqm
## t = 1.1401, df = 197.99, p-value = 0.2556
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -15.56988 58.24303
## sample estimates:
## mean of x mean of y
## 386.7044 365.3678##
## Wilcoxon rank sum test with continuity correction
##
## data: fit_eqm and fit0_eqm
## W = 5445, p-value = 0.2774
## alternative hypothesis: true location shift is not equal to 0Referências
Qualitatively both our procedure and Akaike’s give “a mathematical formulation of the principle of parsimony in model building.” Quantitatively, since our procedure differs from Akaike’s only in that the dimension is multiplied by \(\frac{1}{2}\log n\), our procedure leans more than Akaike’s towards lower-dimensional models (when there are 8 or more observations). For large numbers of observations the procedures differ markedly from each other. If the assumptions we made in Section 2 are accepted, Akaike’s criterion cannot be asymptotically optimal. This would contradict any proof of its optimality, but no such proof seems to have been published, and the heuristics of Akaike [1] and of Tong [4] do not seem to lead to any such proof.↩︎