Regressão

Atividade

Escolher um conjunto de dados do link: UCI Machine Learning Repository

Escolher técnica de Machine Learning Regression e ordene por ano
Escolher uma base do ano de 2019 e praticar:
Utilize a função de regressão:
- Análise do problema
- Variáveis independente e dependente
- lm
- Correlação
- Gráfico de dispersão

###### PACOTES E FUNÇÕES AUXILIARES ######

library("dplyr", warn.conflicts = F)
library("kableExtra", warn.conflicts = F)
library("data.table", warn.conflicts = F)
library("DT", warn.conflicts = F)
library("tidyverse", warn.conflicts = F)
library("formattable", warn.conflicts = F)
library("ggplot2",warn.conflicts = F)
library("GGally",warn.conflicts = F)
library("lattice",warn.conflicts = F)
library("PRROC",warn.conflicts = F)
library("caret",warn.conflicts = F)
library("tidyverse",warn.conflicts = F)
library("profvis",warn.conflicts = F)
library("gridExtra",warn.conflicts = F)


#função retirada do help(pairs)
panel.hist <- function(x, ...)
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(usr[1:2], 0, 1.5) )
  h <- hist(x, plot = FALSE)
  breaks <- h$breaks; nB <- length(breaks)
  y <- h$counts; y <- y/max(y)
  rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...)
}


#função retirada do help(pairs)
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r <- abs(cor(x, y))
  txt <- format(c(r, 0.123456789), digits = digits)[1]
  txt <- paste0(prefix, txt)
  if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
  text(0.5, 0.5, txt, cex = cex.cor * r)
}

Base de dados

Link da base: Gas Turbine CO and NOx Emission Data Set Data Set

O conjunto de dados contém 36733 instâncias de 11 medidas de sensores agregadas ao longo de uma hora (por meio de média ou soma) de uma turbina a gás localizada na região noroeste da Turquia com o objetivo de estudar as emissões de gases de combustão, ou seja, CO e NOx (NO + NO2). Os dados vêm da mesma usina de energia que o conjunto de dados usado para prever o rendimento de energia líquido por hora. Em contrapartida, esses dados são coletados em outra faixa de dados (01.01.2011 - 31.12.2015), inclui parâmetros da turbina a gás (como temperatura de entrada da turbina e pressão de descarga do compressor) além das variáveis ambientais. Observe que as datas não são fornecidas nas instâncias, mas os dados são classificados em ordem cronológica.

Recomendações do site de onde foi extraída a base de dados: Consulte as informações do atributo e o papel relevante para obter detalhes. Siga o protocolo mencionado no artigo (usando os dados dos primeiros três anos para treinamento / validação cruzada e os dois últimos para teste) para reprodutibilidade e comparabilidade dos trabalhos.

O conjunto de dados vai usado para prever o rendimento de energia da turbina (TEY - variável dependente) usando variáveis ambientais como recursos (variáveis independentes).

As explicações sobre as medições do sensor. Abaixo os nomes das variáveis seguido de sua abreviação e unidade:

Temperatura ambiente (AT) C
Pressão ambiente (AP) mbar
Umidade ambiente (AH) (%)
Pressão diferencial do filtro de ar (AFDP) mbar
Pressão de exaustão da turbina a gás (GTEP) mbar
Temperatura de entrada da turbina (TIT) C
Turbina após temperatura (TAT) C
Pressão de descarga do compressor (CDP) mbar
Rendimento de energia da turbina (TEY) MWH
Monóxido de carbono (CO) mg / m3
Óxidos de nitrogênio (NOx) mg / m3

Abaixo o cabeçalho da base:

setwd("~/UPE/EST_COMPUTACIONAL/gas")

df2011 = fread("gt_2011.csv")
df2012 = fread("gt_2012.csv")
df2013 = fread("gt_2013.csv")
df2014 = fread("gt_2014.csv")
df2015 = fread("gt_2015.csv")

df_treino = rbind(df2011, df2012, df2013)

df_teste = rbind(df2014, df2015)

df_treino %>%  
  head  %>% 
  kable("html", escape = F) %>%
  kable_styling("striped", full_width = F)

