library(caret)
library(tidyverse)
base = read_csv("base_treino.csv")Exemplo com dados de energia
O problema
A modelagem das cargas de aquecimento e resfriamento em edificações surge como uma ferramenta fundamental para apoiar o projeto de construções mais eficientes do ponto de vista energético. Tradicionalmente, a obtenção dessas medidas depende de simulações físicas complexas, que demandam tempo, custo computacional e conhecimento especializado. Nesse contexto, o uso de métodos de regressão permite estimar essas cargas a partir de características estruturais do edifício de forma rápida, viabilizando a análise de múltiplos cenários ainda na fase de projeto.
A base de dados
A base de dados está disponível em UC Irvine Machine Learning Repository. As variáveis desta base são:
| Variável | Description |
|---|---|
| X1 | Relative Compactness |
| X2 | Surface Area |
| X3 | Wall Area |
| X4 | Roof Area |
| X5 | Overall Height |
| X6 | Orientation |
| Y1 | Heating Load |
| Y2 | Cooling Load |
A base completa foi dividada em base de treino e base de teste, que estão disponíveis no Google Sala de Aula da turma. A base de teste não contém os valores das variáveis de previsão, que são Y1 e Y2. Estas estão em outro arquivo, target_teste.csv.
Primeiro vamos usar somente a base de treino para escolher o modelo de previsão.
Análise Descritiva
Vamos fazer uma análise descritiva de forma a melhor entender as variáveis da base e a relação entre as variáveis independentes (X) e dependentes (Y).
glimpse(base)Rows: 656
Columns: 10
$ X1 <dbl> 0.98, 0.98, 0.98, 0.98, 0.90, 0.90, 0.90, 0.86, 0.86, 0.86, 0.86, 0…
$ X2 <dbl> 514.5, 514.5, 514.5, 514.5, 563.5, 563.5, 563.5, 588.0, 588.0, 588.…
$ X3 <dbl> 294.0, 294.0, 294.0, 294.0, 318.5, 318.5, 318.5, 294.0, 294.0, 294.…
$ X4 <dbl> 110.25, 110.25, 110.25, 110.25, 122.50, 122.50, 122.50, 147.00, 147…
$ X5 <dbl> 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.…
$ X6 <dbl> 2, 3, 4, 5, 2, 3, 5, 2, 3, 4, 5, 2, 3, 5, 2, 3, 4, 5, 2, 4, 5, 2, 3…
$ X7 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ X8 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ Y1 <dbl> 15.55, 15.55, 15.55, 15.55, 20.84, 21.46, 19.68, 19.50, 19.95, 19.3…
$ Y2 <dbl> 21.33, 21.33, 21.33, 21.33, 28.28, 25.38, 29.60, 27.30, 21.97, 23.4…Vamos começar com uma análise descritiva nas variáveis alvo (Y1 e Y2) e depois seguimos para as variáveis dependentes.
Y1 – Heating Load
A variável Y1 indica a quantidade de energia necessária para aquecer o edifício até uma temperatura confortável. Valores maiores indicam maior demanda por aquecimento.
hist(base$Y1, main="Energia necessária para aquecer o edifício")
Vamos criar uma nova base com as variáveis que serão usadas no modelo. essa nova base será nomeada de base_mod.
base_mod = tibble(Y1 = base$Y1)Y2 – Cooling Load
Quantidade de energia necessária para resfriar o edifício. Valores maiores indicam maior necessidade de climatização para remover calor.
hist(base$Y2, main="Energia necessária para esfriar o edifício")
base_mod = base_mod |> cbind(Y2 = base$Y2)X1 – Relative Compactness
A variável X1 - Relative Compactness - mede o quão “compacto” é o edifício, ou seja, a relação entre volume e área externa. Edifícios mais compactos tendem a perder menos calor, sendo mais eficientes energeticamente. Os valores de X1 representam 12 diferentes tipos de construções.
x = base$X1
table(x)x
0.62 0.64 0.66 0.69 0.71 0.74 0.76 0.79 0.82 0.86 0.9 0.98
56 56 54 54 56 55 57 50 52 56 53 57 Os gráficos abaixo mostram como os valores de Y1 e Y2 variam com os valores de X1.
par(mfrow = c(1,2)) # 1 linha, 2 colunas
plot(x, base$Y1, main = "Y1", ylab = ""); abline(v=0.75,lty=2)
plot(x, base$Y2, main = "Y2", ylab = ""); abline(v=0.75,lty=2)
O gráfico mostra que existe uma mudança de comportamento em Y1 e Y2 para valores de X1 menores ou maiores que \(0,75\). Dessa forma, em vez de mantermos a variável X1 na base vamos mudar para uma variável indicadora.
base_mod = base_mod |> cbind(X1 = as.factor(ifelse(base$X1 < 0.75,0,1)))X2 – Surface Area
A variável X2 - Surface Area - indica a área total da superfície externa do edifício (paredes, telhado e piso). Quanto maior essa área, maior a troca de calor com o ambiente externo.
x = base$X2
par(mfrow = c(1,2)) # 1 linha, 2 colunas
plot(x, base$Y1, main = "Y1", ylab = ""); abline(v=670,lty=2)
plot(x, base$Y2, main = "Y2", ylab = ""); abline(v=670,lty=2)
Novamente vou optar por trocar a variável original por uma indicadora. Essa mudança será salva na base_mod.
base_mod = base_mod |> cbind(
X2 = as.factor(ifelse(base$X2 < 670,"pequena","grande")))X3 – Wall Area
A variável X3 – Wall Area - indica a área total das paredes externas. Influencia diretamente as perdas e ganhos de calor por condução. Vamos manter essa variável original na base_mod.
x = base$X3
par(mfrow = c(1,2)) # 1 linha, 2 colunas
plot(x, base$Y1, main = "Y1", ylab = ""); abline(v=670,lty=2)
plot(x, base$Y2, main = "Y2", ylab = ""); abline(v=670,lty=2)
base_mod = base_mod |> cbind(X3 = base$X3)X4 – Roof Area
A variável X4 - Roof Area - indica a área do telhado da edificação. Afeta a quantidade de calor recebida (radiação solar) e perdida pela cobertura.
x = base$X4
table(x)x
110.25 122.5 147 220.5
57 110 158 331 par(mfrow = c(1,2)) # 1 linha, 2 colunas
plot(base$Y1 ~ x, main = "Y1", ylab = ""); abline(v=180,lty=2)
plot(base$Y2 ~ x, main = "Y2", ylab = ""); abline(v=180,lty=2)
Na base só existem 4 valores distintos para X4 e o comportamento de Y1 e Y2 muda muito para o valor mais alto. Vamos criar uma variável indicadora a partir dela.
base_mod = base_mod |> cbind(
X4 = as.factor(ifelse(base$X4 < 180,"pequena","grande")))X5 – Overall Height
A variável X5 - Overall Height - indica a altura total do edifício. Está relacionada ao número de pavimentos e influencia o volume interno e a circulação de ar.
x = base$X5
table(x)x
3.5 7
331 325 par(mfrow = c(1,2)) # 1 linha, 2 colunas
boxplot(base$Y1 ~ x, main = "Y1", ylab = "")
boxplot(base$Y2 ~ x, main = "Y2", ylab = "")
Pelo gráfico percebemos que prédios altos precisam de mais energia, tanto para aquecer quanto para esfriar. Vamos manter essa variável na base, mas de forma qualitativa: pé direito alto ou baixo.
base_mod = base_mod |> cbind(
X5 = as.factor(ifelse(base$X5 < 4,"baixo","alto")))X6 – Orientation
A variável X6 - Orientation - indica a direção para a qual o edifício está orientado (norte, sul, leste ou oeste). Afeta a incidência de radiação solar e, consequentemente, o aquecimento e resfriamento.
x = base$X6
table(x)x
2 3 4 5
164 154 175 163 par(mfrow = c(1,2)) # 1 linha, 2 colunas
boxplot(base$Y1 ~ x, main = "Y1", ylab = "")
boxplot(base$Y2 ~ x, main = "Y2", ylab = "")
Esta variável tem 4 diferentes valores, que signifcam: 2 - Norte, 3 - Leste, 4- Sul e 5 - Oeste. Apesar dela ser numérica na base, veja que ela representa uma variável qualitativa e não quantitativa. Por esse motivo vamos salvá-la como factor.
O gráfico nos mostra que essa não parece ser uma variável tão importante, o comportamento de Y1 e Y2 não muda muito com a variação de X6. Vamos mantê-la na base, mas pode ser uma opção retirá-la.
base_mod = base_mod |> cbind(
X6 = as.factor(
recode_values(
base$X6,
2 ~ "Norte",
3 ~"Leste",
4 ~"Sul",
5 ~"Oeste")))X7 – Glazing Area
A variável X7 - Glazing Area - indica a proporção da área do edifício composta por superfícies envidraçadas (janelas). Maior área de vidro aumenta a entrada de luz e calor, impactando o consumo energético.
x = base$X7
table(x)x
0 0.1 0.25 0.4
42 205 206 203 par(mfrow = c(1,2)) # 1 linha, 2 colunas
boxplot(base$Y1 ~ x, main = "Y1", ylab = "")
boxplot(base$Y2 ~ x, main = "Y2", ylab = "")
Vamos mater a variável X7 na base com seus valores originais.
base_mod = base_mod |> cbind(X7 = base$X7)X8 – Glazing Area Distribution
A variável X8 – Glazing Area Distribution - idnica a distribuição das janelas nas fachadas do edifício. Indica em quais lados (orientações) o vidro está concentrado, influenciando o ganho solar ao longo do dia.
#X8 - Glazing area distribution
x = base$X8
table(x)x
0 1 2 3 4 5
42 122 127 111 124 130 x = as.factor(x)
boxplot(base$Y1~x, main = "Y1", ylab = "")
boxplot(base$Y2~x, main = "Y1", ylab = "")
Os valores de X8 representam: 0, Sem janelas (ou glazing area = 0); 1 - Vidro distribuído uniformemente em todas as fachadas; 2 - Vidro concentrado na fachada norte; 3 - Vidro concentrado na fachada leste; 4 - Vidro concentrado na fachada sul; 5 - Vidro concentrado na fachada oeste. Com isso, mais uma vez, percebemos que apesar desta variável está guardando um valor numérico, trata-se de uma variável qualitativa. Vamos salvá-la desta forma na base.
Repare que o mais relevante neste caso é a categoria sem janelas para as demais. Mas, a categoria sem janela já é identificada como X7=0. Então, uma possibilidade a ser pensada é a retirada desta variável da base ou a sua modificação para duas categorias: “Sem Vidros” e “Com Vidros”. Isso veremos depois, nos ajustes dos modelos.
base_mod = base_mod |> cbind(X8 = as.factor(
recode_values(
base$X8,
0 ~ "SemVidros",
1 ~ "Uniforme",
2 ~ "Norte",
3 ~ "Leste",
4 ~ "Sul",
5 ~ "Oeste"))
)Com isso chegamos a seguinte base de trabalho:
glimpse(base_mod)Rows: 656
Columns: 10
$ Y1 <dbl> 15.55, 15.55, 15.55, 15.55, 20.84, 21.46, 19.68, 19.50, 19.95, 19.3…
$ Y2 <dbl> 21.33, 21.33, 21.33, 21.33, 28.28, 25.38, 29.60, 27.30, 21.97, 23.4…
$ X1 <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0…
$ X2 <fct> pequena, pequena, pequena, pequena, pequena, pequena, pequena, pequ…
$ X3 <dbl> 294.0, 294.0, 294.0, 294.0, 318.5, 318.5, 318.5, 294.0, 294.0, 294.…
$ X4 <fct> pequena, pequena, pequena, pequena, pequena, pequena, pequena, pequ…
$ X5 <fct> alto, alto, alto, alto, alto, alto, alto, alto, alto, alto, alto, a…
$ X6 <fct> Norte, Leste, Sul, Oeste, Norte, Leste, Oeste, Norte, Leste, Sul, O…
$ X7 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ X8 <fct> SemVidros, SemVidros, SemVidros, SemVidros, SemVidros, SemVidros, S…Vou salvar essa base para não precisar rodar o código novamente e poder recomeçar deste ponto da análise.
saveRDS(base_mod,"base_treino_mod.RDS")Separação da base em treino e validação
Para que possamos analisar o desempenho do modelo de regressão com dados diferentes dos de treinamento, vamos separar a base em treino e validação.
N = nrow(base_mod)
indices = createDataPartition(1:N,p = 0.75)
base_treino = base_mod[-indices[[1]],]
base_valida = base_mod[indices[[1]],]Tratamento das variáveis
Antes de ajustar o modelo de regressão precisamos tratar as variáveis. As variáveis quantitativas devem ser padronizadas e as variáveis qualitativas transformadas em binárias. Para facilitar, vamos fazer a transformação na base de treino e também já na base de teste.
Primeiro a padronização das variáveis dependentes.
m3 = mean(base_treino$X3)
s3 = sd(base_treino$X3)
base_treino$X3 = (base_treino$X3 - m3)/s3
base_valida$X3 = (base_valida$X3 - m3)/s3 m7 = mean(base_treino$X7)
s7 = sd(base_treino$X7)
base_treino$X7 = (base_treino$X7 - m7)/s7
base_valida$X7 = (base_valida$X7 - m7)/s7Agora a padronização das variáveis respostas, Y1 e Y2, considerando que a função de ativação usada será a logística. Como estamos substituíndo os valores na base, e não criando outra base, vamos salvar os valores originais das variáveis alvo para futura comparação dos resultados.
max1 = max(base_treino$Y1)
min1 = min(base_treino$Y1)
Y1_treino = base_treino$Y1
base_treino$Y1 = (base_treino$Y1 - min1)/(max1 - min1)
Y1_valida = base_valida$Y1
base_valida$Y1 = (base_valida$Y1 - min1)/(max1 - min1) max2 = max(base_treino$Y2)
min2 = min(base_treino$Y2)
Y2_treino = base_treino$Y2
base_treino$Y2 = (base_treino$Y2 - min2)/(max2 - min2)
Y2_valida = base_valida$Y2
base_valida$Y2 = (base_valida$Y2 - min2)/(max2 - min2)Por fim vamos usar a função model.marix para criar as variáveis indicadores no tratamento das qualitaivas.
mat_treino = model.matrix(~. , data = base_treino)
mat_valida = model.matrix(~. , data = base_valida)
mat_treino = mat_treino[,-1]
mat_valida = mat_valida[,-1]Modelos de Regressão para Y1
Vamos agora ajustar diferentes modelos de regressão RNA e comparar os desempenhos, tanto na base de treino quanto na base de teste.
Como a nossa base tem 2 variáveis alvo. Aqui vou analisar um modelo para a variável Y1 e deixo como exercício. Aviso importante: precisamos garantir que os valores Y2 não sejam considerados no ajuste do modelo de Y1 e que os valores Y1 não sejam considerados no ajuste do modelo de Y2.
O primeiro modelo será completo (com todas as variáveis dependentes) e com 1 camada oculta com 1 neurônio, chamado de NN1.
library(neuralnet)
NN1_Y1 = neuralnet(
Y1 ~ . ,
hidden = 1,
data = mat_treino[,-2],
lifesign = 'full')Vamos ajustar mais dois modelos com todas as variáveis e 1 camada oculta, NN2 e NN3, com 2 e 3 neruônios na camada oculta, respectivamente.
NN2_Y1 = neuralnet(
Y1 ~ . ,
hidden = 2,
data = mat_treino[,-2],
lifesign = 'full') NN3_Y1 = neuralnet(
Y1 ~ . ,
hidden = 3,
data = mat_treino[,-2],
lifesign = 'full')Vamos replicar os mesmos modelos, agora sem a variável X6 - Orientation. Isso porque na análise descritiva ela não se mostrou muito importante.
colnames(mat_treino) [1] "Y1" "Y2" "X11" "X2pequena" "X3"
[6] "X4pequena" "X5baixo" "X6Norte" "X6Oeste" "X6Sul"
[11] "X7" "X8Norte" "X8Oeste" "X8SemVidros" "X8Sul"
[16] "X8Uniforme" Para isso vamos retirar as colunas 8, 9 e 10 da mat_treino.
NN1_Y1_semX6 = neuralnet(
Y1 ~ . ,
hidden = 1,
data = mat_treino[,-c(2,8,9,10)],
lifesign = 'full') NN2_Y1_semX6 = neuralnet(
Y1 ~ . ,
hidden = 2,
data = mat_treino[,-c(2,8,9,10)],
lifesign = 'full') NN3_Y1_semX6 = neuralnet(
Y1 ~ . ,
hidden = 3,
data = mat_treino[,-c(2,8,9,10)],
lifesign = 'full')Para termos algo com mais camadas, vamos fazer um RNA com 2 camadas ocultas, 3 neurônios na primeira e 2 na segunda, com e sem X6.
NN22_Y1 = neuralnet(
Y1 ~ . ,
hidden = c(2,2),
data = mat_treino[,-2],
lifesign = 'full') NN22_Y1_semX6 = neuralnet(
Y1 ~ . ,
hidden = c(2,2),
data = mat_treino[,-c(2,8,9,10)],
lifesign = 'full')E uma com 2 camadas ocultas, 3 neurônios na primeira e 2 na segunda, com e sem X6.
NN32_Y1 = neuralnet(
Y1 ~ . ,
hidden = c(3,2),
data = mat_treino[,-2],
lifesign = 'full') NN32_Y1_semX6 = neuralnet(
Y1 ~ . ,
hidden = c(3,2),
data = mat_treino[,-c(2,8,9,10)],
lifesign = 'full')Vamos repetir todas as arquiteturas agora considerando o modelo sem X6 e com X8 binário: sem vidro ou outros casos. Para isso vamos manter “X8SemVidros” (coluna 14) e tirar as outras colunas da variável X8: 12, 13, 15 e 16.
colnames(mat_treino) [1] "Y1" "Y2" "X11" "X2pequena" "X3"
[6] "X4pequena" "X5baixo" "X6Norte" "X6Oeste" "X6Sul"
[11] "X7" "X8Norte" "X8Oeste" "X8SemVidros" "X8Sul"
[16] "X8Uniforme" NN1_Y1_semX6X8 = neuralnet(
Y1 ~ . ,
hidden = 1,
data = mat_treino[,-c(2,8,9,10)],
lifesign = 'full')NN2_Y1_semX6X8 = neuralnet(
Y1 ~ . ,
hidden = 2,
data = mat_treino[,-c(2,8,9,10)],
lifesign = 'full')NN3_Y1_semX6X8 = neuralnet(
Y1 ~ . ,
hidden = 3,
data = mat_treino[,-c(2,8,9,10)],
lifesign = 'full')NN22_Y1_semX6X8 = neuralnet(
Y1 ~ . ,
hidden = c(2,2),
data = mat_treino[,-c(2,8,9,10)],
lifesign = 'full')NN32_Y1_semX6X8 = neuralnet(
Y1 ~ . ,
hidden = c(3,2),
data = mat_treino[,-c(2,8,9,10)],
lifesign = 'full')Medidas de erro (R\(^2\))
Modelos completos
Primeiro para Y1 no modelo NN1_Y1.
y1_Y1_treino = NN1_Y1$net.result[[1]][,1]
y1_Y1_treino = y1_Y1_treino*(max1 - min1) + min1
y1_Y1_valida = predict(NN1_Y1,newdata = mat_valida)
y1_Y1_valida = y1_Y1_valida[,1]
y1_Y1_valida = y1_Y1_valida*(max1 - min1) + min1
sse_treino_1_Y1 = sum((y1_Y1_treino - Y1_treino)^2)
sse_valida_1_Y1 = sum((y1_Y1_valida - Y1_valida)^2)
sst_treino_1 = sum((mean(Y1_treino) - Y1_treino)^2)
sst_valida_1 = sum((mean(Y1_treino) - Y1_valida)^2)
(R2_treino_1_Y1 = 1 - sse_treino_1_Y1/sst_treino_1)[1] 0.9460653(R2_valida_1_Y1 = 1 - sse_valida_1_Y1/sst_valida_1)[1] 0.9215985Y1 no modelo NN2_Y1.
y2_Y1_treino = NN2_Y1$net.result[[1]][,1]
y2_Y1_treino = y2_Y1_treino*(max1 - min1) + min1
y2_Y1_valida = predict(NN2_Y1,newdata = mat_valida)
y2_Y1_valida = y2_Y1_valida[,1]
y2_Y1_valida = y2_Y1_valida*(max1 - min1) + min1
sse_treino_2_Y1 = sum((y2_Y1_treino - Y1_treino)^2)
sse_valida_2_Y1 = sum((y2_Y1_valida - Y1_valida)^2)
(R2_treino_2_Y1 = 1 - sse_treino_2_Y1/sst_treino_1)[1] 0.9505097(R2_valida_2_Y1 = 1 - sse_valida_2_Y1/sst_valida_1)[1] 0.9183479Y1 no modelo NN3_Y1.
y3_Y1_treino = NN3_Y1$net.result[[1]][,1]
y3_Y1_treino = y3_Y1_treino*(max1 - min1) + min1
y3_Y1_valida = predict(NN3_Y1,newdata = mat_valida)
y3_Y1_valida = y3_Y1_valida[,1]
y3_Y1_valida = y3_Y1_valida*(max1 - min1) + min1
sse_treino_3_Y1 = sum((y3_Y1_treino - Y1_treino)^2)
sse_valida_3_Y1 = sum((y3_Y1_valida - Y1_valida)^2)
(R2_treino_3_Y1 = 1 - sse_treino_3_Y1/sst_treino_1)[1] 0.9640016(R2_valida_3_Y1 = 1 - sse_valida_3_Y1/sst_valida_1)[1] 0.8881685Y1 no modelo NN22_Y1.
y22_Y1_treino = NN22_Y1$net.result[[1]][,1]
y22_Y1_treino = y22_Y1_treino*(max1 - min1) + min1
y22_Y1_valida = predict(NN22_Y1,newdata = mat_valida)
y22_Y1_valida = y22_Y1_valida[,1]
y22_Y1_valida = y22_Y1_valida*(max1 - min1) + min1
sse_treino_22_Y1 = sum((y22_Y1_treino - Y1_treino)^2)
sse_valida_22_Y1 = sum((y22_Y1_valida - Y1_valida)^2)
(R2_treino_22_Y1 = 1 - sse_treino_22_Y1/sst_treino_1)[1] 0.9581375(R2_valida_22_Y1 = 1 - sse_valida_22_Y1/sst_valida_1)[1] 0.8299462Y1 no modelo NN32_Y1.
y32_Y1_treino = NN32_Y1$net.result[[1]][,1]
y32_Y1_treino = y32_Y1_treino*(max1 - min1) + min1
y32_Y1_valida = predict(NN32_Y1,newdata = mat_valida)
y32_Y1_valida = y32_Y1_valida[,1]
y32_Y1_valida = y32_Y1_valida*(max1 - min1) + min1
sse_treino_32_Y1 = sum((y32_Y1_treino - Y1_treino)^2)
sse_valida_32_Y1 = sum((y32_Y1_valida - Y1_valida)^2)
(R2_treino_32_Y1 = 1 - sse_treino_32_Y1/sst_treino_1)[1] 0.9728928(R2_valida_32_Y1 = 1 - sse_valida_32_Y1/sst_valida_1)[1] 0.9256612Modelos sem X6
Y1 no modelo NN1_Y1_semX6.
y1_Y1_semX6_treino = NN1_Y1_semX6$net.result[[1]][,1]
y1_Y1_semX6_treino = y1_Y1_semX6_treino*(max1 - min1) + min1
y1_Y1_semX6_valida = predict(NN1_Y1_semX6,newdata = mat_valida)
y1_Y1_semX6_valida = y1_Y1_semX6_valida[,1]
y1_Y1_semX6_valida = y1_Y1_semX6_valida*(max1 - min1) + min1
sse_treino_1_Y1_semX6 = sum((y1_Y1_semX6_treino - Y1_treino)^2)
sse_valida_1_Y1_semX6 = sum((y1_Y1_semX6_valida - Y1_valida)^2)
(R2_treino_1_Y1_semX6 = 1 - sse_treino_1_Y1_semX6/sst_treino_1)[1] 0.9439288(R2_valida_1_Y1_semX6 = 1 - sse_valida_1_Y1_semX6/sst_valida_1)[1] 0.9250936Y1 no modelo NN2_Y1_semX6.
y2_Y1_semX6_treino = NN2_Y1_semX6$net.result[[1]][,1]
y2_Y1_semX6_treino = y2_Y1_semX6_treino*(max1 - min1) + min1
y2_Y1_semX6_valida = predict(NN2_Y1_semX6,newdata = mat_valida)
y2_Y1_semX6_valida = y2_Y1_semX6_valida[,1]
y2_Y1_semX6_valida = y2_Y1_semX6_valida*(max1 - min1) + min1
sse_treino_2_Y1_semX6 = sum((y2_Y1_semX6_treino - Y1_treino)^2)
sse_valida_2_Y1_semX6 = sum((y2_Y1_semX6_valida - Y1_valida)^2)
(R2_treino_2_Y1_semX6 = 1 - sse_treino_2_Y1_semX6/sst_treino_1)[1] 0.9534322(R2_valida_2_Y1_semX6 = 1 - sse_valida_2_Y1_semX6/sst_valida_1)[1] 0.925121Y1 no modelo NN3_Y1_semX6.
y3_Y1_semX6_treino = NN3_Y1_semX6$net.result[[1]][,1]
y3_Y1_semX6_treino = y3_Y1_semX6_treino*(max1 - min1) + min1
y3_Y1_semX6_valida = predict(NN3_Y1_semX6,newdata = mat_valida)
y3_Y1_semX6_valida = y3_Y1_semX6_valida[,1]
y3_Y1_semX6_valida = y3_Y1_semX6_valida*(max1 - min1) + min1
sse_treino_3_Y1_semX6 = sum((y3_Y1_semX6_treino - Y1_treino)^2)
sse_valida_3_Y1_semX6 = sum((y3_Y1_semX6_valida - Y1_valida)^2)
(R2_treino_3_Y1_semX6 = 1 - sse_treino_3_Y1_semX6/sst_treino_1)[1] 0.9572166(R2_valida_3_Y1_semX6 = 1 - sse_valida_3_Y1_semX6/sst_valida_1)[1] 0.9271669Y1 no modelo NN22_Y1_semX6.
y22_Y1_semX6_treino = NN22_Y1_semX6$net.result[[1]][,1]
y22_Y1_semX6_treino = y22_Y1_semX6_treino*(max1 - min1) + min1
y22_Y1_semX6_valida = predict(NN22_Y1_semX6,newdata = mat_valida)
y22_Y1_semX6_valida = y22_Y1_semX6_valida[,1]
y22_Y1_semX6_valida = y22_Y1_semX6_valida*(max1 - min1) + min1
sse_treino_22_Y1_semX6 = sum((y22_Y1_semX6_treino - Y1_treino)^2)
sse_valida_22_Y1_semX6 = sum((y22_Y1_semX6_valida - Y1_valida)^2)
(R2_treino_22_Y1_semX6 = 1 - sse_treino_22_Y1_semX6/sst_treino_1)[1] 0.9526964(R2_valida_22_Y1_semX6 = 1 - sse_valida_22_Y1_semX6/sst_valida_1)[1] 0.9313371Y1 no modelo NN32_Y1_semX6.
y32_Y1_semX6_treino = NN32_Y1_semX6$net.result[[1]][,1]
y32_Y1_semX6_treino = y32_Y1_semX6_treino*(max1 - min1) + min1
y32_Y1_semX6_valida = predict(NN32_Y1_semX6,newdata = mat_valida)
y32_Y1_semX6_valida = y32_Y1_semX6_valida[,1]
y32_Y1_semX6_valida = y32_Y1_semX6_valida*(max1 - min1) + min1
sse_treino_32_Y1_semX6 = sum((y32_Y1_semX6_treino - Y1_treino)^2)
sse_valida_32_Y1_semX6 = sum((y32_Y1_semX6_valida - Y1_valida)^2)
(R2_treino_32_Y1_semX6 = 1 - sse_treino_32_Y1_semX6/sst_treino_1)[1] 0.9724241(R2_valida_32_Y1_semX6 = 1 - sse_valida_32_Y1_semX6/sst_valida_1)[1] 0.9584904Modelos sem X6 e X8 modificada
Y1 no modelo NN1_Y1_semX6X8.
y1_Y1_semX6X8_treino = NN1_Y1_semX6X8$net.result[[1]][,1]
y1_Y1_semX6X8_treino = y1_Y1_semX6X8_treino*(max1 - min1) + min1
y1_Y1_semX6X8_valida = predict(NN1_Y1_semX6X8,newdata = mat_valida)
y1_Y1_semX6X8_valida = y1_Y1_semX6X8_valida[,1]
y1_Y1_semX6X8_valida = y1_Y1_semX6X8_valida*(max1 - min1) + min1
sse_treino_1_Y1_semX6X8 = sum((y1_Y1_semX6X8_treino - Y1_treino)^2)
sse_valida_1_Y1_semX6X8 = sum((y1_Y1_semX6X8_valida - Y1_valida)^2)
(R2_treino_1_Y1_semX6X8 = 1 - sse_treino_1_Y1_semX6X8/sst_treino_1)[1] 0.943223(R2_valida_1_Y1_semX6X8 = 1 - sse_valida_1_Y1_semX6X8/sst_valida_1)[1] 0.9252937Y1 no modelo NN2_Y1_semX6X8.
y2_Y1_semX6X8_treino = NN2_Y1_semX6X8$net.result[[1]][,1]
y2_Y1_semX6X8_treino = y2_Y1_semX6X8_treino*(max1 - min1) + min1
y2_Y1_semX6X8_valida = predict(NN2_Y1_semX6X8,newdata = mat_valida)
y2_Y1_semX6X8_valida = y2_Y1_semX6X8_valida[,1]
y2_Y1_semX6X8_valida = y2_Y1_semX6X8_valida*(max1 - min1) + min1
sse_treino_2_Y1_semX6X8 = sum((y2_Y1_semX6X8_treino - Y1_treino)^2)
sse_valida_2_Y1_semX6X8 = sum((y2_Y1_semX6X8_valida - Y1_valida)^2)
(R2_treino_2_Y1_semX6X8 = 1 - sse_treino_2_Y1_semX6X8/sst_treino_1)[1] 0.9495505(R2_valida_2_Y1_semX6X8 = 1 - sse_valida_2_Y1_semX6X8/sst_valida_1)[1] 0.9150365Y1 no modelo NN3_Y1_semX6X8.
y3_Y1_semX6X8_treino = NN3_Y1_semX6X8$net.result[[1]][,1]
y3_Y1_semX6X8_treino = y3_Y1_semX6X8_treino*(max1 - min1) + min1
y3_Y1_semX6X8_valida = predict(NN3_Y1_semX6X8,newdata = mat_valida)
y3_Y1_semX6X8_valida = y3_Y1_semX6X8_valida[,1]
y3_Y1_semX6X8_valida = y3_Y1_semX6X8_valida*(max1 - min1) + min1
sse_treino_3_Y1_semX6X8 = sum((y3_Y1_semX6X8_treino - Y1_treino)^2)
sse_valida_3_Y1_semX6X8 = sum((y3_Y1_semX6X8_valida - Y1_valida)^2)
(R2_treino_3_Y1_semX6X8 = 1 - sse_treino_3_Y1_semX6X8/sst_treino_1)[1] 0.9788383(R2_valida_3_Y1_semX6X8 = 1 - sse_valida_3_Y1_semX6X8/sst_valida_1)[1] 0.9594478Y1 no modelo NN22_Y1_semX6X8.
y22_Y1_semX6X8_treino = NN22_Y1_semX6X8$net.result[[1]][,1]
y22_Y1_semX6X8_treino = y22_Y1_semX6X8_treino*(max1 - min1) + min1
y22_Y1_semX6X8_valida = predict(NN22_Y1_semX6X8,newdata = mat_valida)
y22_Y1_semX6X8_valida = y22_Y1_semX6X8_valida[,1]
y22_Y1_semX6X8_valida = y22_Y1_semX6X8_valida*(max1 - min1) + min1
sse_treino_22_Y1_semX6X8 = sum((y22_Y1_semX6X8_treino - Y1_treino)^2)
sse_valida_22_Y1_semX6X8 = sum((y22_Y1_semX6X8_valida - Y1_valida)^2)
(R2_treino_22_Y1_semX6X8 = 1 - sse_treino_22_Y1_semX6X8/sst_treino_1)[1] 0.9536754(R2_valida_22_Y1_semX6X8 = 1 - sse_valida_22_Y1_semX6X8/sst_valida_1)[1] 0.9304224Y1 no modelo NN32_Y1_semX6X8.
y32_Y1_semX6X8_treino = NN32_Y1_semX6X8$net.result[[1]][,1]
y32_Y1_semX6X8_treino = y32_Y1_semX6X8_treino*(max1 - min1) + min1
y32_Y1_semX6X8_valida = predict(NN32_Y1_semX6X8,newdata = mat_valida)
y32_Y1_semX6X8_valida = y32_Y1_semX6X8_valida[,1]
y32_Y1_semX6X8_valida = y32_Y1_semX6X8_valida*(max1 - min1) + min1
sse_treino_32_Y1_semX6X8 = sum((y32_Y1_semX6X8_treino - Y1_treino)^2)
sse_valida_32_Y1_semX6X8 = sum((y32_Y1_semX6X8_valida - Y1_valida)^2)
(R2_treino_32_Y1_semX6X8 = 1 - sse_treino_32_Y1_semX6X8/sst_treino_1)[1] 0.9776626(R2_valida_32_Y1_semX6X8 = 1 - sse_valida_32_Y1_semX6X8/sst_valida_1)[1] 0.9562661Comparação entre modelos
Modelos completos com 1 camada oculta
Primeiro vamos comparar os modelos com 1 camada oculta e completos.
mat_R2 = matrix(c(
R2_treino_1_Y1,R2_valida_1_Y1,
R2_treino_2_Y1,R2_valida_2_Y1,
R2_treino_3_Y1,R2_valida_3_Y1
),byrow = T,ncol=2)
colnames(mat_R2) = c("Treino","Valida")
rownames(mat_R2) =
c("1","2","3")
mat_R2 Treino Valida
1 0.9460653 0.9215985
2 0.9505097 0.9183479
3 0.9640016 0.8881685barplot(t(mat_R2),beside = T,ylim = c(0,1))
abline(h=max(mat_R2[,1]),lty=2)
abline(h=max(mat_R2[,2]),lty=2)
O modelo com 3 neurônios superou os demais na base de treino mas teve o pior resultado na base de validação. O modelo com 2 neurônios foi ligeramente superios ao com 1 neurônio na base de validação.
Modelos sem X6 com 1 camada oculta
mat_R2 = matrix(c(
R2_treino_1_Y1_semX6,R2_valida_1_Y1_semX6,
R2_treino_2_Y1_semX6,R2_valida_2_Y1_semX6,
R2_treino_3_Y1_semX6,R2_valida_3_Y1_semX6
),byrow = T,ncol=2)
colnames(mat_R2) = c("Treino","Valida")
rownames(mat_R2) =
c("1","2","3")
mat_R2 Treino Valida
1 0.9439288 0.9250936
2 0.9534322 0.9251210
3 0.9572166 0.9271669barplot(t(mat_R2),beside = T,ylim = c(0,1))
abline(h=max(mat_R2[,1]),lty=2)
abline(h=max(mat_R2[,2]),lty=2)
O sobreajuste no modelo com 3 neurônios não foi mais observado nesse caso, sem a variável X6. Mas o desempenho dele foi praticamente igual ao com 2 neurônio e pouco superior ao com 1 neurônio.
Modelos sem X6, com X8 modificado e com 1 camada oculta
mat_R2 = matrix(c(
R2_treino_1_Y1_semX6X8,R2_valida_1_Y1_semX6X8,
R2_treino_2_Y1_semX6X8,R2_valida_2_Y1_semX6X8,
R2_treino_3_Y1_semX6X8,R2_valida_3_Y1_semX6X8
),byrow = T,ncol=2)
colnames(mat_R2) = c("Treino","Valida")
rownames(mat_R2) =
c("1","2","3")
mat_R2 Treino Valida
1 0.9432230 0.9252937
2 0.9495505 0.9150365
3 0.9788383 0.9594478barplot(t(mat_R2),beside = T,ylim = c(0,1))
abline(h=max(mat_R2[,1]),lty=2)
abline(h=max(mat_R2[,2]),lty=2)
Nesse caso o modelo com 3 neurônios apresentou resultados superiores aos demais.
Melhores modelos com 1 camada oculta e diferentes variáveis
mat_R2 = matrix(c(
R2_treino_2_Y1,R2_valida_2_Y1,
R2_treino_2_Y1_semX6,R2_valida_2_Y1_semX6,
R2_treino_3_Y1_semX6X8,R2_valida_3_Y1_semX6X8
),byrow = T,ncol=2)
colnames(mat_R2) = c("Treino","Valida")
rownames(mat_R2) =
c("2","2_semX6","3_semX6X8")
mat_R2 Treino Valida
2 0.9505097 0.9183479
2_semX6 0.9534322 0.9251210
3_semX6X8 0.9788383 0.9594478barplot(t(mat_R2),beside = T,ylim = c(0,1))
abline(h=max(mat_R2[,1]),lty=2)
abline(h=max(mat_R2[,2]),lty=2)
Entre esses 3 o modelo com 3 neurônios na única camada oculta, sem X6 e com X8 modificado teve melhor desempenho.
Comparação Final
Vamos comparar o modelo sem X6, com X8 modificado e com 3 neurônios em 1 camada oculta com os outros de 2 camadas ocultas e 2 neurônios em cada.
mat_R2 = matrix(c(
R2_treino_3_Y1_semX6X8,R2_valida_3_Y1_semX6X8,
R2_treino_22_Y1,R2_valida_22_Y1,
R2_treino_22_Y1_semX6,R2_valida_22_Y1_semX6,
R2_treino_22_Y1_semX6X8,R2_valida_22_Y1_semX6X8
),byrow = T,ncol=2)
colnames(mat_R2) = c("Treino","Valida")
rownames(mat_R2) =
c("3_semX6X8","22","22_semX6","22_semX6X8")
mat_R2 Treino Valida
3_semX6X8 0.9788383 0.9594478
22 0.9581375 0.8299462
22_semX6 0.9526964 0.9313371
22_semX6X8 0.9536754 0.9304224barplot(t(mat_R2),beside = T,ylim = c(0,1))
abline(h=max(mat_R2[,1]),lty=2)
abline(h=max(mat_R2[,2]),lty=2)
O modelo com 1 camada oculta, sem X6 e com X8 modificado pareceu melhor que os demais. Vamos compará-lo com os RNA de 2 camadas ocultas com 3 e 2 neurônios.
mat_R2 = matrix(c(
R2_treino_3_Y1_semX6X8,R2_valida_3_Y1_semX6X8,
R2_treino_32_Y1,R2_valida_32_Y1,
R2_treino_32_Y1_semX6,R2_valida_32_Y1_semX6,
R2_treino_32_Y1_semX6X8,R2_valida_32_Y1_semX6X8
),byrow = T,ncol=2)
colnames(mat_R2) = c("Treino","Valida")
rownames(mat_R2) =
c("3_semX6X8","32","32_semX6","32_semX6X8")
mat_R2 Treino Valida
3_semX6X8 0.9788383 0.9594478
32 0.9728928 0.9256612
32_semX6 0.9724241 0.9584904
32_semX6X8 0.9776626 0.9562661barplot(t(mat_R2),beside = T,ylim = c(0,1))
abline(h=max(mat_R2[,1]),lty=2)
abline(h=max(mat_R2[,2]),lty=2)
Ficaremos com esse modelo: 3 neurônios em 1 camda oculta, sem X6 e com X8 modificado.
Modelo Final com a base de treino completa
Agora que usamos a base de validação para escolher a melhor arquitetura e realizar a seleção das variáveis, vamos ajustá-lo com a base de treino completa.
Primeiro retirando a variável X6 e modificando a X8.
base_mod = base_mod |> select(-X6)
base_mod$X8 = as.factor(ifelse(base_mod$X8 == "SemVidros","SemVidros","Outro"))Padronizando as variáveis quantitativas e criando as binárias para as qualitativas.
m3 = mean(base_mod$X3)
s3 = sd(base_mod$X3)
base_mod$X3 = (base_mod$X3 - m3)/s3
m7 = mean(base_mod$X7)
s7 = sd(base_mod$X7)
base_mod$X7 = (base_mod$X7 - m7)/s7
Y1_treino = base_mod$Y1
max1 = max(Y1_treino)
min1 = min(Y1_treino)
base_mod$Y1 = (base_mod$Y1 - min1)/(max1 - min1)
mat = model.matrix(~. , data = base_mod)
mat = mat[,-1]
head(mat) Y1 Y2 X11 X2pequena X3 X4pequena X5baixo X7
1 0.2572122 21.33 1 1 -0.55649019 1 0 -1.752823
2 0.2572122 21.33 1 1 -0.55649019 1 0 -1.752823
3 0.2572122 21.33 1 1 -0.55649019 1 0 -1.752823
4 0.2572122 21.33 1 1 -0.55649019 1 0 -1.752823
5 0.3998382 28.28 1 1 -0.00169146 1 0 -1.752823
6 0.4165543 25.38 1 1 -0.00169146 1 0 -1.752823
X8SemVidros
1 1
2 1
3 1
4 1
5 1
6 1Treinamento do modelo escolhido.
NN_final = neuralnet(
Y1 ~ . ,
hidden = 3,
data = mat[,-2],
lifesign = "full")
plot(NN_final,rep="best")
Previsão na base de teste
Leitura da base de teste.
base_teste = read_csv("base_teste.csv")Organização da base de teste, da mesma forma que foi feito com a base de treino.
base_mod_teste = tibble(X1 = as.factor(ifelse(base_teste$X1 < 0.75,0,1)))
base_mod_teste = base_mod_teste |> cbind(
X2 = as.factor(ifelse(base_teste$X2 < 670,"pequena","grande")))
base_mod_teste = base_mod_teste |> cbind(X3 = base_teste$X3)
base_mod_teste = base_mod_teste |> cbind(
X4 = as.factor(ifelse(base_teste$X4 < 180,"pequena","grande")))
base_mod_teste = base_mod_teste |> cbind(
X5 = as.factor(ifelse(base_teste$X5 < 4,"baixo","alto")))
base_mod_teste = base_mod_teste |> cbind(X7 = base_teste$X7)
base_mod_teste = base_mod_teste |> cbind(
X8 = as.factor(ifelse(base_teste$X8 == 0,"SemVidros","Outro")
))Padronização das variáveis dependentes quantitativas.
base_mod_teste$X3 = (base_mod_teste$X3 - m3)/s3
base_mod_teste$X7 = (base_mod_teste$X7 - m7)/s7Tratamento nas variáveis categóricas.
mat_teste = model.matrix(~. , base_mod_teste)
mat_teste = mat_teste[,-1]
colnames(mat_teste)[1] "X11" "X2pequena" "X3" "X4pequena" "X5baixo"
[6] "X7" "X8SemVidros"colnames(mat)[1] "Y1" "Y2" "X11" "X2pequena" "X3"
[6] "X4pequena" "X5baixo" "X7" "X8SemVidros"Previsões com o modelo final.
y1_Y1_treino = NN_final$net.result[[1]][,1]
y1_Y1_treino = y1_Y1_treino*(max1 - min1) + min1
y1_Y1_teste = predict(NN_final,newdata = mat_teste)
y1_Y1_teste = y1_Y1_teste[,1]
y1_Y1_teste = y1_Y1_teste*(max1 - min1) + min1Leitura dos valores de Y1 da base de teste para cálculo do R\(^2\).
target_teste = read_csv("target_teste.csv")
Y1_teste = target_teste$Y1Cálculo do R\(^2\) final na base de treino e de teste.
sse_treino_1_Y1 = sum((y1_Y1_treino - Y1_treino)^2)
sse_teste_1_Y1 = sum((y1_Y1_teste - Y1_teste)^2)
sst_treino_1 = sum((mean(Y1_treino) - Y1_treino)^2)
sst_teste_1 = sum((mean(Y1_treino) - Y1_teste)^2)
(R2_treino_1_Y1 = 1 - sse_treino_1_Y1/sst_treino_1)[1] 0.9811842(R2_teste_1_Y1 = 1 - sse_teste_1_Y1/sst_teste_1)[1] 0.9777854