---
title: "Caso S07: Regressão linear simples e múltipla"
output:
  html_document:
    df_print: paged
---

------------------------------------------------------------------------

### Carregamento dos pacotes do R utilizados:


```{r, message=FALSE, warning=FALSE, include = FALSE}
library(readxl)
library(corrplot)
library(ggplot2)
library(dplyr)
library(ggpubr)
library(caret)
```

------------------------------------------------------------------------

### 0. Seleção dos dados:

```{r include = TRUE}
df_RG = read_excel("Df_ProdEnerg.xlsx")
df_RG = data.frame(df_RG)
str(df_RG)
```

------------------------------------------------------------------------

### 1. Análise exploratória e diagramas de dispersão:

```{r}
hist(df_RG$PE, breaks=10, col = "red3", main="Power Output (Y)")
hist(df_RG$AT, breaks=50, col = "brown3", main="Ambient Temperature")
hist(df_RG$V, breaks=10, col = "blue3", main="Exhaust Steam Pressure")
hist(df_RG$AP, breaks=10, col = "green3", main="Atmospheric Pressure")
hist(df_RG$RH, breaks=10, col = "yellow3", main="Relative Humidity")
```

   
  
```{r}
corrplot.mixed(cor(df_RG),upper="ellipse")
```   
   
   
```{r}
df_RG %>% ggplot(aes(x=AT,y=PE)) + geom_point() +
  xlab("Temperatura ambiente (oC)") + ylab("Energia produzida (MW)") +
  geom_smooth()
```

```{r}
df_RG %>% ggplot(aes(x=V,y=PE)) + geom_point() +
  xlab("Pressao do vapor (cm Hg)") + ylab("Energia produzida (MW)")+
  geom_smooth()
```


```{r}
df_RG %>% ggplot(aes(x=AP,y=PE)) + geom_point() +
  xlab("Pressao atmosferica (mbar)") + ylab("Energia produzida (MW)")+
  geom_smooth()
```


```{r}
df_RG %>% ggplot(aes(x=RH,y=PE)) + geom_point() +
  xlab("Umidade relativa (%)") + ylab("Energia produzida (MW)")+
  geom_smooth()
```


___

## Caso Bidimensional:   

------------------------------------------------------------------------

### 2. Definição do modelo geral


```{r}
df_RG %>% ggplot(aes(x=AT,y=PE)) + geom_point() +
  xlab("Temperatura ambiente (oC)") + ylab("Energia produzida (MW)") +
  geom_smooth(method = lm, formula = y ~ x) + stat_regline_equation(label.x = 5, label.y = 420)
```
     
    

------------------------------------------------------------------------

### 3 a 6. Estimação dos parâmetros a precisão da previsão


#### Divisao dos dados em amostras de treinamento e de teste do modelo:

Nesta etapa os dados serao divididos de forma a ficar 80% para treinamento do modelo e 20% para teste do modelo. Assim:


```{r}
set.seed(136)
Index <- createDataPartition(y = df_RG$PE,p = 0.8, list = FALSE)
treino <- df_RG[Index,]
teste <- df_RG[-Index,]
```

#### Criação do modelo:

```{r include = TRUE}
fit.lm <- train(PE ~ AT, data = treino,
                method = "lm",
                trControl =trainControl(method = "cv",number = 10))
fit.lm
summary(fit.lm)
```

- Análise dos resíduos:   
   
```{r}
Prev.teste.lm <- predict(fit.lm,newdata = teste)
predVals <- extractPrediction(list(fit.lm), 
                              testX = teste[,1:4], 
                              testY = teste[,5])
plotObsVsPred(predVals)
postResample(pred = Prev.teste.lm, obs = teste$PE)
```


```{r}
a <- data.frame(y.hat=Prev.teste.lm,y=teste$PE)
cor(a$y.hat,a$y)**2
```

```{r}
fit.lm <- lm(PE ~ ., data = treino)
plot(fit.lm)
library(lmtest) # testes: breusch-godfrey, breuch-pagan
library(tseries) # teste jarque-bera
bgtest(fit.lm) # Teste de autocorrelação (Breuch-Godfrey)
bptest(fit.lm) # Teste de homoscedasticidade (Breuch-Pagan)
jarque.bera.test(fit.lm$residuals) # Teste de normalidade dos resíduos
```
Como hipótese nula de não existência de autocorrelação foi rejeitada, existe evidência de um problema de autocorrelação nos resíduos dessa regressão.
O resultado do teste de heteroscedasticidade Breuch-Pagan mostra que não se pode rejeitar a hipótese nula de homocedasticidade dos resíduos, assim, não há evidência de heteroscedasticidade nos resíduos dessa regressão.
O teste de normalidade Jarque-Bera mostra que não há evidência de não normalidade dos resíduos.


------------------------------------------------------------------------


## Regressão múltipla:   

```{r include = TRUE}
fit.lm <- train(PE ~ ., data = treino,
                method = "lmStepAIC",
                trControl =trainControl(method = "cv",number = 10), 
                trace = FALSE)
fit.lm
summary(fit.lm)
```

```{r}
nova.base <- data.frame(rbind(c(14.96,41.76,1024.07,73.17),
                              c(14.96,41.76,1024.07,74.17)))
colnames(nova.base) <- c("AT","V","AP","RH")
predict(fit.lm,newdata = nova.base)
```

- VIF:

```{r}
car::vif(fit.lm$finalModel)
```
   
- Regressão sem RH:   
   
```{r include = TRUE}
fit.lm <- train(PE ~ . -RH, data = treino,
                method = "lmStepAIC",
                trControl =trainControl(method = "cv",number = 10), 
                trace = FALSE)
fit.lm
summary(fit.lm)
```

- Análise dos resíduos:   
   
```{r}
Prev.teste.lm <- predict(fit.lm,newdata = teste)
predVals <- extractPrediction(list(fit.lm), 
                              testX = teste[,1:4], 
                              testY = teste[,5])
plotObsVsPred(predVals)
postResample(pred = Prev.teste.lm, obs = teste$PE)
```