AT	AP	AH	AFDP	GTEP	TIT	TAT	TEY	CDP	CO	NOX
4.5878	1018.7	83.675	3.5758	23.979	1086.2	549.83	134.67	11.898	0.32663	81.952
4.2932	1018.3	84.235	3.5709	23.951	1086.1	550.05	134.67	11.892	0.44784	82.377
3.9045	1018.4	84.858	3.5828	23.990	1086.5	550.19	135.10	12.042	0.45144	83.776
3.7436	1018.3	85.434	3.5808	23.911	1086.5	550.17	135.03	11.990	0.23107	82.505
3.7516	1017.8	85.182	3.5781	23.917	1085.9	550.00	134.67	11.910	0.26747	82.028
3.8858	1017.7	83.946	3.5824	23.903	1086.0	549.98	134.67	11.868	0.23473	81.748

Análise Descritiva

Observando a natureza das variáveis e o preenchimento delas.

Variavéis quantitativas contínuas. Abaixo suas estatísticas por ordem alfabética:

df_qt_cont <- df_treino %>%  
  select(AT:NOX) %>%
  gather(var_continuas) 


df_qt_cont %>%
  group_by(var_continuas) %>%
  summarise("Qtd" = n(),
            "Missing" = sum(is.na(value)),
            "Media" = mean(value),
            "DP" = sd(value),
            "Min" = min(value),
            "Q1" = quantile(value, 0.25),
            "Mediana" = quantile(value, 0.50),
            "Q3" = quantile(value, 0.75),
            "Max" = max(value)) %>%
  mutate_each(funs(round(.,1)), -var_continuas) %>% 
  kable("html", escape = F) %>%
  kable_styling("striped", full_width = F)

var_continuas	Qtd	Media	DP	Min	Q1	Mediana	Q3	Max
AFDP	22191	4.0	0.8	2.1	3.4	4.1	4.5	7.6
AH	22191	79.6	13.9	27.5	70.3	82.8	90.5	100.2
AP	22191	1012.8	6.4	985.9	1008.8	1012.4	1016.7	1034.2
AT	22191	17.7	7.4	0.3	11.7	17.7	23.7	34.9
CDP	22191	12.1	1.1	9.9	11.4	12.0	12.4	15.1
CO	22191	2.2	2.3	0.0	1.0	1.5	2.5	44.1
GTEP	22191	25.3	4.2	17.9	22.7	25.0	26.8	37.4
NOX	22191	68.8	11.0	27.8	61.5	67.1	74.6	119.9
TAT	22191	545.5	7.7	512.5	542.6	549.9	550.0	550.6
TEY	22191	133.5	16.0	100.2	124.3	133.8	138.6	174.6
TIT	22191	1083.1	16.8	1000.8	1074.6	1088.1	1095.3	1100.8

Análise de Correlação

Uma breve análise de correlação será feita a partir da correlação de Pearson. Isto é possível pois todas as variáveis são de natureza quantitativa numérica.

A partir das análises de correlação, será decidido quais variáveis entrarão ou não no modelo.

Variáveis independentes que possuem alta correlação entre si, não entrarão as duas no modelo para evitar multicolinearidade (redundância). Assim como variáveis que possuem alta correlação com a variável dependente (TEY), também ficarão de fora do modelo de regressão.

(Valores de referência para correlação de Pearson, segundo Silvia E Shimakura 2006.)

Na ilustração abaixo, já podemos observar que as variáveis TAT, CDP, TIT e GTEP possuem correlação muito forte com a variável dependente TEY.

Também existem diversas variáveis independentes que possuem alta correlação entre elas.

ggcorr(df_treino, label=T)

Gráfico de Dispersão

Para melhor visualização, a base foi dividida em partes pelas variáveis independentes e a variável dependente (TEY).

Nessa primeira ilustração de gráficos de dispersão, observa-se que AT e AP são variáveis candidatas para o modelo, pois possuem fraca correlação com a variável dependente TEY e fraca correlação entre si.

Gráfico de dispersão de TEY, AT, AP:

treino_parte1 = df_treino %>% select("TEY", "AT",   "AP")#, "AH",   "AFDP", "GTEP")
#lattice::splom(treino_parte1)
panel.hist <- function(x, ...)
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r <- abs(cor(x, y))
  txt <- format(c(r, 0.123456789), digits = digits)[1]
  txt <- paste0(prefix, txt)
  if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
  text(0.5, 0.5, txt, cex = cex.cor * r)
}
pairs(treino_parte1, diag.panel = panel.hist, upper.panel = panel.cor)

Nessa segunda ilustração de gráficos de dispersão, observa-se que AFDP e AH são variáveis candidatas a entrar num modelo de regressão para explicar TEY devido a correlação moderada e fraca que eles possuem.

