# Pacotes utilizados
library(readxl)
library(tidyverse)
library(lmtest)
library(MASS)
library(gridExtra)
library(forecast)

Exemplo 1

O Comitê Sueco de Análise de Prêmio de Risco no Seguro de Automóveis avaliou a relação entre a quantidade de requerimentos de crédito e o valor gasto para o pagamento de todos os requerimentos (em milhares de coroas suecas) em 54 zonas da Suécia.

  1. Encontre a equação da reta de regressão ajustada. É razoável pensar na hipótese de que quando a quantidade de requerimentos for igual a zero, o valor gasto para o pagamento também será?

  2. Avalie o modelo de regressão sem intercepto.

A partir do gráfico acima parece existir relação linear entre as variáveis \(\textit{quantidade de requerimentos}\) e \(\textit{valor gasto}\). Além disso, note que é intuitivo pensarmos que quando não há requerimentos, o valor gasto será zero. Se traçarmos uma reta para representar a associação entre as variaveis, muito provavelmente ela passará no zero.

Ajustando o modelo de regressão, temos o seguinte resultado.

## 
## Call:
## lm(formula = Y ~ X, data = base.1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -71.693 -18.530   0.818  20.131  62.422 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.6109     7.2548   0.773    0.443    
## X             4.5949     0.3519  13.056   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 31.57 on 52 degrees of freedom
## Multiple R-squared:  0.7662, Adjusted R-squared:  0.7617 
## F-statistic: 170.4 on 1 and 52 DF,  p-value: < 2.2e-16

Neste modelo, a reta de regressão ajustada é \[ \hat{Y} = 5.611 + 4.595X \] Pelo teste \(t\), não rejeitamos a hipótese de que \(\beta_0=0\). Portanto, como esperado, temos evidências de que a reta de regressão passa pela origem. Sendo assim, é preferível que o modelo sem o intercepto seja estimado. Note também que \(MQE=31.57\).

## 
## Call:
## lm(formula = Y ~ X - 1, data = base.1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -71.127 -15.634   1.072  21.351  62.330 
## 
## Coefficients:
##   Estimate Std. Error t value Pr(>|t|)    
## X   4.8142     0.2076   23.19   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 31.45 on 53 degrees of freedom
## Multiple R-squared:  0.9103, Adjusted R-squared:  0.9086 
## F-statistic: 537.7 on 1 and 53 DF,  p-value: < 2.2e-16

No modelo sem interpto, a reta de regressão ajustada é \[ \hat{Y} = 4.814X \]

Note que o valor do \(MQE\) agora é \(31.45\), inferior ao modelo com intercepto. Portanto, neste caso, é preferível utilizar o modelo sem o intercepto.

Exemplo 2

Uma pesquisa coletou dados de consumo de cerveja em São Paulo, em uma área universitária, onde acontecem algumas festas. O estudo avaliou a relação entre a temperatura máxima (em ºC) e o consumo de cerveja (em mil litros) diário em 2015.

  1. Encontre a reta de regressão linear estimada.

  2. Faça a análise dos resíduos.

A partir do gráfico de dispersão entre as variávies conseguimos notar uma relação linear entre as variáveis. Então parece razoável ajustar o modelo de regressão linear para avaliar a relação entre a \(\textit{temperatura máxima}\) e o \(\textit{consumo de cerveja}\).

## 
## Call:
## lm(formula = Consumo ~ Temp_Max, data = base.2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.9116 -2.8451 -0.3342  2.3929  8.6191 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7.97494    1.10459    7.22 3.07e-12 ***
## Temp_Max     0.65485    0.04097   15.98  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.375 on 363 degrees of freedom
## Multiple R-squared:  0.413,  Adjusted R-squared:  0.4114 
## F-statistic: 255.4 on 1 and 363 DF,  p-value: < 2.2e-16

Avaliando o ajuste, os coeficientes \(\beta_0\) e \(\beta_1\) são diferentes de zero e pelo teste \(F\) há associação linear significativa entre \(\textit{temperatura máxima}\) e o \(\textit{consumo de cerveja}\), apesar de apenas \(41,3\%\) da variabilidade do consumo de cerveja ser explicado pela regressão. A reta ajustada é: \[\hat{Y} = 7.975 + 0.655 X\]

Agora, uma vez ajustado o modelo, iniciaremos as análises do resíduos. As funções residuals(), stdres() e studres() calculam os resíduos usuais, os resíduos padronizados e os resíduos estudentizados, respectivamente. Todas essas funções pertencem ao pacote MASS.

Os valores ajustados podem ser encontrados a partir da função fitted(). Portanto, abaixo, calculamos essas medidas para iniciar as análises gráficas dos resíduos.

