Prática no R! Roteiro 5 - Regressão Linear
Elaborado por Luis Felipe Bortolatto da Cunha
21 de setembro de 2020
1. Introdução
Este roteiro tem como objetivo auxiliar na análise de regressão linear no software R.
Para executar a regressão linear vamos utilizar dados demográficos e de consumo de água de 2010, para uma amostra de 4.417 municípios, extraídos do Censo Demográfico (IBGE) e Sistema Nacional de Informações sobre Saneamento (SNIS), para investigar em que medida o consumo de água está relacionado com a renda, conforme análise apresentada por Carmo et al., 2013.
Uma base de dados pré-processada está disponível no endereço abaixo: https://1drv.ms/u/s!AjettDH-3Gbni9kKwbI1eaRcG3MRRw?e=SdtiNN
1.1. Instalando e abrindo os pacotes
O software estatístico R já conta com funções básicas que permitem executar uma análise de regressão linear. Mas além das funções básicas, vamos usar os pacotes:
tidyverse: para adicionar o operador%>%e a funçãoselect()performance: para comparar modelosbroom: para salvar as estatísticas resultantes do modelo
Se você ainda não possui esses pacotes instalados, é necessário executar o comando abaixo para instalar.
Após a instalação, você pode importá-los com o uso da função library().
1.2. Importando os dados
A base de dados pode ser importada conforme as instruções do Roteiro 2.
Como ela está hospedada na nuvem, também pode ser importada com o endereço web como argumento, ao invés do endereço local, conforme o exemplo abaixo. É importante lembrar que a importação pelo endereço web exige conexão com a internet!
1.3. Explorando os dados
A função names() exibe os nomes das variáveis.
## [1] "ID_IBGE" "DOMICIL" "REDE" "PROPREDE" "ID_SNIS" "NOME_MUN"
## [7] "UF" "REGIAO" "PIB" "RENDAPITA" "GINI" "IDH"
## [13] "IDH_CLASS" "GE012" "AG001" "AG020" "AG022" "CONSUMO1"
## [19] "CONSUMO2"Como é possível ver após a execução do comando, os nomes das variáveis estão codificados. Mas a tabela abaixo apresenta uma descrição de cada variável.
| Código | Descrição |
|---|---|
| ID_IBGE | Código IBGE (7 dígitos) |
| DOMICIL | Quantidade de Domicílios |
| REDE | Quantidade de Domicílios com Acesso à Rede Geral de Água |
| PROPREDE | Proporção de Domicílios com com Acesso à Rede Geral de Água (REDE/DOMICIL) |
| ID_SNIS | Código IBGE (6 dígitos) |
| NOME_MUN | Nome do Município |
| UF | Unidade da Federação |
| REGIAO | Região do País |
| PIB | Produto Interno Bruto 2010 |
| RENDAPITA | Renda per Capita 2010 |
| GINI | Índice GINI 2010 |
| IDH | Índice de Desenvolvimento Humano 2010 |
| IDH_CLASS | Classificação do Índice de Desenvolvimento Humano 2010: Muito Alto >= 0,9; Alto >= 0,8; Médio >= 0,5; Baixo < 0,5. |
| GE012 | População Total Residente no Município |
| AG001 | População Total Atendida com Abastecimento de Água |
| AG020 | Volume Micromedido nas Economias Residenciais Ativas de Agua - 1.000 m3/ano |
| AG022 | Quantidade de Economias Residenciais Ativas Micromedidas |
| CONSUMO1 | Consumo de Água per capita - População Total - m3/ano (AG020/GE012) |
| CONSUMO2 | Consumo de Água per capita - População Atendida - m3/ano (AG020/AG001) |
2. Regressão linear simples
Com o código e descrição das variáveis já é possível executar uma regressão linear simples.
A regressão linear é uma ferramenta estatística que permite explorar e inferir a relação de uma variável dependente com variáveis independentes, a partir da fórmula:
\(Y = \beta_0 + \beta_1 X_1\)
Onde:
- Y = variável dependente (resposta/saída)
- X = variável independente (indicadora/explicativa/preditora)
- \(\beta_0\) = coeficiente do intercepto
- \(\beta_1\) = coeficiente da inclinação
Para executar uma regressão linear simples no R usaremos a função lm(), salvando os resultados no objeto modelo1. São necessários os argumentos: formula que indica as variaveis seguindo a estrutura “Y ~ X”, data que indica a base de dados e adicionalmente foi especificado o argumento na.action = na.exclude para excluir os valores faltantes (NA).
2.1. Análise do ajuste global e coeficientes do modelo
É possível exibir o ajuste global e coeficientes do modelo executando a função summary().
##
## Call:
## lm(formula = CONSUMO1 ~ RENDAPITA, data = agua_rede1, na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -57.650 -7.564 -0.668 6.631 164.515
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.25200 0.45548 9.335 <2e-16 ***
## RENDAPITA 0.04100 0.00082 49.997 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.14 on 4415 degrees of freedom
## (1149 observations deleted due to missingness)
## Multiple R-squared: 0.3615, Adjusted R-squared: 0.3614
## F-statistic: 2500 on 1 and 4415 DF, p-value: < 2.2e-16Primeiro sobre o ajuste. O R² (coeficiente de determinação) do modelo é 0,3615, que podemos interpretar como a medida que a variável preditora X explica a variação em Y. O R² ajustado é 0,314, sendo uma medida alternativa ao R² que penaliza a inclusão de variáveis independentes (X) pouco explicativas. O R² e R² ajustado indicam um bom ajuste do modelo.
O Teste F é uma medida do quanto o modelo melhorou na previsão de valores comparado com o nível de não precisão do modelo. O modelo apresenta um valor alta de estatística F (2500) e, ainda mais importante, é que o p-valor (0,00000000000000022) é menor que o nível de significância (0,05), ou seja, existe uma alta probabilidade de que os resultados do modelo não representem um erro amostral.
Sobre os coeficientes beta, a estimativa do \(\beta_0\) (coeficiente do intercepto) é 4,252 e do \(\beta_1\) (coeficiente da inclinação) é 0,041. A tabela de coeficientes também apresenta o erro padrão, estatística t e p-valor, que indicam que os resultados são significativos até para um nível de significância de 0,001.
Portanto, podemos inferir que, em média, um aumento de R$ 1,00 na renda per capita está associada a um aumento de 0,041 m³/ano no consumo de água (41 litros/ano).
2.2. Análise dos resíduos
O modelo de regressão linear assume algumas hipóteses sobre os dados de entrada para que seu resultado seja significativo.
- O modelo é linear nos parâmetros.
- A amostragem é aleatória.
- Variação amostral da variável independente.
- Média condicional do erro igual a zero.
- O erro tem a mesma variância para qualquer valor da variável explicativa.
Se o modelo for adequado, os resíduos devem refletir as propriedades impostas pelo termo de erro do modelo. Portanto, a análise dos resíduos se faz necessária para avaliar a adequação do modelo.
O modelo salvo permite a visualização de seis gráficos que auxiliam na análise dos resíduos e hipóteses.
Cada um desses gráficos pode ser acessado com a função plot(), com o modelo criado como o primeiro argumento e um segundo argumento which = com o número do gráfico a ser acessado.
O primeiro gráfico exibe a relação entre os resíduos e os valores ajustados. Ele permite verificar a não-linearidade do modelo. Um bom modelo possui os valores distribuidos em torno da linha de resíduos igual a zero (linha pontilhada). O modelo aplicado (linha vermelha) desvia um pouco do zero, principalmente nos extremos dos valores ajustados.
O segundo gráfico (QQ-plot) exibe os resíduos normalizados e os quantis teóricos da curva normal, ou seja, verifica a hipótese de normalidade dos resíduos. O ideal é que as observações sigam a linha pontilhada.
O terceiro gráfico é útil para verificar a hipótese de homocedasticidade dos resíduos. No modelo ideal, os pontos estão distribuídos uniformemente ao redor da linha vermelha.
O quarto gráfico é útil para detectar observações extremas, que possuem alta influência no modelo. As observações mais influentes possuem um valor de distância de Cook maior.
O quinto gráfico exibe a influência pelos resíduos padronizados. Neste gráfico, as observações mais distantes possuem maior influência.
O sexto gráfico também apresenta uma análise da distância de Cook, permitindo identificar pontos influentes, mas em termos da alavancagem.
2.3. Diagrama de dispersão com ajuste
Após a aplicação do modelo, é possível criar um diagrama de dispersão com a linha de tendência (ajuste). Para criar este gráfico, é só executar a função abline() após a criação do gráfico.
plot(x = agua_rede1$RENDAPITA,
y = agua_rede1$CONSUMO1,
xlab = "Renda per capita",
ylab = "Consumo de água per capita (m³/ano)")
abline(modelo1, col = "red") 
Também é possível criar esse gráfico com o pacote ggplot2, conforme indica o código abaixo:
ggplot(data = agua_rede1, aes(x = RENDAPITA, y = CONSUMO1)) +
geom_point() +
geom_smooth(data = lm(formula = CONSUMO1 ~ RENDAPITA, data = agua_rede1), method = "lm", col = "red", se = FALSE) +
theme_bw() +
xlab("Renda per capita") +
ylab("Consumo de água per capita (m³/ano)")
3. Regressão linear múltipla
É possível executar uma regressão linear múltipla de forma semelhante à regressão linear simples.
Executando o código abaixo é possível investigar se PROPREDE também influencia em CONSUMO1, além de RENDAPITA. Vamos salvar este modelo como `modelo2.
3.1. Resumo do modelo
Para acessar o ajuste e resumo do modelo, use a função summary().
##
## Call:
## lm(formula = CONSUMO1 ~ RENDAPITA + PROPREDE, data = agua_rede1,
## na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -43.399 -5.934 0.250 5.043 159.926
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.825e+01 6.982e-01 -26.15 <2e-16 ***
## RENDAPITA 2.892e-02 7.722e-04 37.45 <2e-16 ***
## PROPREDE 4.023e+01 1.032e+00 39.00 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.33 on 4414 degrees of freedom
## (1149 observations deleted due to missingness)
## Multiple R-squared: 0.5251, Adjusted R-squared: 0.5249
## F-statistic: 2440 on 2 and 4414 DF, p-value: < 2.2e-164. Comparando modelos
Uma forma de comparar modelos é usando o pacote performance. Para comparar os dois modelos (regressão linear simples e regressão linear múltipla), vamos usar a função compare_performance() e salvar a comparação como o objeto tabela_modelos.
## # Comparison of Model Performance Indices
##
## Model | Type | AIC | BIC | R2 | R2_adjusted | RMSE | BF | Performance_Score
## -------------------------------------------------------------------------------------------------
## modelo2 | lm | 33983.66 | 34009.24 | 0.52 | 0.52 | 11.33 | 1.33e+282 | 100.00%
## modelo1 | lm | 35289.28 | 35308.46 | 0.36 | 0.36 | 13.13 | 1.00 | 0.00%
##
## Model modelo2 (of class lm) performed best with an overall performance score of 100.00%.O modelo 2 (regressão linear múltipla) possui um R² e R² ajustado maior que o modelo 1 (regressão linear simples).
Se quiser salvar essa tabela, é possível fazê-lo com a função write.csv2().
Além das métricas apresentadas, também é importante avaliar as hipóteses e redíduos do segundo modelo, conforme as instruções apresentadas na seção 2.2.
5. Exportando os resultados do modelo
Ao executar a regressão linear simples ou múltipla, são computadas diversas estatísticas relacionadas a cada observação, por exemplo os valores previstos para Y e os resíduos padronizados.
Existem diversas formas de salvar essas estatísticas, mas recomendamos o uso da função augment() do pacote broom. Ela exige dois argumentos: x que especifica o modelo e data que especifica os dados usados para gerar o modelo.
É possível visualizar a base de dados nova com a função head().
## # A tibble: 6 x 25
## ID_IBGE DOMICIL REDE PROPREDE ID_SNIS NOME_MUN UF REGIAO PIB RENDAPITA
## <int> <int> <int> <dbl> <int> <chr> <chr> <chr> <dbl> <dbl>
## 1 3548807 50492 50472 1.00 354880 "S\xe3o~ SP Sudes~ 6.69e6 2009.
## 2 4214904 1312 412 0.314 421490 "Rio Fo~ SC Sul 3.90e4 1571.
## 3 3303302 169237 164768 0.974 330330 "Niter\~ RJ Sudes~ 5.83e6 1951.
## 4 3205309 108515 107715 0.993 320530 "Vit\xf~ ES Sudes~ 9.27e6 1821.
## 5 3547304 31610 28728 0.909 354730 "Santan~ SP Sudes~ 1.09e6 1798.
## 6 4205407 147437 137984 0.936 420540 "Floria~ SC Sul 4.28e6 1770.
## # ... with 15 more variables: GINI <dbl>, IDH <dbl>, IDH_CLASS <chr>,
## # GE012 <int>, AG001 <int>, AG020 <dbl>, AG022 <int>, CONSUMO1 <dbl>,
## # CONSUMO2 <dbl>, .fitted <dbl>, .resid <dbl>, .std.resid <dbl>, .hat <dbl>,
## # .sigma <dbl>, .cooksd <dbl>Foram adicionadas à base de dados original seis novas variáveis:
- .fitted = Valor ajustado ou previsto
- .resid = A diferença entre os valores observados e ajustados
- .std.resid = Resíduos padronizados
- .hat = Diagonal da matriz de influência
- .sigma = Desvio padrão residual estimado a observação correspondente é excluída do modelo
- .cooksd = Distância de Cook
É possível exportar essa nova tabela com a função write.csv2. Ela será usada nas próximas aulas para visualizar a distribuição espacial dos resíduos padronizados.
6. Regressão linear múltipla com o método stepwise
A regressão stepwise consiste em adicionar e remover iterativamente preditores (X) no intuito de encontro o subconjunto de variáveis que resulta no melhor desempenho, que é o modelo que melhor reduz o erro de predição (AIC).
Ela pode ser feita de três formas diferentes:
- Forward selection: começa sem preditores no modelo e adiciona iterativamente os preditores com maior poder preditivo, até não existir melhoria estatisticamente significativa.
- Backward selection: começa com todos os preditores no modelo e remove iterativamente os preditores com menor poder preditivo, resultando em um modelo onde todos os preditores são estatisticamente significativos.
- Stepwise selection: é uma combinação da forward selection e backward selection, sendo o método mais eficaz.
Para mostrar como funciona esse método, vamos aplicar a regressão do tipo stepwise selection em um conjunto de variáveis da base de dados. As variáveis escolhidas foram: RENDAPITA, PROPREDE, PIB, IDH e GINI. Para aplicar o método, primeiro é necessário criar o modelo com todas as variáveis e depois aplicar a função step() para obter o modelo com os melhores resultados.
modelo3 <- lm(formula = CONSUMO1 ~ RENDAPITA + PROPREDE + PIB + IDH + GINI, data = agua_rede1, na.action = na.exclude)##
## Call:
## lm(formula = CONSUMO1 ~ RENDAPITA + PROPREDE + PIB + IDH + GINI,
## data = agua_rede1, na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -43.474 -5.665 0.301 4.870 157.762
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.378e+00 2.468e+00 -3.395 0.000693 ***
## RENDAPITA 2.593e-02 1.085e-03 23.905 < 2e-16 ***
## PROPREDE 3.900e+01 1.031e+00 37.837 < 2e-16 ***
## PIB -7.854e-08 5.744e-08 -1.367 0.171574
## IDH 7.347e+00 2.721e+00 2.700 0.006955 **
## GINI -2.517e+01 2.858e+00 -8.809 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.19 on 4411 degrees of freedom
## (1149 observations deleted due to missingness)
## Multiple R-squared: 0.5367, Adjusted R-squared: 0.5362
## F-statistic: 1022 on 5 and 4411 DF, p-value: < 2.2e-16## Start: AIC=21343.2
## CONSUMO1 ~ RENDAPITA + PROPREDE + PIB + IDH + GINI
##
## Df Sum of Sq RSS AIC
## - PIB 1 234 552927 21343
## <none> 552693 21343
## - IDH 1 914 553607 21349
## - GINI 1 9723 562416 21418
## - RENDAPITA 1 71603 624296 21879
## - PROPREDE 1 179380 732073 22583
##
## Step: AIC=21343.08
## CONSUMO1 ~ RENDAPITA + PROPREDE + IDH + GINI
##
## Df Sum of Sq RSS AIC
## <none> 552927 21343
## + PIB 1 234 552693 21343
## - IDH 1 956 553884 21349
## - GINI 1 10325 563252 21423
## - RENDAPITA 1 72685 625612 21887
## - PROPREDE 1 179171 732098 22581##
## Call:
## lm(formula = CONSUMO1 ~ RENDAPITA + PROPREDE + IDH + GINI, data = agua_rede1,
## na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -43.474 -5.711 0.267 4.874 157.771
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.066709 2.457560 -3.282 0.00104 **
## RENDAPITA 0.025644 0.001065 24.083 < 2e-16 ***
## PROPREDE 38.931694 1.029642 37.811 < 2e-16 ***
## IDH 7.509447 2.718485 2.762 0.00576 **
## GINI -25.702970 2.831787 -9.077 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.19 on 4412 degrees of freedom
## (1149 observations deleted due to missingness)
## Multiple R-squared: 0.5366, Adjusted R-squared: 0.5361
## F-statistic: 1277 on 4 and 4412 DF, p-value: < 2.2e-16Neste caso, a aplicação do método stepwise excluiu a variável PIB, que não representava ganhos na estatística AIC (critério de informação de Akaike) e, coincidentemente, não era significativa (pois o p-valor > 0,05).
Mas é necessário cautela ao aplicar o método stepwise, pois ele permite não pensar no problema o qual buscamos a solução. Além disso, ele pode apresentar: (1) viés na estimativa dos parâmetros, (2) o problema inerente (mas frequentemente esquecido) do teste de múltiplas hipóteses e (3) foco inadequado na busca de um único modelo. Tomando o exemplo apresentado, sabemos que RENDAPITA, PIB, IDH e GINI são altamente correlacionados e não deveriam ser considerados variáveis explicativas de um fenômeno conjuntamente.