Neural Networks

Data preparation

Let’s use the Boston dataset from the MASS package, and suppose we want to predict the median value of houses (last variable) based on the other available variables.

Separate the data into train and test datasets, with 75% of the data in the training dataset, sampled on our dependent variable medv:

data <- MASS::Boston

index <- createDataPartition(data$medv, p = 0.75, list = FALSE, times = 1)
train <- data[index,]
test  <- data[-index,]

Linear regression model for comparison

Calculate the linear regression model on the training set.

lm.fit <- lm(medv ~ ., data=train)
sjPlot::tab_model(lm.fit)

Make predictions on the testing set, and calculate the mean standard error:

pr.lm <- predict(lm.fit, test)
MSE.lm <- sum((pr.lm - test$medv)^2)/nrow(test)
MSE.lm

## [1] 22.89701

Set up a neural network

norm.data <- data %>% 
  # calculate the pre-process parameters from the dataset
  preProcess(method = c("range")) %>% 
  # transform the dataset based on the pre-processing parameters
  predict(data)

# summarize the transformed dataset to check that it worked out well
summary(norm.data) %>% t() %>% knitr::kable()

# rebuild training and testing datasets, but normalized
norm.train <- norm.data[index,]
norm.test  <- norm.data[-index,]

nn <- neuralnet(medv ~ ., data = norm.train, hidden=c(5,3), linear.output=TRUE)
plot(nn, rep = "best")

pr.nn <- predict(nn, newdata = norm.test %>% select(-medv))

pr.nn_ <- pr.nn * (max(data$medv) - min(data$medv)) + min(data$medv)
test.r <- norm.test$medv * (max(data$medv) - min(data$medv)) + min(data$medv)

MSE.nn <- sum((test.r - pr.nn_)^2)/nrow(norm.test)
MSE.nn

## [1] 12.12083

data.frame(
  real = test$medv,
  `Neuronal Net` = pr.nn_,
  `Linear Model` = pr.lm
) %>% 
  gather(key = "variable", value = "Prediction", -real) %>% 
  ggplot() +
    aes(x = real, y = Prediction) +
    geom_abline() +
    geom_point() +
    facet_wrap(~variable) +
    labs(x = "Observed value",
         title = "Performance comparison between methods",
         caption = "Straight line is the 45-degree line") +
    ggthemes::theme_tufte(base_size = 16)

Repeat the analysis many times

# repeat analysis 10 times
MSE <- map_df(1:10, ~{
  
  # create datasets
  index <- createDataPartition(data$medv, p = 0.75, list = FALSE, times = 1)
  train <- data[index, ]
  test  <- data[-index, ]
  
  # the linear regression
  lm.fit <- lm(medv ~ ., data=train)
  pr.lm <- predict(lm.fit, test)
  MSE.lm <- sum((pr.lm - test$medv)^2)/nrow(test)
  
  # normalize datasets
  norm.data <- data %>% 
    preProcess(method = c("range")) %>% 
    predict(data)
  norm.train <- norm.data[index,]
  norm.test  <- norm.data[-index,]
  
  # train the neural net
  nn <- neuralnet(medv ~ ., data = norm.train, hidden=c(5,3), linear.output=TRUE)
  
  # make prediction
  pr.nn <- predict(nn, newdata = norm.test %>% select(-medv))
  
  # de-normalize and compare to observed data
  pr.nn_ <- pr.nn * (max(data$medv) - min(data$medv)) + min(data$medv)
  test.r <- norm.test$medv * (max(data$medv) - min(data$medv)) + min(data$medv)

  MSE.nn <- sum((test.r - pr.nn_)^2)/nrow(norm.test)
  
  return(data.frame(MSE.lm, MSE.nn))
  
  })

MSE %>% knitr::kable()

MSE %>% map_dbl(mean) %>% t() %>% knitr::kable()

The neural net indeed has systematically better predictive power than the linear regression.

Compare different network architectures

The internal structure of the network generally affects its performance. The total number of neurons in the network is not the only thing that matters. Let us explore a few network configurations with 8 neurons:

• 4 hidden layers, 2 neurons each
• 3 hidden layers, 2-4-2
• 2 hidden layers, 4-4
• 1 hidden layer, 8 neurons

# create normalized dataset
norm.data <- data %>% 
  preProcess(method = c("range")) %>% 
  predict(data)

# repeat 100 times on 100 different test/train datasets
MSE <- map(1:100, ~{

  # create training and testing datasets
  index <- createDataPartition(data$medv, p = 0.75, list = FALSE, times = 1)
  norm.train <- norm.data[index,]
  norm.test  <- norm.data[-index,]
  
  # try out different network architectures
  MSE.nn.vector <- 
    list(
      c(2,2,2,2),
      c(2,4,2),
      c(4,4),
      c(8)
      ) %>% 
    map_dbl(function(architecture) {
      
      # train the neural net
      nn <- neuralnet(medv ~ ., data = norm.train, 
                      hidden=architecture, linear.output=TRUE)
  
      # make prediction
      pr.nn <- predict(nn, newdata = norm.test %>% select(-medv))
  
      # de-normalize and compare to observed data
      pr.nn_ <- pr.nn * (max(data$medv) - min(data$medv)) + min(data$medv)
      test.r <- norm.test$medv * (max(data$medv) - min(data$medv)) + min(data$medv)

      MSE.nn <- sum((test.r - pr.nn_)^2)/nrow(norm.test)  
      
      return(MSE.nn)
    })

  return(MSE.nn.vector)
  
  })

# turn list into dataframe
MSE <- do.call(rbind.data.frame, MSE) %>% mutate_all(round,3)
colnames(MSE) <- c("a2222", "a242", "a44", "a8")

MSE %>% DT::datatable()

Calculate the average performance for each architecture:

MSE.means <- data.frame(
    mean  = map_dbl(MSE, mean),
    error = map_dbl(MSE, ~{qt(0.975, df = nrow(MSE)-1) * sd(.x)/sqrt(nrow(MSE))})
  ) %>% rownames_to_column("architecture")

MSE.means %>% 
  ggplot() +
    aes(x = architecture, y = mean, label = round(mean,1)) +
    geom_bar(stat="identity") +
    geom_errorbar(aes(ymin=mean-error, ymax=mean+error), width=.2) +
    geom_label() +
    labs(title = "Performance differences between different network architectures",
         caption = "Error bars show the 95% confidence intervals") +
    ggthemes::theme_tufte(base_size = 16)