Gráfico de dispersão de TEY, AH, AFDP:

treino_parte2 = df_treino %>% select("TEY","AH","AFDP")#,   "GTEP")
#lattice::splom(treino_parte2)
panel.hist <- function(x, ...)
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r <- abs(cor(x, y))
  txt <- format(c(r, 0.123456789), digits = digits)[1]
  txt <- paste0(prefix, txt)
  if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
  text(0.5, 0.5, txt, cex = cex.cor * r)
}
pairs(treino_parte2, diag.panel = panel.hist, upper.panel = panel.cor)

Nessa terceira ilustração de gráficos de dispersão, observa-se que GTEP, TIT não são variáveis candidatas a entrar num modelo de regressão para explicar TEY devido a correlação muito forte que elas possuem entre si e também com a variável dependente.

Gráfico de dispersão de TEY, GTEP, TIT:

treino_parte3 = df_treino %>% select("TEY", "GTEP", "TIT")
#lattice::splom(treino_parte3)
panel.hist <- function(x, ...)
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r <- abs(cor(x, y))
  txt <- format(c(r, 0.123456789), digits = digits)[1]
  txt <- paste0(prefix, txt)
  if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
  text(0.5, 0.5, txt, cex = cex.cor * r)
}
pairs(treino_parte3, diag.panel = panel.hist, upper.panel = panel.cor)

Nessa quarta ilustração de gráficos de dispersão, observa-se que TAT, CDP não são variáveis candidatas a entrar num modelo de regressão para explicar TEY devido a correlação forte que elas possuem entre si e também com a variável dependente.

Gráfico de dispersão de TEY, TAT, CDP:

treino_parte4 = df_treino %>% select("TEY", "TAT", "CDP") #,    "CO",   "NOX")
#lattice::splom(treino_parte4)
panel.hist <- function(x, ...)
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r <- abs(cor(x, y))
  txt <- format(c(r, 0.123456789), digits = digits)[1]
  txt <- paste0(prefix, txt)
  if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
  text(0.5, 0.5, txt, cex = cex.cor * r)
}
pairs(treino_parte4, diag.panel = panel.hist, upper.panel = panel.cor)

Nessa quinta e última ilustração de gráficos de dispersão, observa-se que apenas NOX é uma variável candidata a entrar num modelo de regressão para explicar TEY devido a correlação fraca que ela possui com a variável dependente. A variável CO possui correlação moderada com a variável dependente, por isso poderá entrar no modelo de regressão.

Gráfico de dispersão de TEY, CO, NOX:

treino_parte5 = df_treino %>% select("TEY", "CO",   "NOX")
#lattice::splom(treino_parte5)
panel.hist <- function(x, ...)
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r <- abs(cor(x, y))
  txt <- format(c(r, 0.123456789), digits = digits)[1]
  txt <- paste0(prefix, txt)
  if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
  text(0.5, 0.5, txt, cex = cex.cor * r)
}
pairs(treino_parte5, diag.panel = panel.hist, upper.panel = panel.cor)

Regressão

No modelo de regressão foram usadas as variáveis AT , AP , AH, NOX, AFDP e CO para explicar TEY.

ajuste_lm <- lm(TEY ~ AT + AP + AH + NOX + AFDP + CO, data = df_treino)

Segue abaixo:

Chamada do modelo;
Medidas resumo dos resíduos;
Tabela de coeficientes, desvios padrão e testes T para hipótese nula de parâmetros iguais a zero;
Média dos quadrados do resíduo e os respectivos graus de liberdade, R^2 ajustado da regressão, Estatística F para qualidade do ajuste (comparação com o modelo com apenas o intercepto).

summary(ajuste_lm)

## 
## Call:
## lm(formula = TEY ~ AT + AP + AH + NOX + AFDP + CO, data = df_treino)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -33.486  -7.202   0.223   5.478  96.989 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 92.534064  11.649460   7.943 2.06e-15 ***
## AT          -0.832036   0.015110 -55.065  < 2e-16 ***
## AP           0.037934   0.011203   3.386  0.00071 ***
## AH          -0.145979   0.006354 -22.973  < 2e-16 ***
## NOX         -0.060375   0.008280  -7.292 3.16e-13 ***
## AFDP         9.683505   0.089481 108.219  < 2e-16 ***
## CO          -2.715593   0.033386 -81.339  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.657 on 22184 degrees of freedom
## Multiple R-squared:  0.637,  Adjusted R-squared:  0.6369 
## F-statistic:  6489 on 6 and 22184 DF,  p-value: < 2.2e-16

