Tutorial 4.B: Redes Neurais Artificiais com o pacote Keras (Python)

Este tutorial é baseado nos exemplos do livro An Introduction to Statistical Learning, de Gareth James, Daniela Witten, Trevor Hastie e Rob Tibshirani (https://www.statlearning.com/), em grande parte adaptados do livro Deep Learning with R, de François Chollet e J.J. Allaire.

Antes de começarmos com este tutorial, precisamos considerar alguns pontos relacionados ao pacote keras. Esse pacote é originalmente uma biblioteca de alto nível para construir e treinar modelos de redes neurais em Python.

A versão do keras para o R fornece uma interface para a biblioteca keras do Python, que por sua vez faz interface com outros pacotes do Python, como o TensorFlow, muito popular para construção de redes neurais. Ou seja, a versão do keras para o R é um wrapper: ele atua como uma “ponte”, permitindo que você escreva um script em R para interagir com bibliotecas para redes neurais do Python nos bastidores.

Para que essa “ponte” funcione, você precisa ter uma instalação funcional do Python com a biblioteca keras instalada. O Conda é uma ferramenta popular para gerenciar ambientes e pacotes, incluindo Python e suas bibliotecas (e miniconda é simplemente sua verão mais leve). O pacote keras do R facilita a instalação e o gerenciamento dessas dependências do Python através do Conda.

A vantagem de se dar a esse trabalho todo é que o pacote keras é um dos melhores pacotes disponíveis (e talvez o mais popular) para se trabalhar com redes neurais, e sua implementação é muito menos complexa comparada ao pacote torch que utilizamos no tutorial anterior. (Curiosidade: O pacote torch se tornou disponível como uma alternativa ao pacote keras para deep learning em R. Embora o torch não exija uma instalação do Python, sua implementação atual é consideravelmente mais lenta do que a do keras). Vamos então tentar instalar o keras e suas dependências com o script abaixo. Note que você não precisa rodar a linha reticulate::install_miniconda() se miniconda já estiver instalado em seu computador, e você precisa instalar o pacote keras apenas uma vez (e a instalação pode demorar um pouco):

# install.packages("keras")
# reticulate::install_miniconda()
# keras::install_keras(method = "conda", python_version = "3.10")

Rede de Camada Única

Como no tutorial anterior, também começamos carregando e separando os dados que usaremos (Hitters) em um conjunto de treino e um conjunto de teste:

library(ISLR2)
Gitters <- na.omit(Hitters)
n <- nrow(Gitters)
set.seed(13)
ntest <- trunc(n / 3)
testid <- sample(1:n, ntest)

Primeiro, ajustaremos um modelo linear. Observe o uso do comando with() abaixo: o primeiro argumento é um data frame, e o segundo é uma expressão que pode se referir aos elementos do dataframe pelo nome. Neste caso, o data frame corresponde aos dados de teste e a expressão calcula o erro absoluto médio de predição nesses dados:

lfit <- lm(Salary ~ ., data = Gitters[-testid, ])
lpred <- predict(lfit, Gitters[testid, ])
with(Gitters[testid, ], mean(abs(lpred - Salary)))
[1] 254.6687

Em seguida, normalizamos as variáveis independentes criando um novo objeto chamado x e mantemos os valores originais de nossa variável dependente em um novo objeto chamado y:

x <- scale(model.matrix(Salary ~ . - 1, data = Gitters))
y <- Gitters$Salary

A primeira linha do script acima utiliza a função model.matrix(), que produz a mesma matriz utilizada por lm() (o -1 omite o intercepto). Note que essa função converte automaticamente fatores em variáveis dummies (i.e., com \(0\)s e \(1\)s). A função scale() padroniza a matriz para que cada coluna tenha média \(0\) e variância \(1\):

library(glmnet)
Loading required package: Matrix
Loaded glmnet 4.1-8
cvfit <- cv.glmnet(x[-testid, ], y[-testid],
    type.measure = "mae")
cpred <- predict(cvfit, x[testid, ], s = "lambda.min")
mean(abs(y[testid] - cpred))
[1] 252.2994

Para ajustar a rede neural, primeiro definimos sua estrutura, dessa vez, utilizando o pacote keras (essa parte do script pode demorar um pouco para rodar):

library(keras)
reticulate::use_condaenv(condaenv = "r-tensorflow")
modnn <- keras_model_sequential() %>%
  layer_dense(units = 50, activation = "relu", input_shape = ncol(x)) %>%
  layer_dropout(rate = 0.4) %>%
  layer_dense(units = 1)

Criamos um simples objeto em R (modelo) chamado modnn e adicionamos detalhes sobre suas camadas de maneira sequencial utilizando a função keras_model_sequential().

Note que essa rede é igual a rede do tutorial anterior (onde utilizamos o pacote torch ao inves do keras):

  • O objeto modnn possui \(1\) única camada oculta (hidden layer), definida como layer_dense(units = 50, activation = "relu", input_shape = ncol(x)), com \(50\) unidades (nodes) ou neurônios (neurons) e uma função de ativação (activation) do tipo ReLU ("relu").

  • Em seguida, temos uma “camada” de self$dropout (dropout layer), onde aleatórios \(40\%\) dos \(50\) neurônios da camada anterior são definidos como zero durante cada iteração do algoritmo de descida de gradiente estocástico (stochastic gradient descent, ou SGD).

  • Finalmente, a camada de saída possui apenas \(1\) neurônio (ou seja, convergindo de \(50\) para \(1\): layer_dense(units = 1)), sem função de ativação, indicando que o modelo fornece uma saída quantitativa única.

Em seguida, adicionamos alguns “detalhes” ao modnn que controlam o algoritmo de ajuste, como a função de perda baseada no erro quadrático médio (loss = "mse"). O algoritmo acompanhará erro absoluto médio (MAE, em inglês) nos dados de treinamento e nos dados de validação (metrics = list("mean_absolute_error")):

modnn %>% compile(loss = "mse",
    optimizer = optimizer_rmsprop(),
    metrics = list("mean_absolute_error")
   )

Na linha anterior, o operador pipe passa modnn como o primeiro argumento para compile(). A função compile() não altera o objeto modnn, mas comunica as especificações da nossa rede à instância correspondente em Python.

Agora ajustaremos o modelo com o script abaixo. Fornecemos os dados de treinamento e dois parâmetros de ajuste, epochs e batch_size. Utilizar batch_size = 32 significa que, a cada passo do SGD, o algoritmo seleciona aleatoriamente 32 observações de treinamento para o cálculo do gradiente. Novamente, lembre-se que uma época (epoch) corresponde ao número de passos/iterações necessários para processar/passar \(n\) observações pela rede. Como o conjunto de treinamento tem \(n=176\), uma época é \(176/32=5,5\) passos do SGD.

A função fit() no script abaixo possui um argumento chamado validation_data. Como vimos anteriormente, esses dados não são usados no ajuste, mas podem ser utilizados para acompanhar o progresso do modelo—neste caso, reportando o erro absoluto médio (val_mean_absolute_error). Aqui, fornecemos os dados de teste para que possamos ver o erro tanto dos dados de treinamento quanto dos dados de teste conforme as epoch avançam. (Para ver mais opções de ajuste, digite ?fit.keras.engine.training.Model.)

Note que estamos utilizando epochs = 10 no script abaixo, mas o pacote keras é tão eficiente que poderiamos facilemente aumentar esse número para \(500\) ou mais. (Na verdade, recomendamos que você utilize, pelo menos, epochs = 200 para ver a “magia” da optimização em curso):

history <- modnn %>% fit(
    x[-testid, ], y[-testid], epochs = 10, batch_size = 32,
    validation_data = list(x[testid, ], y[testid])
  )
Epoch 1/10
6/6 - 1s - loss: 455815.9688 - mean_absolute_error: 532.7981 - val_loss: 554588.8125 - val_mean_absolute_error: 538.6646 - 573ms/epoch - 95ms/step
Epoch 2/10
6/6 - 0s - loss: 455384.1875 - mean_absolute_error: 532.5001 - val_loss: 554157.8750 - val_mean_absolute_error: 538.3992 - 44ms/epoch - 7ms/step
Epoch 3/10
6/6 - 0s - loss: 455081.6875 - mean_absolute_error: 532.2415 - val_loss: 553741.8125 - val_mean_absolute_error: 538.1532 - 49ms/epoch - 8ms/step
Epoch 4/10
6/6 - 0s - loss: 454652.5312 - mean_absolute_error: 531.9547 - val_loss: 553328.5625 - val_mean_absolute_error: 537.9091 - 34ms/epoch - 6ms/step
Epoch 5/10
6/6 - 0s - loss: 454349.0000 - mean_absolute_error: 531.7136 - val_loss: 552900.2500 - val_mean_absolute_error: 537.6608 - 49ms/epoch - 8ms/step
Epoch 6/10
6/6 - 0s - loss: 453988.5312 - mean_absolute_error: 531.4333 - val_loss: 552485.4375 - val_mean_absolute_error: 537.4127 - 41ms/epoch - 7ms/step
Epoch 7/10
6/6 - 0s - loss: 453416.4062 - mean_absolute_error: 531.0616 - val_loss: 552017.6875 - val_mean_absolute_error: 537.1510 - 25ms/epoch - 4ms/step
Epoch 8/10
6/6 - 0s - loss: 453107.1250 - mean_absolute_error: 530.8090 - val_loss: 551580.4375 - val_mean_absolute_error: 536.8915 - 56ms/epoch - 9ms/step
Epoch 9/10
6/6 - 0s - loss: 452747.2812 - mean_absolute_error: 530.6142 - val_loss: 551114.6875 - val_mean_absolute_error: 536.6298 - 44ms/epoch - 7ms/step
Epoch 10/10
6/6 - 0s - loss: 452138.6875 - mean_absolute_error: 530.2269 - val_loss: 550633.3125 - val_mean_absolute_error: 536.3476 - 49ms/epoch - 8ms/step

Podemos plotar o history para exibir o erro absoluto médio para os dados de treinamento e de teste. Note que o gráfico não fará muito sentido se você tiver rodado o script acima com um número baixo de epochs… Para uma melhor estética, instale o pacote ggplot2 antes de chamar a função plot(). O ggplot2 é um pacote extremamente popular e flexível para a criação de gráficos no R. Caso não tenha instalado o ggplot2, o script abaixo ainda será executado sem erros, mas o gráfico será menos atraente:

# install.packages("ggplot2")
plot(history)

Finalmente, fazemos a predição com o modelo final e avaliamos seu desempenho nos dados de teste. Novamente, note você tiver rodado o modelo com um número baixo de epochs, provavelmente os resulados serão relativamente ruins… Devido ao uso do SGD, os resultados variam ligeiramente a cada ajuste. Infelizmente, a função set.seed() do R não garante resultados idênticos (já que o ajuste é feito em Python), então seus resultados serão ligeiramente diferentes dos resultados reportados aqui.

npred <- predict(modnn, x[testid, ])
3/3 - 0s - 70ms/epoch - 23ms/step
mean(abs(y[testid] - npred))
[1] 536.3476

Rede Multicamadas

O pacote keras vem com diversos conjuntos de dados de exemplo, incluindo os dados MNIST (Modified National Institute of Standards and Technology database), que utilizamos no tutorial anterior. Nosso primeiro passo é carregar os dados MNIST. A função dataset_mnist() é fornecida para essa finalidade:

mnist <- dataset_mnist()
x_train <- mnist$train$x
g_train <- mnist$train$y
x_test <- mnist$test$x
g_test <- mnist$test$y
dim(x_train)
[1] 60000    28    28
dim(x_test)
[1] 10000    28    28

Como visto anteriormente, há \(60.000\) imagens nos dados de treinamento e \(10.000\) nos dados de teste. As imagens são \(28 \times 28\) e, nesse caso, estão armazenadas como um array tridimensional, então precisamos remodelá-las em uma matriz. Além disso, precisamos codificar os rótulos de classe usando one-hot encoding (i.e., transformação em \(0\)s e \(1\)s). Felizmente, o keras possui muitas funções integradas que fazem isso por nós:

x_train <- array_reshape(x_train, c(nrow(x_train), 784))
x_test <- array_reshape(x_test, c(nrow(x_test), 784))
y_train <- to_categorical(g_train, 10)
y_test <- to_categorical(g_test, 10)

Já mencionamos que redes neurais são um tanto sensíveis à escala das entradas, e que as images originais do MNIST contém valores entre \(0\) e \(255\). Então, fazemos novamente o reescalonamento para o intervalo unitário (entre \(0\) e \(1\)):

x_train <- x_train / 255
x_test <- x_test / 255

Agora estamos prontos para ajustar nossa rede neural (note que o script abaixo é bem mais simples que o script do tutorial anterior baseado no pacote torch):

modelnn <- keras_model_sequential()
modelnn %>%
  layer_dense(units = 256, activation = "relu", input_shape = c(784)) %>%
  layer_dropout(rate = 0.4) %>%
  layer_dense(units = 128, activation = "relu") %>%
  layer_dropout(rate = 0.3) %>%
  layer_dense(units = 10, activation = "softmax")

Vamos interpretar a estrutura da nossa rede:

  • A primeira camada (input_shape) vai de \(28 \times 28 = 784\) neurônios de entrada para uma camada oculta de \(256\) neurônios, e utiliza uma função de ativação ReLU (activation = "relu").

  • Em seguida, os dados são passados através de uma “camada” layer_dropout() para realizar uma regularização de dropout onde \(40\%\) dos \(256\) neurônios restantes são aleatóriamente definidos como zero a cada iteração da rede.

  • A segunda camada oculta vem em seguida, convergindo \(256\) neurônios de entrada em \(128\) neurônios de saída, também utilizando uma função de ativação ReLU (activation = "relu").

  • Em seguida, temos outra “camada” de dropout, dessa vez definindo aleatóriamente \(30\%\) dos \(128\) neurônios restantes como zero a cada iteração da rede.

  • A camada final, ou camada de saída (output layer), converge \(128\) neurônios de entrada em \(10\) neurônios de saída. Note que essa camada utiliza uma função de ativação "softmax", que é a função mais utilizada para o problemas de classificação multiclasse: ela computa probabilidades para cada uma das \(10\) classes consideradas para a classificação, e a classe com a maior probabilidade é a “vencedora”.

Finalmente, usamos summary() para resumir o modelo e para garantir que tudo esteja correto:

summary(modelnn)
Model: "sequential_1"
________________________________________________________________________________
 Layer (type)                       Output Shape                    Param #     
================================================================================
 dense_4 (Dense)                    (None, 256)                     200960      
 dropout_2 (Dropout)                (None, 256)                     0           
 dense_3 (Dense)                    (None, 128)                     32896       
 dropout_1 (Dropout)                (None, 128)                     0           
 dense_2 (Dense)                    (None, 10)                      1290        
================================================================================
Total params: 235146 (918.54 KB)
Trainable params: 235146 (918.54 KB)
Non-trainable params: 0 (0.00 Byte)
________________________________________________________________________________

Curiosidade: Note que os nomes das camadas, como dropout_1 e dense_4, incluem esses underlines “_“. Eles são não tem importancia para nós: são gerados automaticamente de forma sequencial pelo keras e funcionam como identificadores para cada camada. Essa convenção faz parte de como a interface do keras para R interage com o Python.

Em seguida, adicionamos detalhes ao modelo para especificar o algoritmo de ajuste. Ajustamos o modelo minimizando a função de entropia cruzada (loss = "categorical_crossentropy"):

modelnn %>% compile(loss = "categorical_crossentropy",
    optimizer = optimizer_rmsprop(), metrics = c("accuracy")
  )

Agora estamos prontos para começar. A etapa final é fornecer os dados de treinamento e ajustar o modelo (novamente, aumente o número de epochs para melhorar o ajuste):

system.time(
  history <- modelnn %>%
      fit(x_train, y_train, epochs = 10, batch_size = 128,
        validation_split = 0.2)
)
Epoch 1/10
375/375 - 2s - loss: 0.4324 - accuracy: 0.8698 - val_loss: 0.1614 - val_accuracy: 0.9532 - 2s/epoch - 5ms/step
Epoch 2/10
375/375 - 1s - loss: 0.2033 - accuracy: 0.9389 - val_loss: 0.1149 - val_accuracy: 0.9663 - 1s/epoch - 3ms/step
Epoch 3/10
375/375 - 1s - loss: 0.1578 - accuracy: 0.9536 - val_loss: 0.1050 - val_accuracy: 0.9705 - 1s/epoch - 4ms/step
Epoch 4/10
375/375 - 2s - loss: 0.1275 - accuracy: 0.9616 - val_loss: 0.0991 - val_accuracy: 0.9721 - 2s/epoch - 6ms/step
Epoch 5/10
375/375 - 2s - loss: 0.1148 - accuracy: 0.9657 - val_loss: 0.0955 - val_accuracy: 0.9740 - 2s/epoch - 6ms/step
Epoch 6/10
375/375 - 3s - loss: 0.1059 - accuracy: 0.9672 - val_loss: 0.0849 - val_accuracy: 0.9762 - 3s/epoch - 7ms/step
Epoch 7/10
375/375 - 3s - loss: 0.0953 - accuracy: 0.9718 - val_loss: 0.0860 - val_accuracy: 0.9777 - 3s/epoch - 7ms/step
Epoch 8/10
375/375 - 3s - loss: 0.0909 - accuracy: 0.9733 - val_loss: 0.0824 - val_accuracy: 0.9787 - 3s/epoch - 9ms/step
Epoch 9/10
375/375 - 3s - loss: 0.0829 - accuracy: 0.9743 - val_loss: 0.0853 - val_accuracy: 0.9778 - 3s/epoch - 9ms/step
Epoch 10/10
375/375 - 4s - loss: 0.0791 - accuracy: 0.9756 - val_loss: 0.0862 - val_accuracy: 0.9768 - 4s/epoch - 10ms/step
   user  system elapsed 
  12.84    0.80   25.28 
plot(history, smooth = TRUE) #use `smooth = TRUE` para remove a linha de tendência

A saída é um relatório de progresso sobre o ajuste do modelo, agrupado por época. Isso é muito útil, já que em conjuntos de dados grandes o processo de ajuste pode levar tempo.

Aqui especificamos novamente uma divisão de validação de \(20\%\). Ou seja, o treinamento da rede foi realizado em \(80\%\) das \(60.000\) imagens disponíveis nos dados de treinamento (consulte ?fit.keras.engine.training.Model para todos os argumentos de ajuste opcionais). O SGD utilizou lotes de \(128\) observações no cálculo do gradiente. Ou seja, fazendo as contas, vemos que uma epoch corresponde a \(375\) iterações.

Novamente, para obter o erro da rede baseado nos dados de teste, primeiro escrevemos uma função simples accuracy() que compara as classes previstas e reais, e depois a usamos para avaliar nossas previsões:

accuracy <- function(pred, truth)
  mean(drop(as.numeric(pred)) == drop(truth))
modelnn %>% predict(x_test) %>% k_argmax() %>% accuracy(g_test)
313/313 - 1s - 785ms/epoch - 3ms/step
[1] 0.9786

Semelhante ao pacote torch, também podemos usar o keras para ajustar um modelo linear ou logístico através de uma rede neural. Novamente, apenas precisamos de uma camada de entrada e uma camada de saída, e omitimos as camadas ocultas (hidden layers):

modellr <- keras_model_sequential() %>%
  layer_dense(input_shape = 784, units = 10,
       activation = "softmax")
summary(modellr)
Model: "sequential_2"
________________________________________________________________________________
 Layer (type)                       Output Shape                    Param #     
================================================================================
 dense_5 (Dense)                    (None, 10)                      7850        
================================================================================
Total params: 7850 (30.66 KB)
Trainable params: 7850 (30.66 KB)
Non-trainable params: 0 (0.00 Byte)
________________________________________________________________________________

Ajustamos o modelo da mesma forma que antes (aumente o número de epochs para melhorar o ajuste):

modellr %>% compile(loss = "categorical_crossentropy",
     optimizer = optimizer_rmsprop(), metrics = c("accuracy"))
modellr %>% fit(x_train, y_train, epochs = 5,
      batch_size = 128, validation_split = 0.2)
Epoch 1/5
375/375 - 2s - loss: 0.6726 - accuracy: 0.8339 - val_loss: 0.3591 - val_accuracy: 0.9040 - 2s/epoch - 5ms/step
Epoch 2/5
375/375 - 1s - loss: 0.3537 - accuracy: 0.9028 - val_loss: 0.3124 - val_accuracy: 0.9126 - 1s/epoch - 4ms/step
Epoch 3/5
375/375 - 1s - loss: 0.3186 - accuracy: 0.9112 - val_loss: 0.2926 - val_accuracy: 0.9192 - 1s/epoch - 4ms/step
Epoch 4/5
375/375 - 1s - loss: 0.3024 - accuracy: 0.9151 - val_loss: 0.2839 - val_accuracy: 0.9213 - 1s/epoch - 3ms/step
Epoch 5/5
375/375 - 2s - loss: 0.2926 - accuracy: 0.9187 - val_loss: 0.2793 - val_accuracy: 0.9231 - 2s/epoch - 4ms/step
modellr %>% predict(x_test) %>% k_argmax() %>% accuracy(g_test)
313/313 - 1s - 822ms/epoch - 3ms/step
[1] 0.9209

Exercício

Agora é sua vez de ajustar modelos de redes neurais artificiais a um novo conjunto de dados. Você pode fazer isso utilizando seus próprios dados, dados disponíveis na internet ou em pacotes do R. Outra opção é utilizar os dados sobre desmatamento municipal em 2004, disponíveis em https://thaleswest.wixsite.com/home/tutorials (note que o site contém a descrição das variáveis deste conjunto de dados).

Siga os seguintes passos:

  • Primeiro, faça uma inspeção geral dos dados a serem utilizados.
  • Selecione uma variável dependente e ao menos \(15\) variáveis independentes para ajustar sua rede de classificação ou regressão.
  • Realize alguns testes modificando a estrutura da sua rede e seus hiperparâmetros.