Capitolo 22 Reti Neurali (Neural Networks: NN)

22.1 Introduzione

La letteratura statistica sulle tecniche di previsione/modellazione sta crescendo molto rapidamente. Ogni giorno compaiono nuove tecniche per progetti di data-mining e business intelligence, e gli sviluppi più recenti permettono al ricercatore di ottenere previsioni sempre più precise.
In questo capitolo vorrei introdurre molto brevemente le basi dell’analisi delle reti neurali per problemi di regressione.

22.2 Esempio: Boston housing

Questo esempio userà i dati di Boston housing già visti negli esempi precedenti.
The aim of study is to find a prediction model to assess the probability of die for each passenger based on its Age, Gender, and Class of accomodation.

Rivediamo un riepilogo del dataset:

summary(bostonhousing)

##       crim                zn             indus       chas         nox        
##  Min.   : 0.00632   Min.   :  0.00   Min.   : 0.46   0:471   Min.   :0.3850  
##  1st Qu.: 0.08204   1st Qu.:  0.00   1st Qu.: 5.19   1: 35   1st Qu.:0.4490  
##  Median : 0.25651   Median :  0.00   Median : 9.69           Median :0.5380  
##  Mean   : 3.61352   Mean   : 11.36   Mean   :11.14           Mean   :0.5547  
##  3rd Qu.: 3.67708   3rd Qu.: 12.50   3rd Qu.:18.10           3rd Qu.:0.6240  
##  Max.   :88.97620   Max.   :100.00   Max.   :27.74           Max.   :0.8710  
##        rm             age              dis              rad        
##  Min.   :3.561   Min.   :  2.90   Min.   : 1.130   Min.   : 1.000  
##  1st Qu.:5.886   1st Qu.: 45.02   1st Qu.: 2.100   1st Qu.: 4.000  
##  Median :6.208   Median : 77.50   Median : 3.207   Median : 5.000  
##  Mean   :6.285   Mean   : 68.57   Mean   : 3.795   Mean   : 9.549  
##  3rd Qu.:6.623   3rd Qu.: 94.08   3rd Qu.: 5.188   3rd Qu.:24.000  
##  Max.   :8.780   Max.   :100.00   Max.   :12.127   Max.   :24.000  
##       tax           ptratio            b              lstat            medv      
##  Min.   :187.0   Min.   :12.60   Min.   :  0.32   Min.   : 1.73   Min.   : 5.00  
##  1st Qu.:279.0   1st Qu.:17.40   1st Qu.:375.38   1st Qu.: 6.95   1st Qu.:17.02  
##  Median :330.0   Median :19.05   Median :391.44   Median :11.36   Median :21.20  
##  Mean   :408.2   Mean   :18.46   Mean   :356.67   Mean   :12.65   Mean   :22.53  
##  3rd Qu.:666.0   3rd Qu.:20.20   3rd Qu.:396.23   3rd Qu.:16.95   3rd Qu.:25.00  
##  Max.   :711.0   Max.   :22.00   Max.   :396.90   Max.   :37.97   Max.   :50.00

E quindi carichiamo i package per l’analisi

library(nnet)
library(dplyr)

Il grafico che segue mostra le relazioni tra tutte le colonne del data frame:

library(ggplot2)
library(GGally)
ggpairs(bostonhousing)

Figura 17.1: Relazioni tra le variabili

Appaiono alcune relazioni; soprattutto la variabile dipendente medv sembra relazionarsi con diverse delle altre variabili (indipendenti).

Prepariamo quindi dapprima i dati per l’analisi

set.seed(1000)
bostonhousing$lstat <- log(bostonhousing$lstat)
bostonhousing$rm <- bostonhousing$rm^2
bostonhousing$chas <- factor(bostonhousing$chas, levels = 0:1, labels = c("no", "yes"))
bostonhousing$rad <- factor(bostonhousing$rad, ordered = TRUE)

bostonhousing <- bostonhousing %>%
  dplyr::select(medv, age, lstat, rm, zn, indus, chas, nox, age,
         dis, rad, tax, crim, b, ptratio)