As figuras abaixo apresentam os resíduos versus os valores ajustados e a covariável. Como visto em aula, no caso da regressão simples, estes dois gráficos são equivalentes. De maneira geral, todos os pontos variam, de forma aleatória, em torno do zero. Entretanto, se avaliarmos as figuras, há resíduos que estão no extremo do gráfico (em termos do eixo \(Y\)) que desviam um pouco da constância da variabilidade. De qualquer forma, parece não haver quebra da suposição de homogeneidade da variância.

A partir dessas figuras também é possível verificar que não há quebra da suposição de linearidade entre as variáveis analisadas.

A função bptest(), do pacote lmtest, realiza o teste de Breusch-Pagan (BP). As hipóteses testadas são: \[ H_0: \textrm{variância dos erros é constante}\\ H_1: \textrm{variância dos erros não é constante} \]

## 
##  studentized Breusch-Pagan test
## 
## data:  fit.2
## BP = 3.7723, df = 1, p-value = 0.05211

Como comentado anteriormente, parece não haver padrão no gráfico dos resíduos acima, indicando de que a hipótese de homocedasticidade é satisfeita. O teste de BP apresentou \(p-valor=0.0521\), não rejeitando a hipótese de homocedasticidade aos níveis \(\alpha<5\%\). Provavelmente, aqueles poucos pontos extremos estão fazendo com que o teste não rejeite a hipótese nula para níveis maiores de \(\alpha\). Sendo assim, o teste também nos fornece indícios de que a suposição de homocedasticidade é satisfeita.

Como os dados são coletados diariamente, faz sentindo avaliarmos se existe autocorrelação nos erros. O dado é temporal e espera-se que o consumo de cerveja seja maior nos fins de semana e, além disso, o valor observado num dia devem ter relação com o valor observado no dia anterio. Ao ajustar um modelo de regressão aos dados de série temporal, é comum encontrar autocorrelação nos resíduos. Nesse caso, o modelo estimado viola o pressuposto de não haver autocorrelação nos erros, e as previsões podem ser ineficientes. Portanto, devemos sempre olhar para um gráfico ACF dos resíduos.

A função checkresiduals(), do pacote forecast, apresenta o gráfico do resíduos versos a ordem de coleta, o gráfico ACF e o histograma dos resíduos. Além disso, essa função ainda tem a opção de apresentar dois testes de autocorrelação (Ljung-Box ou Breusch-Godfrey). Avaliando a figura da ACF é visíveil que existe dependência na ordem de coleta (temporal).

O teste de Durbin-Watson avalia a relação \[ \epsilon_t=\rho \epsilon_{t-1} + \nu_t \] com \(t\) sendo o índice que indica a ordem da coleta (geralmente é o tempo). Então, as hipóteses do teste são: \[ H_0: \rho=0\\ H_1: \rho \neq 0 \] Observe que se \(\rho=0\), então \(\epsilon_t=\nu_t\), i.e, os erros não são autocorrelacionados. A função dwtest(), do pacote lmtest realiza o teste de Durbin-Watson. Para este conjunto de dados, \(p-valor=0.0013\), indicando que rejeitamos a hipótese de que \(\rho=0\), ou seja, há evidências de que o erros são autocorrelacionados (também havíamos chegado nessa conclusão a partir da análise gráfica).

## 
##  Durbin-Watson test
## 
## data:  fit.2
## DW = 1.6883, p-value = 0.001253
## alternative hypothesis: true autocorrelation is greater than 0

Por fim, a figura abaixo apresenta o gráfico quantil-quantil. Observe que os dados parecem não ter distribuição normal. Os pontos iniciais desviam muito da reta identidade.

A função shapito.test(), do pacote stats, realiza o teste de normalidade. As hipóteses testadas são: \[ H_0: \textrm{os erros são normais}\\ H_1: \textrm{os erros não são normais} \] Observe que o \(p-valor \approx 0\), rejeitando a hipótese de normalidade (como esperado).

## 
##  Shapiro-Wilk normality test
## 
## data:  res.e
## W = 0.97723, p-value = 1.625e-05

Vimos que o modelo de regressão normal parece não ser o modelo ideal para representar a relação entre \(\textit{temperatura máxima}\) e o \(\textit{consumo de cerveja}\) numa região universitária de São Paulo, em 2015. Como alguns presupostos do modelo foram violados, precisamos encontrar alternativas para modelar esses dados. As vezes,uma simples transformação na variável resposta é suficiente para que a suposição de normalidade passe a ser válida. Entretanto, os dados são claramente correlacionados no tempo e, talvez, um modelo de série temporal seja mais adequado.