Leitura dos dados

base <-read_excel("~/Library/CloudStorage/OneDrive-Pessoal/Documentos/7. Especialização/1. Disciplinas/REGRESSÃO NÃO-LINEAR/Aula1/2022_Aula01/2022_Lista01/baseDEC_Compensacoes_CEMIG2018.xlsx")
#View(base)

Analise exploratoria da base de dados

Variável resposta de interesse: Compensacoes.Pagas

packHV::hist_boxplot(base$Compensacoes.Pagas, main="", col="light blue", xlab="Compensacoes.Pagas")

Análise do gráfico entre a variável resposta e as demais presentes no banco de dados

Variáveis analisadas
Variáveis
Volume.chuva
Descargas.atm
Vento
umidade
temperatura
AG1
AE1
AE2
VC1
DS1
DS2
AR1
grupos2
Unidades.Consumidoras 
 par(mfrow=c(5,3))
   plot(Compensacoes.Pagas ~ Volume.chuva, data=base, pch=19, col="dark blue")
   plot(Compensacoes.Pagas ~ Descargas.atm, data=base, pch=19, col="dark blue")
   plot(Compensacoes.Pagas ~ Vento, data=base, pch=19, col="dark blue")
   plot(Compensacoes.Pagas ~ umidade, data=base, pch=19, col="dark blue")
   plot(Compensacoes.Pagas ~ temperatura, data=base, pch=19, col="dark blue")
   plot(Compensacoes.Pagas ~ AG1, data=base, pch=19, col="dark blue")
   plot(Compensacoes.Pagas ~ AE1, data=base, pch=19, col="dark blue")
   plot(Compensacoes.Pagas ~ AE2, data=base, pch=19, col="dark blue")
   plot(Compensacoes.Pagas ~ VC1, data=base, pch=19, col="dark blue")
   plot(Compensacoes.Pagas ~ DS1, data=base, pch=19, col="dark blue")
   plot(Compensacoes.Pagas ~ DS2, data=base, pch=19, col="dark blue")
   plot(Compensacoes.Pagas ~ AR1, data=base, pch=19, col="dark blue")
   plot(Compensacoes.Pagas ~ as.factor(grupos2), data=base, pch=19, col="dark blue", xlab="Grupos")
   plot(Compensacoes.Pagas ~ Unidades.Consumidoras, data=base, pch=19, col="dark blue")

Matriz de correlacao

Será plotada a matriz de correlação somente com as variáveis desejadas

variaveis = c("Compensacoes.Pagas","Volume.chuva", "Descargas.atm","Vento","umidade","temperatura","AG1","AE1","AE2","VC1","DS1","DS2","AR1","grupos2","Unidades.Consumidoras")

base1 = base[variaveis]

mx <- cor(base1, method = "spearman")
 corrplot::corrplot(mx, method = 'square', type="upper", 
          order="AOE", col= RColorBrewer::brewer.pal(n=8, name="RdYlBu"),diag=FALSE)

É possível observar que Compensacoe.Pagas tem correlações mais forte com:

  • DS1

  • DS2

  • AE2

  • Unidades.Consumidoras

Ajuste de modelos univariados

 # Ajuste com as variaveis de interesse
 saida <- exploreR::masslm(base1, dv.var = "Compensacoes.Pagas")

 # Ordena modelo pelo decrescente do  R quadrado
 saida <- saida[order(saida$R.squared, decreasing = TRUE),]
 
 # Nomeia cabecalho da tabela
 names(saida) <- c("Variavel", "Coeficiente", "p-valor", "R2")
 
 # Tabela
 knitr::kable(saida, digits = 2, row.names = FALSE,
              format.args = list(scientific = FALSE))
Variavel Coeficiente p-valor R2
DS2 130500.00 0.00 0.37
AE1 105400.00 0.00 0.24
AE2 105100.00 0.00 0.24
DS1 90050.00 0.00 0.18
AR1 89050.00 0.00 0.17
AG1 78680.00 0.00 0.14
Unidades.Consumidoras 2.33 0.00 0.07
Descargas.atm 17.33 0.00 0.06
grupos2 67810.00 0.01 0.03
temperatura 18630.00 0.12 0.01
VC1 -12930.00 0.32 0.00
umidade -2611.00 0.39 0.00
Vento 3067.00 0.93 0.00
Volume.chuva -4.26 0.98 0.00

Ajuste de modelos univariados - escala logarítmica

# nova base
base_log = as.data.frame(log(base1))

#substituir NaN por zero
base_log[is.na(base_log)] = 0

#View(base_log)

saida_log =  exploreR::masslm(base_log, dv.var = "Compensacoes.Pagas")
saida_log <- saida_log[order(saida_log$R.squared, decreasing = TRUE),]

# Tabela
names(saida_log) <- c("Variavel", "Coeficiente", "p-valor", "R2")
knitr::kable(saida_log, digits = 2, row.names = FALSE,
              format.args = list(scientific = FALSE))