train <- sample(nrow(bostonhousing), 400)
bh_train <- bostonhousing[train,]
bh_test <-  bostonhousing[-train,]

E ora applichiamo la rete neurale. Il package R nnet fornisce una funzione nnet() che permette di adattare una rete neurale con un solo strato nascosto ai dati. nnet() produce in output un oggetto di classe nnet.formula e nnet.

La sintassi per costruire un oggetto nnet è molto simile a quella usata per i modelli lineari (generalizzati). Si noti come per riprodurre l’analisi con le reti neurali sia necessario impostare il seme per il generatore dei numeri casuali:

set.seed(100)

nn0 <- nnet(medv ~ ., data = bh_train,
            size=1, linout=TRUE, skip=TRUE, MaxNWts=10000, trace=FALSE, maxit=500)

Il parametro size della chiamata nnet() qui sopra specifica il numero di unità (neuroni) dello strato nascosto.

Ovviamente possiamo produrre un grafico dei valori osservati Vs. previsti, per l’insieme di test e per quello di training:

data_gr <- bh_train %>%
  mutate(set="train") %>%
  bind_rows(bh_test %>% mutate(set="test"))

data_gr$fit <- predict(nn0, data_gr)

mse <- data_gr %>%
  filter(set=="test") %>%
  summarise(mse = mean((fit-medv)^2)) %>%
  pull()

print(mse)

## [1] 22.60025

ggp <- ggplot(data = data_gr, mapping = aes(x=fit, y=medv)) +
  geom_point(aes(colour=set), alpha=0.6) +
  geom_abline(slope=1, intercept = 0) +
  geom_smooth(method = "lm", se = FALSE, aes(colour=set), alpha=0.6)
print(ggp)

Con un errore quadratico medio per l’insieme di test pari a 22.6002468; il doppio delle foreste casuali.
Possiamo quindi provare ad aumentare in maniera importante il numero di neuroni dello strato nascorso ed il numero massimo di iterazioni per la stima dei parametri; vediamo se questo permetterà di ottenere risultati migliori:

set.seed(100)

nn0 <- nnet(medv ~ ., data = bh_train,
            size=5, linout=TRUE, skip=TRUE, MaxNWts=10000, trace=FALSE, maxit=1000)

Produciamo un grafico dei valori osservati Vs. previsti, per l’insieme di test e per quello di training:

data_gr <- bh_train %>%
  mutate(set="train") %>%
  bind_rows(bh_test %>% mutate(set="test"))

data_gr$fit <- predict(nn0, data_gr)

mse <- data_gr %>%
  filter(set=="test") %>%
  summarise(mse = mean((fit-medv)^2)) %>%
  pull()

print(mse)

## [1] 21.04216

ggp <- ggplot(data = data_gr, mapping = aes(x=fit, y=medv)) +
  geom_point(aes(colour=set), alpha=0.6) +
  geom_abline(slope=1, intercept = 0) +
  geom_smooth(method = "lm", se = FALSE, aes(colour=set), alpha=0.6)
print(ggp)

L’errore quadratico medio per l’insieme di test è pari a 21.0421595, migliore del precedente.
Ovviamente si potrebbe approfondire la ricerca dei parametri ottimi per il modello a reti neurali, ma non lo faremo qui

22.3 Esempio: Boston housing deep learning

In questo esempio applicheremo una rete neurale “profonda” (“deep”) ai dati di Boston housing per prevedere il valore degli stabili a partire dalle variabili indipendenti.
Useremo la libreria keras.
I dati dell’analisi sono sempre relativi al valore delle case a Boston, come nell’esempio precedente, tuttavia le variabili indipendenti, nonché gli insiemi di training e test, sono diversi.
I risultati di questo modello, quindi, non sono confrontabili con quelli del modello di rete neurale precednte.

Cominciamo quindi col richiamare la libreria keras e caricare/preparare i dati dell’analisi:

# https://tensorflow.rstudio.com/keras/articles/tutorial_basic_regression.html
library(keras)
set.seed(100)

