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

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 \(60\) zonas da Suécia.

Encontre a equação da reta de regressão ajustada e faça análise dos resíduos.
Alguma transformação precisa ser feita nos dados?

# Lendo a base
base.1 <- read_excel("ex1.xlsx")

# Gráfico de dispersão entre as variáveis
ggplot(base.1, aes(x = X, y = Y)) +
  geom_point(size=2) +
  labs(x="Requerimentos",y="Valor total gasto (milhares SEK)",
       title="Requerimentos segundo \n valor total gasto",
       caption="Fonte: Comitê Sueco de Análise de Prêmio de Risco Seguro Automóvel.")+
  theme(plot.title=element_text(hjust=0.5))

A partir do gráfico acima parece existir relação linear entre as variáveis \(\textit{quantidade de requerimentos}\) e \(\textit{valor gasto}\).

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

# Ajuste do modelo
fit.1 = lm(Y~X,data=base.1)
summary(fit.1)

## 
## Call:
## lm(formula = Y ~ X, data = base.1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -89.150 -25.009   0.858  26.106  82.333 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  18.7090     7.6072   2.459   0.0169 *  
## X             3.5083     0.2953  11.881   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 36.62 on 58 degrees of freedom
## Multiple R-squared:  0.7088, Adjusted R-squared:  0.7038 
## F-statistic: 141.2 on 1 and 58 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 a \(\textit{quantidade de requerimentos}\) e o \(\textit{valor gasto}\). Além disso, \(70,9\%\) da variabilidade do valor gasto no pagamento dos requerimentos é explicado pela regressão. A reta ajustada é

\[ \hat{Y} = 18.710 + 3.508X \]

Agora, uma vez ajustado o modelo, faremos a análise dos resíduos para verificação das suposições do modelo. A seguir calculamos as medidas para iniciar as análises gráficas dos resíduos.

# Residuos
res.1 <- residuals(fit.1)
res.p.1 <- stdres(fit.1)
res.e.1 <- studres(fit.1)

# Valor ajustado
y.hat.1 <- fitted(fit.1)

A figura abaixo apresenta os resíduos versus os valores ajustados. Note que os pontos têm formato de losango e não apresentam aleatoriedade em torno do zero. Isso nos da evidências de que a variância dos erros não é constante.

# Gráficos dos resíduos versus valores ajustados
D.1 <- tibble(y.hat.1,res.1,res.p.1,res.e.1)

ggplot(D.1, aes(x = y.hat.1,y = res.e.1)) +
       geom_point(size=2) +
       labs(x='valor gasto ajustado',y='Resíduos estudentizados',
            title='Resíduos estudentizados versus \n valor gasto ajustado')+
       theme(plot.title=element_text(hjust=0.5)) +
       geom_hline(aes(yintercept = 0), col="red")

Fazendo o teste de Breusch-Pagan para testar a hipótese nula de que a variância dos erros é constante, vemos que a suposição de homocedasticidade é, de fato, rejeitada \((p-valor=0.0046)\).

bptest(fit.1)

## 
##  studentized Breusch-Pagan test
## 
## data:  fit.1
## BP = 8.0208, df = 1, p-value = 0.004624

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

# QQ-plot
ggplot(D.1,aes(sample=res.e.1)) + 
       stat_qq() +
       stat_qq_line(col=4, lwd=1.1) +
       labs(x='quantis teóricos',y='quantis amostrais',title='QQ-plot') +
       theme(plot.title=element_text(hjust=0.5))

Pelo teste de Shapiro-Wilk, o \(p-valor = 0.4434\), não rejeitando a hipótese de normalidade.

shapiro.test(res.e.1)

## 
##  Shapiro-Wilk normality test
## 
## data:  res.e.1
## W = 0.98035, p-value = 0.4434

Vimos que o modelo de regressão ajustado não respeita a suposição de constância da variância. Podemos verificar se alguma transformção nos dados é necessária para contornar esse problema.