Variavel Coeficiente p-valor R2
Unidades.Consumidoras 0.63 0.00 0.30
Descargas.atm 0.14 0.00 0.07
grupos2 0.38 0.02 0.02
Volume.chuva 0.47 0.07 0.01
umidade -0.83 0.36 0.00
AG1 0.06 0.43 0.00
AR1 0.05 0.50 0.00
Vento -0.13 0.56 0.00
AE2 0.04 0.63 0.00
VC1 0.03 0.69 0.00
AE1 0.01 0.92 0.00
DS1 0.00 0.97 0.00
temperatura -0.04 0.97 0.00
DS2 0.00 1.00 0.00

Ajuste de modelo com regressão múltipla

modelo1 = lm(Compensacoes.Pagas ~ ., data=base1)
summary(modelo1)
## 
## Call:
## lm(formula = Compensacoes.Pagas ~ ., data = base1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -328615  -53674  -10267   26291 1863863 
## 
## Coefficients: (1 not defined because of singularities)
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           -3.594e+05  5.597e+05  -0.642   0.5214    
## Volume.chuva          -1.453e+01  1.665e+02  -0.087   0.9305    
## Descargas.atm          6.485e+00  4.329e+00   1.498   0.1354    
## Vento                  6.209e+03  2.757e+04   0.225   0.8220    
## umidade               -1.262e+02  3.346e+03  -0.038   0.9699    
## temperatura            2.345e+04  1.488e+04   1.576   0.1163    
## AG1                   -1.048e+05  2.259e+04  -4.640 5.60e-06 ***
## AE1                    2.291e+04  2.945e+04   0.778   0.4374    
## AE2                    6.625e+03  2.824e+04   0.235   0.8147    
## VC1                           NA         NA      NA       NA    
## DS1                    8.798e+04  3.832e+04   2.296   0.0225 *  
## DS2                    1.581e+05  2.220e+04   7.123 1.11e-11 ***
## AR1                   -2.554e+03  1.730e+04  -0.148   0.8828    
## grupos2                2.674e+04  3.031e+04   0.882   0.3785    
## Unidades.Consumidoras -2.983e+00  1.385e+00  -2.153   0.0322 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 160500 on 252 degrees of freedom
## Multiple R-squared:  0.4632, Adjusted R-squared:  0.4356 
## F-statistic: 16.73 on 13 and 252 DF,  p-value: < 2.2e-16

Análise de resíduos e do ajuste

hist_boxplot(residuals(modelo1))

 par(mfrow=c(2,2))
   plot(modelo1, pch=17, col=" dark blue", lty=0)

 # Teste de Normalidade dos erros
 shapiro.test(residuals(modelo1))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modelo1)
## W = 0.55304, p-value < 2.2e-16

Há evidências de que os resíduos do modelo não é normal.

Ajuste do modelo com seleção de variáveis

   modelo_final = step(lm(Compensacoes.Pagas  ~ ., data = base1), trace = FALSE)
   summary(modelo_final)
## 
## Call:
## lm(formula = Compensacoes.Pagas ~ Descargas.atm + temperatura + 
##     AG1 + DS1 + DS2 + Unidades.Consumidoras, data = base1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -342819  -56821  -10286   27575 1874411 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           -5.550e+05  2.318e+05  -2.394 0.017359 *  
## Descargas.atm          5.336e+00  3.784e+00   1.410 0.159734    
## temperatura            3.241e+04  9.415e+03   3.442 0.000673 ***
## AG1                   -9.162e+04  1.765e+04  -5.189 4.26e-07 ***
## DS1                    8.951e+04  3.296e+04   2.715 0.007064 ** 
## DS2                    1.688e+05  1.981e+04   8.521 1.32e-15 ***
## Unidades.Consumidoras -2.557e+00  1.204e+00  -2.124 0.034610 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 158800 on 259 degrees of freedom
## Multiple R-squared:  0.4597, Adjusted R-squared:  0.4472 
## F-statistic: 36.73 on 6 and 259 DF,  p-value: < 2.2e-16
  par(mfrow = c(2,2))
   plot(modelo_final, col ="red", pch = 18 )

Previsao do modelo

base1$fit = predict(modelo_final,newdata = base1)
 plot(fit ~ Compensacoes.Pagas , data=base1, pch=19, col="blue",
      xlab = "Variavel resposta", ylab ="Predito")

Consideracões finais e conclusão

O modelo final de regressão linear múltipla para estimar o comportamento médio da variável Compensacoes.Pagas identificou 6 variáveis preditoras que não apresentam comportamento linear. O coeficiente de determinação doo modelo é de 45,97%, indicando um modelo com baixa explicação do comportamento médio das compensações pagas para os clientes da CEMIG.

Referêcia bibliografica

Harrison, David, and Rubinfeld, Daniel L., Hedonic Housing Prices and the Demand for Clean Air, Journal of Environmental Economics and Management, Volume 5, (1978), 81-102.