boston_housing <- dataset_boston_housing()

c(train_data, train_labels) %<-% boston_housing$train
c(test_data, test_labels) %<-% boston_housing$test

dim(train_data)

## [1] 404  13

dim(test_data)

## [1] 102  13

# Aggiungiamo i nomi di colonna per maggiore leggibilità
library(tibble)

train_df <- as_tibble(train_data)
colnames(train_df) <- c('CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 
                  'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT')
# Mostra le prime righe delle "feature" (variabili indipendenti)
train_df

## # A tibble: 404 x 13
##      CRIM    ZN INDUS  CHAS   NOX    RM   AGE   DIS   RAD   TAX PTRATIO     B
##     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl> <dbl>
##  1 1.23     0    8.14     0 0.538  6.14  91.7  3.98     4   307    21    397.
##  2 0.0218  82.5  2.03     0 0.415  7.61  15.7  6.27     2   348    14.7  395.
##  3 4.90     0   18.1      0 0.631  4.97 100    1.33    24   666    20.2  376.
##  4 0.0396   0    5.19     0 0.515  6.04  34.5  5.99     5   224    20.2  397.
##  5 3.69     0   18.1      0 0.713  6.38  88.4  2.57    24   666    20.2  391.
##  6 0.284    0    7.38     0 0.493  5.71  74.3  4.72     5   287    19.6  391.
##  7 9.19     0   18.1      0 0.7    5.54 100    1.58    24   666    20.2  397.
##  8 4.10     0   19.6      0 0.871  5.47 100    1.41     5   403    14.7  397.
##  9 2.16     0   19.6      0 0.871  5.63 100    1.52     5   403    14.7  169.
## 10 1.63     0   21.9      0 0.624  5.02 100    1.44     4   437    21.2  397.
## # … with 394 more rows, and 1 more variable: LSTAT <dbl>

# Mostra i primi 10 valori della risposta
train_labels[1:10] 

##  [1] 15.2 42.3 50.0 21.1 17.7 18.5 11.3 15.6 15.6 14.4

# Normalizza dati di training delle variabili indipendenti
train_data <- scale(train_data) 

# Ricava medie e deviazioni standard dal set di training per normalizzare il set di test
col_means_train <- attr(train_data, "scaled:center") 
col_stddevs_train <- attr(train_data, "scaled:scale")
test_data <- scale(test_data, center = col_means_train, scale = col_stddevs_train)

A questo punto creiamo una funzione che genera e compila la rete.
La rete è fatta dai neuroni di input, due strati intermedi entrambi di 64 neuroni, e lo strato di output, con un solo neurone.

# Crea una funzione per impostare la rete neurale con tutti i suoi parametri
build_model <- function() {
  
  model <- keras_model_sequential() %>%
    layer_dense(units = 64, activation = "relu",
                input_shape = dim(train_data)[2]) %>%
    layer_dense(units = 64, activation = "relu") %>%
    layer_dense(units = 1)
  
  model %>% compile(
    loss = "mse",
    optimizer = optimizer_rmsprop(),
    metrics = list("mean_squared_error")
  )
  
  return(model)
}

# Costruisce il modello
model <- build_model()
# Riepilogo del modello
summary(model)

## Model: "sequential"
## __________________________________________________________________________________
## Layer (type)                        Output Shape                     Param #      
## ==================================================================================
## dense (Dense)                       (None, 64)                       896          
## __________________________________________________________________________________
## dense_1 (Dense)                     (None, 64)                       4160         
## __________________________________________________________________________________
## dense_2 (Dense)                     (None, 1)                        65           
## ==================================================================================
## Total params: 5,121
## Trainable params: 5,121
## Non-trainable params: 0
## __________________________________________________________________________________

Cerchiamo i parametri che meglio adattano la rete neurale ai dati di training. Quindi vediamo quali sono i valori della funzione di perdita per la rete stimata tramite valutazione sui dati di validazione.
Alla fine valutiamo il valore di mse per i dati di test e facciamo il grafico dei valori osservati Vs. valori previsti per insieme di training e insieme di test.