A transformação Box-Cox, que consiste em transformar os dados de acordo com a expressão \[ Y'= \begin{cases} \frac{Y^\lambda - 1}{\lambda}, \ \ \ \lambda \neq 0 \\ log(Y), \ \lambda = 0 \end{cases} \]

A função boxcox(), do pacote MASS, calcula o parâmetro \(\lambda\) que maximiza a log-verossimilhança.

boxcox(fit.1, lambda=seq(0, 1, 1/10))

Note que há um intervalo de confiança para o valor de \(\lambda\). Lembre que transformando o dado, todas as análises e concluções serão feitas sob a variável transformada. Assim, consideramos \(\lambda=1/2\) para facilitar a transformação inversa para a escala original dos dados. Para facilicar a interpretação, podemos testar a transformação potência diretamente na variável, isto é, \(Y^\lambda\). Neste caso, ajustaremos o modelo onde a variável resposta é \(\sqrt{Y}\).

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

# Ajuste do modelo
fit.1a = lm(Y^(1/2)~X,data=base.1)
summary(fit.1a)

## 
## Call:
## lm(formula = Y^(1/2) ~ X, data = base.1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0115 -1.8749  0.4001  1.3263  3.8752 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.03191    0.43922   11.46  < 2e-16 ***
## X            0.18407    0.01705   10.80 1.68e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.114 on 58 degrees of freedom
## Multiple R-squared:  0.6677, Adjusted R-squared:  0.662 
## F-statistic: 116.6 on 1 and 58 DF,  p-value: 1.679e-15

\[ \widehat{Y^{1/2}} = 5.032 + 0.184X \]

Agora, para esse novo modelo, faremos a análise dos resíduos para verificação das suposições. A seguir calculamos as medidas para iniciar as análises gráficas dos resíduos.

# Residuos
res.1a <- residuals(fit.1a)
res.p.1a <- stdres(fit.1a)
res.e.1a <- studres(fit.1a)

# Valor ajustado
y.hat.1a <- fitted(fit.1a)

A figura abaixo apresenta os resíduos versus os valores ajustados. Agora, parece que os pontos estão distribuídos de forma aleatória e centrados no zero. Isso nos da evidências de que a variância dos erros é constante.

# Gráficos dos resíduos versus valores ajustados
D.1a <- tibble(y.hat.1a,res.1a,res.p.1a,res.e.1a)

ggplot(D.1a, aes(x = y.hat.1a,y = res.e.1a)) +
       geom_point(size=2) +
       labs(x='valor gasto ajustado',y='Resíduos estudentizados',
            title='Resíduos estudentizados versus \n valor gasto ajustado')+
       theme(plot.title=element_text(hjust=0.5)) +
       geom_hline(aes(yintercept = 0), col="red")

Fazendo o teste de Breusch-Pagan para testar a hipótese nula de que a variância dos erros é constante, vemos que a suposição de homocedasticidade não é rejeitada \((p-valor=0.6567)\).

bptest(fit.1a)

## 
##  studentized Breusch-Pagan test
## 
## data:  fit.1a
## BP = 0.19755, df = 1, p-value = 0.6567

A partir do gráfico quantil-quantil os pontos desviam levemente da reta identidade nos extremos.

# QQ-plot
ggplot(D.1a,aes(sample=res.e.1a)) + 
       stat_qq() +
       stat_qq_line(col=4, lwd=1.1) +
       labs(x='quantis teóricos',y='quantis amostrais',title='QQ-plot') +
       theme(plot.title=element_text(hjust=0.5))

Pelo teste de Shapiro-Wilk rejeita-se a hipótese de normalidade para qualquer nível de significância \(\alpha > 5\%\). Por exemplo, se utilizarmos \(\alpha=1\%\), não rejeitamos a hipótese de que os erros têm distribuição normal, ou seja, aqueles pequenos desvios vistos no gráfico \(QQ-plot\) parecem ser aceitáveis.

shapiro.test(res.e.1a)

## 
##  Shapiro-Wilk normality test
## 
## data:  res.e.1a
## W = 0.9605, p-value = 0.04976

Sendo assim, a transformação utilizada fez com que o novo modelo ajustado respeitasse todos os pressupostos. Nem sempre a transformação será suficiente para contornar os problemas encontrados. Às vezes, o modelo normal pode realmente não ser o mais adequado para representar o comportamento dos dados e, assim, a busca por outros modelos torna-se essencial.

Exemplo 2

Uma pesquisa coletou dados de preços de imóveis no Condado de King, estado de Washington, EUA. Existem vários fatores que influenciam o preço de um imóvel. Nesse exemplo, o interesse é verificar a relação do tamanho do imóvel (em \(m^2\)) e o valor de venda (em dólares). Foram considerados \(21007\) imóveis com mais de \(60m^2\) e menos de \(400m^2\), em \(2014\) e \(2015\).

# Lendo a base
base.2 <- read_excel("ex2.xlsx")

# Gráfico de Dispersão entre as variáveis
ggplot(base.2, aes(x=metro,y=price)) +
  geom_point(size=2,col="darkblue") +
  labs(x='Tamanho (em m²)',y='Preço (em $)',
       title='Preço dos imóveis no Condado de King \n segundo o tamanho (em m²)',
       caption='Fonte: Kaggle') +
  theme(plot.title=element_text(hjust=0.5))

A partir do gráfico de dispersão entre as variávies conseguimos notar relação linear entre as variáveis. De fato, em média, espera-se que quanto maior é o imóvel, mais alto seja seu valor. Então parece razoável ajustar o modelo de regressão linear.

# Ajuste do modelo
fit.2 <- lm(price~metro,data=base.2)
summary(fit.2)

## 
## Call:
## lm(formula = price ~ metro, data = base.2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -624065 -139115  -24494   97661 2355384 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 38373.03    4302.54   8.919   <2e-16 ***
## metro        2542.51      21.49 118.303   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 223800 on 21005 degrees of freedom
## Multiple R-squared:  0.3999, Adjusted R-squared:  0.3998 
## F-statistic: 1.4e+04 on 1 and 21005 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 o \(\textit{tamanho}\) e o \(\textit{preço}\), apesar de apenas \(40,0\%\) da variabilidade do preço de venda ser explicado pela variabilidade do tamanho do imóvel. A reta ajustada é: \[\hat{Y} = 38373.03 + 2542.51 X\]

Agora, uma vez ajustado o modelo, iniciaremos as análises do resíduos.

# Residuos
res.2 <- residuals(fit.2)
res.p.2 <- stdres(fit.2)
res.e.2 <- studres(fit.2)

# Valor ajustado
y.hat.2 <- fitted(fit.2)

A figura abaixo apresenta os resíduos versus os valores ajustados. Observe que o grafico tem formato de funil, indicando que a variância não é constante. Quando variamos o nível de \(\hat{Y}\), a variância aumenta. Portanto, a suposição de homocedasticidade parece não ser válida.

# Gráficos dos resíduos
D.2 = tibble(y.hat.2,res.2,res.p.2,res.e.2)

ggplot(D.2, aes(x = y.hat.2,y = res.e.2)) +
       geom_point(size=2) +
       labs(x='Preço dos imóveis ajustado',y='Resíduos estudentizados',
            title='Resíduos estudentizados versus \n preço dos imóveis ajustado')+
       theme(plot.title=element_text(hjust=0.5)) +
       geom_hline(aes(yintercept = 0), col="red")

Fazendo o teste de Breusch-Pagan para testar a hipótese nula de que a variância dos erros é constante, vemos que a suposição de homocedasticidade é, de fato, rejeitada \((p-valor \approx 0)\).

bptest(fit.2)

## 
##  studentized Breusch-Pagan test
## 
## data:  fit.2
## BP = 1419.7, df = 1, p-value < 2.2e-16

Por fim, a figura abaixo apresenta o gráfico quantil-quantil. Observe que a distribuição parece ser assimétrica.

# Verificação de normalidade
g.3 <- ggplot(D.2,aes(sample=res.e.2)) + 
             stat_qq() +
             stat_qq_line(col=4, lwd=1.1) +
             labs(x='quantis teóricos',y='quantis amostrais',title='QQ-plot') +
             theme(plot.title=element_text(hjust=0.5))

g.4 <- ggplot(D.2,aes(x=res.e.2)) + 
             geom_histogram(color="white", fill="gray40") +
             labs(x='Resíduos estudentizados',y='Frequência',title='Histograma dos resíduos')+
             theme(plot.title=element_text(hjust=0.5))

grid.arrange(g.3,g.4,ncol=2)

O teste de Shapiro-Wilk tem limitação quanto à quantidade de dados. Quando fornecemos mais dados, as chances de a hipótese nula ser rejeitada se torna maior. Mesmo desvios muito pequenos da normalidade podem ser detectados, levando à rejeição da hipótese nula, embora para fins práticos os dados sejam “normais”. Por isso, a função shapiro.test() restringe o limite de observações avaliadas em 5000. Assim, como temos uma amostra de tamanho 21007, utilizaremos o teste de Kolmogorov-Smirnov para testar a hipótese de normalidade. A função utilizada, netse caso, é lillie.test().

lillie.test(res.e.2)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  res.e.2
## D = 0.082546, p-value < 2.2e-16

O modelo de regressão normal parece não ser o modelo ideal para representar a relação entre as variáveis avaliadas. Os pressupostos de homocedasticidade e normalidade são violados. Visivelmente, o distribuição é assimétrica e, nesse caso, a transformação log() pode ser utilizada. O modelo de regressão normal que tem como variável resposta o \(log(Y)\) é conhecido como modelo log-normal.

Avaliando a transformação a partir da função boxcox(), parece razoável usar o log(). Neste caso, ajustaremos o modelo onde a variável resposta é \(log(Y)\).

boxcox(fit.2, lambda=seq(-0.1,0.1, 1/10))

Agora, precisamos ajustar o modelo com a variável transformada e verificar seus pressupostos.

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

# Ajuste do modelo
fit.2a = lm(log(price)~metro,data=base.2)
summary(fit.2a)

## 
## Call:
## lm(formula = log(price) ~ metro, data = base.2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.2492 -0.2838  0.0142  0.2591  1.2595 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.219e+01  7.204e-03  1692.4   <2e-16 ***
## metro       4.450e-03  3.598e-05   123.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3747 on 21005 degrees of freedom
## Multiple R-squared:  0.4213, Adjusted R-squared:  0.4213 
## F-statistic: 1.529e+04 on 1 and 21005 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 as variáveis

Agora, para esse novo modelo, faremos a análise dos resíduos para verificação das suposições. A seguir calculamos as medidas para iniciar as análises gráficas dos resíduos.

# Residuos
res.2a <- residuals(fit.2a)
res.p.2a <- stdres(fit.2a)
res.e.2a <- studres(fit.2a)

# Valor ajustado
y.hat.2a <- fitted(fit.2a)

A figura abaixo apresenta os resíduos versus os valores ajustados. Agora, parece que os pontos estão distribuídos de forma aleatória em torno de zero, com poucos pontos fora do padrão na parte inferior esquerda. Isso nos da evidências de que a variância dos erros é constante.

# Gráficos dos resíduos versus valores ajustados
D.2a <- tibble(y.hat.2a,res.2a,res.p.2a,res.e.2a)

ggplot(D.2a, aes(x = y.hat.2a,y = res.e.2a)) +
       geom_point(size=2) +
       labs(x='preço do imóvel ajustado',y='Resíduos estudentizados',
            title='Resíduos estudentizados versus \n preço do imóvel ajustado')+
       theme(plot.title=element_text(hjust=0.5)) +
       geom_hline(aes(yintercept = 0), col="red")

Fazendo o teste de Breusch-Pagan para testar a hipótese nula de que a variância dos erros é constante, vemos que a suposição de homocedasticidade não é rejeitada \((p-valor=0.3456)\).

bptest(fit.2a)

## 
##  studentized Breusch-Pagan test
## 
## data:  fit.2a
## BP = 0.88968, df = 1, p-value = 0.3456

A partir do gráfico quantil-quantil é possível visualizar que os pontos desviam da reta identidade nos extremos.

# QQ-plot
g.1 <- ggplot(D.2a,aes(sample=res.e.2a)) + 
            stat_qq() +
            stat_qq_line(col=4, lwd=1.1) +
            labs(x='quantis teóricos',y='quantis amostrais',title='QQ-plot') +
            theme(plot.title=element_text(hjust=0.5))
g.2 <- ggplot(D.2,aes(x=res.e.2a)) + 
            geom_histogram(color="white", fill="gray40") +
            labs(x='Resíduos estudentizados',y='Frequência',title='Histograma dos resíduos')+
            theme(plot.title=element_text(hjust=0.5))

grid.arrange(g.1,g.2,ncol=2)

Pelo teste de Lilliefors rejeita-se a hipótese de normalidade. Entretanto, pelo histograma dos resíduos, parece que a distribuição tem o comportamento próximo da normal.

lillie.test(res.e.2a)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  res.e.2a
## D = 0.027299, p-value < 2.2e-16

A transformação utilizada fez com que o novo modelo ajustado respeitasse a hipótese de homocedasticidade. Apesar do teste de Lilliefors rejeitar a hipótese de normalidade, parece razoavel considerar que os erros seguem distribuição aproximadamente normal, a partir da análise gráfica. Nem sempre a transformação será suficiente para contornar os problemas encontrados. Neste caso, o modelo normal pode realmente não ser o mais adequado para representar o comportamento dos dados.

Modelos Lineares I

Modelos Lineares I

Rafael Erbisti

Rafael Erbisti

Beatriz Mello

Laboratório 04

Laboratório 04

Exemplo 1

Exemplo 2