Projeto Integrador - Análise de Regressão Linear (CE071)
Helen Lourenço
Maria Eduarda Mochinski
Vitor Kroeff
Resumo
Como projeto integrador da disciplina de Análise de Regressão Linear, e com auxílio do software R, coletamos e tratamos dados do Instituto Paranaense de Desenvolvimento Econômico e Social (IPARDES), com a intenção de explicar a variação da probabilidade de sobrevivência até os sessenta anos, por município do Paraná, através de um modelo de regressão linear. Com algumas limitações, utilizamos técnicas de diagnóstico e medidas corretivas trabalhadas em sala de aula para otimizar a performance do ajuste.
Introdução
Com o poder de escolher as variáveis mais adequadas para o modelo, para incorporar o tamanho do município no ajuste, utilizamos o logaritmo da população de cada um deles, com o objetivo de amenizar o efeito de escala. Além disso, decidimos trabalhar com indicadores de educação, desenvolvimento econômico e urbanização das cidades. Nos primeiros tópicos do documento, será possível entender o comportamento individual de cada uma dessas variáveis (especificadas em Dados), além de visualizar a relação individual das mesmas com a resposta. Já em Ajuste, temos um modelo aditivo utilizando todas as covariáveis, que, nas próximas sessões, será analisado criteriosamente. Como estamos trabalhando com uma regressão linear, precisamos ter cuidado e garantir que certas suposições sejam satisfeitas. Assim, precisamos analisar o comportamento dos erros e resíduos do modelo, além de procurar por observações atípicas e checar a multicolinearidade entre as variáveis.
Dados
Disponíveis em: http://www.ipardes.gov.br/imp/
Variável resposta (por município):
- sob: probabilidade de uma criança recém-nascida viver até os 60 anos * 100
Covariáveis (por município):
mat: proporção de matrículas na educação básica
vul: proporção de vulneráveis à pobreza
urb: grau de urbanização
gini: índice de gini
renda: renda média domiciliar per capita
log_pop: logaritmo da população estimada
Primeiras linhas da base:
head(dados)## sob mat vul urb gini renda log_pop
## 1 83.48 25.15 19.70 85.33 0.53 870.59 16.18
## 2 80.80 23.36 29.50 73.83 0.44 533.79 8.98
## 3 82.28 26.84 50.03 32.31 0.53 430.79 8.83
## 4 79.82 24.76 36.56 34.12 0.48 516.96 9.08
## 5 84.60 23.19 18.33 95.82 0.43 629.58 11.49
## 6 78.56 31.96 48.16 49.58 0.58 504.11 8.24Análise Descritiva
summary(dados)## sob mat vul urb
## Min. :77.37 Min. :17.97 Min. : 4.48 Min. : 9.35
## 1st Qu.:81.16 1st Qu.:23.64 1st Qu.:19.56 1st Qu.: 57.42
## Median :82.64 Median :25.41 Median :27.20 Median : 72.57
## Mean :82.49 Mean :25.56 Mean :28.27 Mean : 69.23
## 3rd Qu.:83.92 3rd Qu.:27.17 3rd Qu.:35.31 3rd Qu.: 85.12
## Max. :86.49 Max. :40.92 Max. :63.97 Max. :100.00
## NA's :35 NA's :35
## gini renda log_pop
## Min. :0.3300 Min. : 274.1 Min. : 7.270
## 1st Qu.:0.4300 1st Qu.: 507.7 1st Qu.: 8.585
## Median :0.4700 Median : 582.1 Median : 9.320
## Mean :0.4658 Mean : 600.7 Mean : 9.576
## 3rd Qu.:0.5000 3rd Qu.: 673.5 3rd Qu.:10.085
## Max. :0.6600 Max. :1536.4 Max. :16.180
## NA's :35 NA's :35Presença de valores nulos, distância entre os pontos mínimo e máximo nas variáveis vul e urb (não preocupantes por conta da distribuição dos quartis), e possível simetria (com exceção da variável urb).
Boxplot
Além de um intervalo de observações pequeno e leve assimetria na resposta, o gráfico não indica valores atípicos.
Histogramas
Variável resposta:
A resposta apresenta um pouco de simetria, mas acaba tendo observações bem distribuídas em um pequeno intervalo. Com um range pequeno (78:86), esse comportamento é esperado.
Covariáveis:
As variáveis explicativas expressam certa assimetria (menos significativa em mat e gini). urb se destaca com um formato de “escada”.
Gráficos de Dispersão
Com os gráficos acima, vemos que a variável vul apresenta uma relação linear mais forte com a resposta, seguida por urb e gini, que, apesar disso, acabam demonstrando um padrão mais aleatório e preocupante. mat também é assim, mas tem pontos menos dispersos, e renda e log_pop podem indicar não linearidade.
Ajuste do Modelo
ajuste <- lm(sob ~ mat + vul + urb + gini + renda + log_pop, data = dados)
summary(ajuste)##
## Call:
## lm(formula = sob ~ mat + vul + urb + gini + renda + log_pop,
## data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.0531 -1.1858 -0.0399 1.1634 3.4851
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 80.3981299 1.1096712 72.452 < 2e-16 ***
## mat -0.0354469 0.0299301 -1.184 0.23700
## vul -0.0784317 0.0258647 -3.032 0.00259 **
## urb 0.0039112 0.0059757 0.655 0.51316
## gini 5.9646972 3.0647114 1.946 0.05234 .
## renda 0.0002582 0.0017610 0.147 0.88349
## log_pop 0.2160577 0.0933653 2.314 0.02118 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.603 on 393 degrees of freedom
## (35 observations deleted due to missingness)
## Multiple R-squared: 0.2961, Adjusted R-squared: 0.2853
## F-statistic: 27.55 on 6 and 393 DF, p-value: < 2.2e-16Com o nosso ajuste inicial, atingimos um coeficiente de determinação de 0,2961, que pode ser considerado baixo, e um p-valor muito pequeno (ideal). Além disso, as estimativas de mínimos quadrados para as variáveis urb e renda indicam uma contribuição muito baixa das mesmas no ajuste.
Diagnóstico do Ajuste
No gráfico Residuals vs Fitted, observamos pontos dispersos aleatóriamente e próximos da linha de tendência, que é quase constante em torno de zero. Além de validar a não presença de padrões sistemáticos e outliers, nada indica que os erros não têm variância constante.
Já no Q-Q Residuals, conseguimos atestar a suposição de normalidade.
No Scale-Location, alternativa ao Residuals vs Fitted, que se baseia nos resíduos padronizados, também não encontramos evidências de que os erros não têm variância constante, por conta de seu padrão aleatório.
Por fim, o gráfico Cook’s distance registra valores elevados para as observações 1, 168, e 273, possíveis observações influentes, que precisam ser verificadas.
Observações Influentes
Para verificar a influência das observações 1, 168, e 273 no ajuste do modelo, vamos ajustar um novo modelo sem as mesmas e checar o seu desempenho.
Removendo as três observações:
ajuste2 <- lm(sob ~ mat + vul + urb + gini + renda + log_pop, data = dados[-c(1, 168, 273),])
summary(ajuste2)##
## Call:
## lm(formula = sob ~ mat + vul + urb + gini + renda + log_pop,
## data = dados[-c(1, 168, 273), ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.0565 -1.1623 -0.0048 1.1833 3.4784
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 80.3474534 1.1202725 71.721 < 2e-16 ***
## mat -0.0510922 0.0314581 -1.624 0.10515
## vul -0.0741283 0.0259871 -2.853 0.00457 **
## urb 0.0053042 0.0060737 0.873 0.38303
## gini 6.0614808 3.0559926 1.983 0.04802 *
## renda 0.0003863 0.0017667 0.219 0.82701
## log_pop 0.2272021 0.1003460 2.264 0.02411 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.598 on 390 degrees of freedom
## (35 observations deleted due to missingness)
## Multiple R-squared: 0.2996, Adjusted R-squared: 0.2888
## F-statistic: 27.8 on 6 and 390 DF, p-value: < 2.2e-16compareCoefs(ajuste, ajuste2)## Calls:
## 1: lm(formula = sob ~ mat + vul + urb + gini + renda + log_pop, data =
## dados)
## 2: lm(formula = sob ~ mat + vul + urb + gini + renda + log_pop, data =
## dados[-c(1, 168, 273), ])
##
## Model 1 Model 2
## (Intercept) 80.40 80.35
## SE 1.11 1.12
##
## mat -0.0354 -0.0511
## SE 0.0299 0.0315
##
## vul -0.0784 -0.0741
## SE 0.0259 0.0260
##
## urb 0.00391 0.00530
## SE 0.00598 0.00607
##
## gini 5.96 6.06
## SE 3.06 3.06
##
## renda 0.000258 0.000386
## SE 0.001761 0.001767
##
## log_pop 0.2161 0.2272
## SE 0.0934 0.1003
## Com pouca variação no \(R^2\) e no p-valor dos testes, e vendo que os coeficientes da regressão mudam muito pouco, com erros padrões quase idênticos, concluímos que a variação na qualidade do ajuste do modelo com a exclusão dessas observações é mínima, e decidimos continuar com o ajuste inicial.
Testes
Apesar de não recomendados em substituição à análise gráfica, decidimos performar alguns testes que dizem respeito aos resíduos do ajuste.
Shapiro-Wilk
Hipótese nula: resíduos têm distribuição normal
shapiro.test(rstudent(ajuste))##
## Shapiro-Wilk normality test
##
## data: rstudent(ajuste)
## W = 0.99224, p-value = 0.03544Com a rejeição da hipótese nula, ao nível de significância de 5%, os resíduos apresentam algum desvio da normalidade. Para verificar a magnitude desse desvio, vamos analisar graficamente a distribuição dos resíduos.
Apesar do que foi apontado pelo teste, acreditamos, por conta do gráfico acima, que esse desvio não é significativo.
ncvTest
Hipótese nula: resíduos têm variância constante
ncvTest(ajuste)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 5.482556, Df = 1, p = 0.019207Com a rejeição da hipótese nula, ao nível de significância de 5%, a variância dos resíduos não é nula. Assim como no caso acima, com base na análise gráfica da distribuição dos resíduos, não consideramos que esse seja um fator problemático.
Durbin-Watson
Hipótese nula: resíduos não apresentam autocorrelação
dwtest(ajuste)##
## Durbin-Watson test
##
## data: ajuste
## DW = 2.0166, p-value = 0.564
## alternative hypothesis: true autocorrelation is greater than 0Ao nível de significância de 5%, não temos evidências para rejeitar a hipótese nula. Sendo assim, os resíduos não apresentam autocorrelação.
Box-Cox
Com a intenção de verificar possíveis transformações pertinentes aos dados, vamos utilizar o método de Box-Cox.
b1 <- boxCox(ajuste, lambda = seq(0,1,0.1))
b1$x[which.max(b1$y)]## [1] 1Como \(\lambda = 1\), a escala original parece ser a mais indicada para o ajuste, não precisando de nenhuma transformação.
Multicolinearidade
Matriz de correlação:
A relação linear negativa entre renda e vul, e urb e vul, pode ser preocupante. Interpretação: Maior renda, menor vulnerabilidade à pobreza (o que faz sentido). Maior urbanização, menor vulnerabilidade à pobreza.
VIF do ajuste:
vif(ajuste)## mat vul urb gini renda log_pop
## 1.123326 13.622051 2.272173 4.759197 10.583778 1.671268Com dois valores um pouco acima de 10, temos que nos preocupar com as variáveis renda e vul.
Para avaliar o efeito da multicolinearidade, vamos perturbar levemente a variável resposta, adicionando erros aleatórios com distribuição Normal(0,1) a cada valor de y.
ajuste3 <- lm(sob + rnorm(435, 0, 1) ~ . , data = dados)
summary(ajuste3)##
## Call:
## lm(formula = sob + rnorm(435, 0, 1) ~ ., data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.683 -1.308 -0.003 1.345 4.460
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 80.9695724 1.2618733 64.166 < 2e-16 ***
## mat -0.0518404 0.0340353 -1.523 0.128529
## vul -0.0989268 0.0294123 -3.363 0.000845 ***
## urb 0.0002527 0.0067953 0.037 0.970358
## gini 7.5886661 3.4850661 2.177 0.030039 *
## renda -0.0010400 0.0020025 -0.519 0.603789
## log_pop 0.2846414 0.1061713 2.681 0.007650 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.823 on 393 degrees of freedom
## (35 observations deleted due to missingness)
## Multiple R-squared: 0.2568, Adjusted R-squared: 0.2455
## F-statistic: 22.64 on 6 and 393 DF, p-value: < 2.2e-16compareCoefs(ajuste, ajuste2)## Calls:
## 1: lm(formula = sob ~ mat + vul + urb + gini + renda + log_pop, data =
## dados)
## 2: lm(formula = sob ~ mat + vul + urb + gini + renda + log_pop, data =
## dados[-c(1, 168, 273), ])
##
## Model 1 Model 2
## (Intercept) 80.40 80.35
## SE 1.11 1.12
##
## mat -0.0354 -0.0511
## SE 0.0299 0.0315
##
## vul -0.0784 -0.0741
## SE 0.0259 0.0260
##
## urb 0.00391 0.00530
## SE 0.00598 0.00607
##
## gini 5.96 6.06
## SE 3.06 3.06
##
## renda 0.000258 0.000386
## SE 0.001761 0.001767
##
## log_pop 0.2161 0.2272
## SE 0.0934 0.1003
## Como o \(R^2\) do ajuste3 acabou sendo inferior ao do ajuste original, e as estimativas pontuais geradas pelos dois modelos são muito parecidas, não temos indicação significativa de multicolinearidade.
Resíduos Parciais
Os resíduos parciais permitem avaliar a relação entre a resposta e uma particular covariável ajustado o efeito das demais covariáveis.
Com a figura acima, fica claro que vul (vulnerabilidade à pobreza), uma vez que ajustado o efeito das demais covariáveis, apresenta uma relação linear (nesse caso, negativa) com a resposta. As outras covariáveis não têm um ajuste satisfatório, sobretudo urb e renda, com linhas de tendência quase constantes em zero.
Por fim, retirando as variáveis urb e renda do modelo, temos uma performance muito similar.
ajuste4 <- lm(sob ~ mat + vul + gini + log_pop, data = dados)
summary(ajuste4)##
## Call:
## lm(formula = sob ~ mat + vul + gini + log_pop, data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.0839 -1.1665 -0.0113 1.1955 3.5235
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 80.633538 1.014116 79.511 < 2e-16 ***
## mat -0.032176 0.029407 -1.094 0.274546
## vul -0.085946 0.008635 -9.954 < 2e-16 ***
## gini 6.187284 1.755549 3.524 0.000474 ***
## log_pop 0.238894 0.082648 2.891 0.004058 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.6 on 395 degrees of freedom
## (35 observations deleted due to missingness)
## Multiple R-squared: 0.2953, Adjusted R-squared: 0.2882
## F-statistic: 41.38 on 4 and 395 DF, p-value: < 2.2e-16Resultados e discussão
Com um coeficiente de determinação de 0,29, a variação na probabilidade de sobrevivência até os sessenta anos não foi bem explicada pelas variáveis escolhidas. Com a utilização do método de Box-Cox, esperávamos encontrar uma transformação que melhorasse a performance do nosso modelo, mas não foi o caso. Com isso, além da inclusão de novas variáveis explicativas, modelos não lineares podem ser uma alternativa. Por fim, mesmo com um ajuste não satisfatório, conseguimos atender as suposições necessárias para um modelo de regressão linear.