history <- model %>% fit(
  train_data,
  train_labels,
  epochs = 200,
  validation_split = 0.2,
  verbose = 0
)

library(ggplot2)

plot(history, metrics = "mean_squared_error", smooth = FALSE) +
  coord_cartesian(ylim = c(0, 20))

# Il parametro patience è l'ammontare di epoche per verificare miglioramenti nel
# campione di validazione.
early_stop <- callback_early_stopping(monitor = "val_loss", patience = 20)

model <- build_model()
history <- model %>% fit(
  train_data,
  train_labels,
  epochs = 40,
  validation_split = 0.2,
  verbose = 0,
  callbacks = list(early_stop)
)

# Il numero di epoche (cicli di ottimizzazione) per ottenere il migliore
# errore di validazione è circa pari a 40:
plot(history, metrics = "mean_squared_error", smooth = FALSE) +
  coord_cartesian(xlim = c(0, 50), ylim = c(0, 20))

c(loss, mse) %<-% (model %>% evaluate(test_data, test_labels, verbose = 0))
print(paste("mse:", mse))

## [1] "mse: 24.2888965606689"

test_predictions <- model %>% predict(test_data)
train_predictions <- model %>% predict(train_data)

data_gr <- data.frame(fit=c(test_predictions, train_predictions),
                      medv=c(test_labels, train_labels),
                      set=c(rep("test", nrow(test_predictions)), rep("train", nrow(train_predictions))))
ggp <- ggplot(data = data_gr, mapping = aes(x=fit, y=medv)) +
  geom_point(aes(colour=set), alpha=0.6) +
  geom_abline(slope=1, intercept = 0) +
  geom_smooth(method = "lm", se = FALSE, aes(colour=set), alpha=0.6)
print(ggp)

Il valore di mse per l’insieme di test è pari a 24.29.

22.4 Alcuni cenni di teoria sulle Reti Neurali

Le reti neurali possono essere considerate un tipo di regressione non linerare che prende un insieme di input (variabili esplicative), le trasforma e pesa con un insieme di unità nascoste e strati nascosti per produrre poi un insieme di output o previsioni (a loro volta trasformate).

La fgura che segue è un esempio di una rete neurale feed forward composta di quattro input, uno strato intermedio nascosto che contiene tre unità, e uno strato di output che contiene due output.

Gli output dei nodi in uno strato sono gli input dello strato successivo. Gli input di ciascun nodo sono combinati usando una combinazione lineare pesata. Il risultato è quindi usualmente modificata tramite una funnzione non lineare prima di essere passata in output. Per esempio, gli input nel neurone nascosto \(j\) nella figura precedente sono combinati per dare

\(z_j=b_j+\sum_{i=1}^4 w_{i,j} x_i\).

Nello strato nascosto, questo è quindi modificato usando una funzione “di attivazione” non lineare quale la sigmoide (come detto prima, in questo esempio questa operazione non è eseguita),

\(\phi(z)=\dfrac{1}{1+e^{-z}}\),

per produrre l’input dello strato successivo. Questo permette alla formula di ridurre l’effetto dei valori di input estremi, e quindi rendere la rete più robusta alla presenza di outlier. Inoltre, le funzioni di attivazione aggiungono una componente non lineare alla rete, che altrimenti risulterebbe in una semplice funzione lineare, per quanto complicata.

I valori dei parametri \(b_1\),\(b_2\),\(b_3\) e \(w_{1,1}, \cdots ,w_{4,3}\) sono “appresi” a partire dai dati.

I pesi usualmente assumono inizialmente valori casuali, e sono quindi aggiornati usando i dati osservati. Di conseguenza, si ha un elemento di casualità nelle previsioni prodotte da una rete neurale; di conseguenza la rete è usualmente addestrata diverse volte usando punti di partenza differenti, producendo quindi una media dei diversi risultati.

Il numero di strati nascosti, nonché il numero di nodi in ogni strato nascosto, deve essere specificato prima dell’analisi.