Observa-se que todas as variáveis do modelo foram signifcativas.

Unindo na mesma base a variável dependente e a predição, calcula-se o AUC.

predicao_treino = predict.lm(ajuste_lm, df_treino)
predicao_teste = predict.lm(ajuste_lm, df_teste)


predicao_treino = as.data.frame(predicao_treino)
predicao_teste = as.data.frame(predicao_teste)


df_teste_predict = cbind.data.frame(df_teste, predicao_teste)
df_treino_predict = cbind.data.frame(df_treino, predicao_treino)

Para ambas as bases, treinamento e teste, o AUC resultou em cerca de 50%.

BASE DE TREINAMENTO:

# ROC Curve  treino  
roc_treino <- roc.curve(scores.class0 = df_treino_predict$TEY, scores.class1 = df_treino_predict$predicao_treino, curve = T)
plot(roc_treino)

BASE DE TESTE:

# ROC Curve teste   
roc_teste <- roc.curve(scores.class0 = df_teste_predict$TEY, scores.class1 = df_teste_predict$predicao_teste, curve = T)
plot(roc_teste)

Para melhorar o overfitting do modelo, talvez uma validação cruzada resolveria esse problema. Vejamos a seguir:

Validação Cruzada

Ajustando e treinando o modelo com a validação cruzada, temos as seguintes saídas:

ajuste_lm2 <- train(TEY ~ AT + AP + AH + NOX + CO, 
      data = df_treino,
      method = "knn",
      preProcess = c("center", "scale"),
      tuneLength = 10, 
      trControl = trainControl(method="cv", number=10))

ajuste_lm2

## k-Nearest Neighbors 
## 
## 22191 samples
##     5 predictor
## 
## Pre-processing: centered (5), scaled (5) 
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 19972, 19971, 19973, 19972, 19974, 19972, ... 
## Resampling results across tuning parameters:
## 
##   k   RMSE      Rsquared   MAE     
##    5  6.271735  0.8469389  3.656946
##    7  6.298297  0.8456197  3.771124
##    9  6.367264  0.8423062  3.884542
##   11  6.421369  0.8397728  3.979463
##   13  6.478993  0.8370593  4.063056
##   15  6.537991  0.8342359  4.137740
##   17  6.584276  0.8320805  4.201025
##   19  6.638215  0.8294831  4.263267
##   21  6.684270  0.8272352  4.314847
##   23  6.724292  0.8253023  4.363731
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was k = 5.

Resultou em k = 5, o número de folds necessários.

# Realizando a predição para os dados de teste e treinamento

predictions_teste <- predict(ajuste_lm2, df_teste)

predictions_treino <- predict(ajuste_lm2, df_treino)

BASE DE TREINAMENTO:

# ROC Curve  treino  
roc_treino_2 <- roc.curve(scores.class0 = df_treino$TEY, scores.class1 = predictions_treino, curve = T)
plot(roc_treino_2)

BASE DE TESTE:

# ROC Curve teste   
roc_teste_2 <- roc.curve(scores.class0 = df_teste_predict$TEY, scores.class1 = predictions_teste, curve = T)
plot(roc_teste_2)

Conclusão

Como sugerido, foram usados 2011, 2012 e 2013 para treinamento e 2014 e 2015 para teste.

Mesmo com a seleção de variáveis pela correlação de Pearson e utilizando validação cruzada, os resultados não foram satisfatórios devido a um provável compartamento diferente que a base de teste possui em relação a base de treino.

Para solucionar este problema, deve-se usar métodos como Regressão Stepwise para selecionar variáveis onde o modelo aprendesse melhor o comportamento de 2011, 2012 e 2013, afim de refletir uma boa predição para 2014 e 2015.

Visto que estamos falando de anos, o compartamento de treino pode sim ser bem diferente da de teste. Uma outra solução que eu aplicaria, seria juntar todas as bases de 2011 a 2015 e daí separar uma parte para treino e outra para teste. Como por exemplo, 75% da base para treinamento e 25% da base para teste, de forma que mantivessem o mesmo compartamentos gerais da variável dependente.

Regressão

Renata Kelly Marcelino dos Santos

Estatística Computacional

Atividade

Base de dados

Análise Descritiva

Análise de Correlação

Gráfico de Dispersão

Regressão

Validação Cruzada

Conclusão