class: center, middle, inverse, title-slide # Diagnóstico dos Modelos de Regressão ## Como verificar pressupostos dos modelos de Regressão ### Steven Dutt Ross --- # Construindo um modelo de regressão ## Vamos construir um modelo para prever as vendas com base nas despesas com publicidade em mídias do youtube. ```r modelo <- lm(vendas ~ youtube, data = marketing) modelo ``` ``` ## ## Call: ## lm(formula = vendas ~ youtube, data = marketing) ## ## Coefficients: ## (Intercept) youtube ## 8.43911 0.04754 ``` A equação do nosso modelo de regressão é: `$$y = 8,43 + 0,07*x$$`, isto é, `$$vendas = 8,43 + 0,047*youtube$$`. --- # O que são as suposições do modelo ? O **modelo de regressão linear faz várias suposições** sobre os dados disponíveis. vamos descreve quais são as **suposições do modelo** e fornece meios para o diagnóstico desses **pressupostos da regressão** na linguagem de programação R. Depois de realizar uma regressão, **você deve sempre verificar** se o modelo funciona bem para os dados disponíveis. Ao desenvolver modelos de regressão, buscamos inspecionar a **significância dos coeficientes da regressão**, bem como o R2 que nos **indica quão bem o modelo de regressão linear se ajusta aos dados**. Todavia, não podemos ficar somente nesta análise. É necessário realizar etapas adicionais para avaliar o quão bem o modelo se ajusta aos dados e **ter certeza que não violamos nenhum pressuposto da regressão**. --- # Por que verificar pressupostos em modelos de regressão? O modelo de **regressão linear** faz a suposição de que a relação entre os preditores (x) e a variável resposta é **linear**. Isso pode não ser verdade. O relacionamento pode ser polinomial ou logarítmico. Os dados também podem conter **outliers**. Isso pode afetar o resultado da regressão. Portanto, devemos gerar um diagnóstico do modelo de regressão construído para detectar possíveis problemas e verificar se as suposições (pressupostos) feitas pelo modelo de regressão linear são atendidas ou não. Para fazer isso, geralmente examinamos a **distribuição de erros residuais**, que podem fornecer outras informações sobre seus dados. --- # Valores ajustados e residuais Os valores ajustados (ou previstos) são os valores de y que você esperaria para os valores de x fornecidos de acordo com o modelo de regressão construído (ou visualmente, a linha de regressão reta de melhor ajuste). No nosso modelo, para um gasto em publicidade no youtube, o valor das vendas ajustado (previsto) seria de vendas = 8,44 + 0,0048 * youtube. A partir do gráfico de dispersão abaixo, pode ser visto que nem todos os pontos de dados caem exatamente na linha de regressão estimada. Isso significa que, para um determinado orçamento de publicidade do YouTube, os valores de venda observados podem ser diferentes dos valores de venda previstos. A diferença é chamada de erros residuais, representados por linhas vermelhas verticais. <img src="C:/Users/Steven/Documents/GitHub/Regressao_Diagnostico/linear-regression.png" width="152" style="display: block; margin: auto;" /> --- Podemos aumentar os dados para adicionar valores ajustados e resíduos usando a função *augment()* do *[broom package]*. Vamos guardar a saída em um onjeto chamado *model.diag.metrics* porque ela contém várias métricas úteis para diagnósticos de regressão. ```r library(broom) ``` ``` ## Warning: package 'broom' was built under R version 3.5.2 ``` ```r model.diag.metrics <- augment(modelo) head(model.diag.metrics) ``` ``` ## # A tibble: 6 x 9 ## vendas youtube .fitted .se.fit .resid .hat .sigma .cooksd .std.resid ## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 26.5 276. 21.6 0.385 4.96 0.00970 3.90 0.00794 1.27 ## 2 12.5 53.4 11.0 0.431 1.50 0.0122 3.92 0.000920 0.387 ## 3 11.2 20.6 9.42 0.502 1.74 0.0165 3.92 0.00169 0.449 ## 4 22.2 182. 17.1 0.277 5.12 0.00501 3.90 0.00434 1.31 ## 5 15.5 217. 18.8 0.297 -3.27 0.00578 3.91 0.00205 -0.839 ## 6 8.64 10.4 8.94 0.525 -0.295 0.0180 3.92 0.0000534 -0.0762 ``` --- Entre as colunas da tabela, existem: * youtube: o orçamento de publicidade do youtube investido * vendas: os valores de venda observados * .fitted: os valores de venda ajustados * .resid: os erros residuais No gráfico a seguir temos o erro residual (na cor vermelha) entre os valores observados e a reta de regressão ajustada. Cada segmento vertical vermelho representa um erro residual entre um valor de venda observado e o correspondente valor previsto (ou seja, ajustado). --- ``` ## Warning: package 'ggplot2' was built under R version 3.5.3 ``` <img src="index_files/figure-html/unnamed-chunk-4-1.png" width="672" /> --- # Pressupostos dos modelos de regressão Os modelos de regressão linear fazem várias suposições sobre os dados, como: 1. **Linearidade dos dados**: A relação entre o preditor (x) e a variável resposta (y) é linear. 2. **Normalidade dos resíduos**: Os erros residuais são considerados normalmente distribuídos (seguem a distribuição normal). 3. **Homogeneidade da variância**. Presume-se que os resíduos tenham uma variância constante (homocedasticidade) 4. **Independência** entre cada um dos termos de erro dos resíduos. Devemos verificar se essas suposições são verdadeiras ou não toda vez que . Problemas potenciais incluem: 1. Não linearidade do resultado - relações preditoras 2. Heterocedasticidade: Variância não constante dos termos de erro. 3. Presença de valores influentes nos dados. Outliers: valores extremos na variável resposta (y) Pontos de alta alavancagem: valores extremos na variável preditores (x) Todas essas suposições e problemas potenciais podem ser verificados produzindo alguns gráficos de diagnóstico dos erros residuais. --- # Diagnóstico dos Modelos de Regressão ## Gráficos de diagnóstico Os gráficos de diagnósticos de regressão podem ser criados facilmente com o R. ```r par(mfrow = c(2, 2)) plot(modelo) # ou par(mfrow = c(2, 2)) library(ggfortify) autoplot(modelo) ``` --- ``` ## Warning: package 'ggfortify' was built under R version 3.5.3 ``` <img src="index_files/figure-html/unnamed-chunk-6-1.png" width="672" /> --- Os gráficos de diagnóstico mostram os resíduos de quatro maneiras diferentes: 1. **Residuals vs Fitted**: Usado para verificar as suposições de relacionamento linear. Uma linha horizontal, sem padrões distintos, é uma indicação de um relacionamento linear, o que é bom. 2. **Normal Q-Q**: Usado para examinar se os resíduos são normalmente distribuídos. É bom se os pontos residuais seguirem a linha reta. 3. **Scale-Location**: Usado para verificar a homogeneidade da variância dos resíduos (homocedasticidade). A linha horizontal com pontos igualmente espalhados é uma boa indicação de homocedasticidade. Este não é o caso em nosso exemplo, onde temos um problema de heteroscedasticidade. 4. **Residuals vs Leverage**: Usado para identificar casos influentes, ou seja, valores extremos que podem influenciar os resultados da regressão quando incluídos ou excluídos da análise. --- # Linearidade dos dados O pressuposto de linearidade pode ser verificada inspecionando o gráfico *Residuals vs Fitted*: ```r plot(modelo, 1) ``` <img src="index_files/figure-html/unnamed-chunk-7-1.png" width="672" /> Idealmente, o gráfico residual não mostrará nenhum padrão ajustado. Ou seja, a linha vermelha deve ser aproximadamente horizontal em zero. A presença de um padrão pode indicar um problema com algum aspecto do modelo linear. No nosso exemplo, **não há padrão no gráfico dos resíduos**. Isso sugere que podemos assumir uma relação linear entre os preditores e a variável resposta. Note que, se o gráfico dos resíduos indicar uma relação NÃO linear nos dados, então uma abordagem simples é usar **transformações não-lineares**, como log(x), sqrt(x) e x^2, no modelo de regressão. --- ### Homocedasticidade (Homogeneidade de variância) ```r plot(modelo, 3) ``` <img src="index_files/figure-html/unnamed-chunk-8-1.png" width="672" /> --- # Homocedasticidade ### (Homogeneidade de variância) Nesse gráfico seria bom se podemos ver uma linha horizontal com pontos igualmente espalhados. No nosso modelo, esse não é o caso. Pode-se observar que a variabilidade (variância) dos pontos residuais aumenta com o valor da variável de resposta, sugerindo variâncias não constantes nos erros dos resíduos (ou heteroscedasticidade). --- Uma solução possível para reduzir o problema da heterocedasticidade é usar uma transformação de log ou raiz quadrada da variável de resposta (y). ```r modelo2 <- lm(log(vendas) ~ youtube, data = marketing) plot(modelo2, 3) ``` <img src="index_files/figure-html/unnamed-chunk-9-1.png" width="672" /> --- ### Normalidade dos resíduos O gráfico *QQNorm* pode ser usado para verificar visualmente o pressuposto de normalidade. O gráfico dos resíduos deve seguir aproximadamente a linha reta. Todos os pontos caem aproximadamente ao longo desta linha de referência, então podemos assumir a normalidade. --- ```r plot(modelo, 2) ``` <img src="index_files/figure-html/unnamed-chunk-10-1.png" width="672" /> --- # Normalidade dos resíduos  --- # Outliers Um *outlier* é um ponto que tem um valor extremo. A presença de *outliers* pode afetar a interpretação do modelo, pois aumenta o erro. Os *outliers* podem ser identificados examinando os resíduos padronizados. Resíduos padronizados podem ser interpretados como o número de erros-padrão longe da reta de regressão. Observações cujos resíduos padronizados são maiores que 3 em valor absoluto são possíveis outliers. --- ```r plot(modelo, 5) ``` <img src="index_files/figure-html/unnamed-chunk-12-1.png" width="672" /> O gráfico acima destaca os 3 pontos mais extremos (linhas 26, 36 e 179 do banco de dados), com um resíduo padronizado abaixo de -2. No entanto, não há valores discrepantes que excedam 3 desvios padrão, o que é bom. --- background-image: url(https://raw.githubusercontent.com/DATAUNIRIO/Regressao_Diagnostico/master/checkpoint.png) --- Descobri a pouco tempo um novo pacote do R que ajuda na construção do diagnóstico. O pacote gvlma, realiza uma validação global de pressupostos do modelo linear, bem como avaliações separadas de assimetria, curtose e heterocedasticidade. ```r # Global test of model assumptions library(gvlma) ``` ``` ## Warning: package 'gvlma' was built under R version 3.5.2 ``` ```r gvmodel <- gvlma(modelo) summary(gvmodel) ``` ``` ## ## Call: ## lm(formula = vendas ~ youtube, data = marketing) ## ## Residuals: ## Min 1Q Median 3Q Max ## -10.0632 -2.3454 -0.2295 2.4805 8.6548 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 8.439112 0.549412 15.36 <2e-16 *** ## youtube 0.047537 0.002691 17.67 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 3.91 on 198 degrees of freedom ## Multiple R-squared: 0.6119, Adjusted R-squared: 0.6099 ## F-statistic: 312.1 on 1 and 198 DF, p-value: < 2.2e-16 ## ## ## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS ## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM: ## Level of Significance = 0.05 ## ## Call: ## gvlma(x = modelo) ## ## Value p-value Decision ## Global Stat 5.3240 0.25564 Assumptions acceptable. ## Skewness 0.2619 0.60885 Assumptions acceptable. ## Kurtosis 0.4070 0.52352 Assumptions acceptable. ## Link Function 3.6906 0.05472 Assumptions acceptable. ## Heteroscedasticity 0.9646 0.32602 Assumptions acceptable. ``` --- background-image: url(https://raw.githubusercontent.com/DATAUNIRIO/Regressao_Diagnostico/master/exemplo_gvlma.png) --- class: center, middle # Análises Adicionais --- ## Multicolinearidade Multicolinearidade consiste em um problema comum em regressões, no qual as variáveis independentes possuem relações lineares. Em outras palavras, as variáveis explicativas são correlacionadas. O índício mais claro da existência da multicolinearidade é quando o R² é alto, mas nenhum dos coeficientes da regressão é estatisticamente significativo segundo a estatística t convencional. ### Detectando a multicolinearidade A função vif() do pacote *car* pode ser utilizada para detectar a multicolinearidade em um modelo de regressão. Valores acima de 5 são considerados multicolineariedade fraca. Valores acima de 10 são considerados multicolineariedade forte. ```r modelo3 <- lm(vendas ~ youtube + facebook + jornal, data = marketing) car::vif(modelo3) ``` ``` ## youtube facebook jornal ## 1.004611 1.144952 1.145187 ``` --- ## Multicolinearidade No modelo a variável com maior valor foi o Jornal com VIF = 1,145187 indicando que não temos problema de multicolineariedade. Existem várias abordagens possíveis quando estamos enfrentando o problema de Multicolinearidade. A mais simples é retirar do modelo uma das variáveis independentes. --- background-image: url(https://raw.githubusercontent.com/DATAUNIRIO/Regressao_Diagnostico/master/checkpoint.png) --- # Uma outra abordagem para verificar os pressupostos Existe uma outra abordagem além da análise gráfica para verificar os pressupostos:os **testes de hipóteses** Por exemplo, para verificar a Homocedasticidade existem testes que são úteis para estabelecer a presença ou ausência de heterocedasticidade. Vamos esstudar um deles. O teste de Breush-Pagan. Nesse segmento, vamos utilizar a base de dados *cars*. --- # O teste Breusch-Pagan H0: não há heterocedasticidade. H1: há heterocedasticidade. ```r lmMod <- lm(dist ~ speed, data=cars) library(lmtest) bptest(lmMod, studentize = TRUE) ``` ``` ## ## studentized Breusch-Pagan test ## ## data: lmMod ## BP = 3.2149, df = 1, p-value = 0.07297 ``` o teste tem um p-valor menor que um nível de significância de 0,05, portanto podemos rejeitar a hipótese nula de que a variância dos resíduos é constante (Homocedastica) e inferir que a heterocedasticidade está de fato presente. --- # Como corrigir? Transformação da variável resposta: a transformação Box-Cox. ### Transformação Box-Cox A transformação Box-cox é uma transformação matemática da variável para aproximá-la de uma distribuição normal. Muitas vezes, fazer uma transformação box-cox da variável Y resolve o problema. ```r distBCMod <- caret::BoxCoxTrans(cars$dist) cars <- cbind(cars, dist_new=predict(distBCMod, cars$dist)) #head(cars) lmMod_bc <- lm(dist_new ~ speed, data=cars) bptest(lmMod_bc) ``` ``` ## ## studentized Breusch-Pagan test ## ## data: lmMod_bc ## BP = 0.011192, df = 1, p-value = 0.9157 ``` Com um p-valor de 0,91, não rejeitamos a hipótese nula (que a variância dos resíduos é constante e, portanto, inferimos que os resíduos são homocedásticos. --- # Teste de Normalidade H0: os dados seguem uma distribuição normal. H1: os dados NÃO seguem uma distribuição normal. ```r residuos<-residuals(lmMod_bc) shapiro.test(residuos) ``` ``` ## ## Shapiro-Wilk normality test ## ## data: residuos ## W = 0.97332, p-value = 0.3143 ``` o teste tem um p-valor (0,3143) maior que o nível de significância de 0,05, portanto não rejeitamos a hipótese nula (normalidade dos dados). --- # O teste de Bonferroni para outlier ```r library(car) ``` ``` ## Warning: package 'car' was built under R version 3.5.3 ``` ``` ## Loading required package: carData ``` ```r outlierTest(lmMod_bc) ``` ``` ## No Studentized residuals with Bonferonni p < 0.05 ## Largest |rstudent|: ## rstudent unadjusted p-value Bonferonni p ## 23 3.148775 0.0028478 0.14239 ``` --- class: center, middle # Análises Avançadas --- # Análises Avançadas Uma violação comum das suposições na regressão MQO é a suposição de homocedasticidade. Essa suposição exige que o termo de erro tenha variância constante em todos os valores das variáveis independentes. Quando essa suposição falha, os erros padrão de nossas estimativas de regressão MQO são inconsistentes. Até agora vimos duas abordagens para identificar a heterocedasticidade (gráfico e teste) e duas formas para corrigi-la (trasnformações simples e transformações de box-cox). Mas esses não são os unicos métodos. vou apresentar mais dois métodos: 1. Mínimos Quadrados Generalizados. 2. Erro padrão consistente a heterocedasticidade (Heteroskedasticity-consistent standard errors). --- # Mínimos Quadrados Generalizados Até agora, temos lidado com a heterocedasticidade sob o ponto de vista do MQO. Todavia, se tivessemos conhecimento da matriz de variância-covariância do termo de erro, então poderíamos fazer um modelo heteroscedástico se tornar um modelo homocedástico. Como não temos, vamos utilizar a estimativa da matriz. ```r # Feasible Generalized Least Square GLS lmMod <- lm(dist ~ speed, data=cars) fgls <- lm(dist ~ speed, data=cars, weights = 1/lmMod$fitted.values^2) summary(fgls) ``` ``` ## ## Call: ## lm(formula = dist ~ speed, data = cars, weights = 1/lmMod$fitted.values^2) ## ## Weighted Residuals: ## Min 1Q Median 3Q Max ## -2.08780 -0.26997 -0.03402 0.25451 2.23779 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -6.4508 1.6875 -3.823 0.00038 *** ## speed 3.0780 0.3304 9.317 2.43e-12 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.5932 on 48 degrees of freedom ## Multiple R-squared: 0.6439, Adjusted R-squared: 0.6365 ## F-statistic: 86.81 on 1 and 48 DF, p-value: 2.434e-12 ``` --- # Erros-padrão consistentes a heterocedasticidade Podemos calcular erros padrão consistentes de heterocedasticidade de forma fácil. ```r # Método original library(sandwich) coeftest(lmMod, vcov = vcovHC(lmMod)) ``` ``` ## ## t test of coefficients: ## ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -17.57909 5.93180 -2.9635 0.004722 ** ## speed 3.93241 0.42754 9.1978 3.636e-12 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` ```r # Método mais simples (Robust standard errors for regression models) library(sjstats) ``` ``` ## Warning: package 'sjstats' was built under R version 3.5.3 ``` ``` ## ## Attaching package: 'sjstats' ``` ``` ## The following object is masked from 'package:broom': ## ## bootstrap ``` ```r robust(lmMod) ``` ``` ## term estimate std.error statistic p.value ## 1 (Intercept) -17.579095 5.9318033 -2.963533 4.722042e-03 ## 2 speed 3.932409 0.4275372 9.197816 3.635819e-12 ``` --- Fontes: [STHDA](http://www.sthda.com/english/articles/39-regression-model-diagnostics/161-linear-regression-assumptions-and-diagnostics-in-r-essentials/) [how-to-detect-heteroscedasticity](https://www.r-bloggers.com/how-to-detect-heteroscedasticity-and-rectify-it/) [Rpubs](https://rpubs.com/cbw1243/gls) [Robust](https://www.brodrigues.co/blog/2018-07-08-rob_stderr/) [FGLS](http://www3.grips.ac.jp/~yamanota/Lecture_Note_10_GLS_WLS_FGLS.pdf)