Info

Objective

The purpose of writing graded lab reports is to help students to stay on track and to provide summative feedback. Each lab report is just 1% of the total course mark. Please do not cheat - it is not worth it!

Your task

Solve the practical questions, knit your document into a PDF and submit to NTULearn before the deadline. The deadline is very tight because the task is simple. We are sure that everyone is capable to do it by themselves and we want to discourage taking someone else’s report and writing it with your own words.

Deadline

6 Oct 2026, midnight

Libraries

We will work with a dataset of RMS Titanic passengers. The response variable is Survived indicating whether a passenger died or survived the sinking. We will rename Survived to Y and change it to a factor variable with possible values survived and died.

Here, we load libraries, data and set the random seed. Replace the number “1729” with the numeric part of your matric no.

library(tidyverse)   # data manipulation
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.1     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.3     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(fastDummies) # dummy variables
library(patchwork)
library(janitor)
## 
## Attaching package: 'janitor'
## 
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test
library(titanic)     # Titanic data
library(caret)       # confusionMatrix
## Loading required package: lattice
## 
## Attaching package: 'caret'
## 
## The following object is masked from 'package:purrr':
## 
##     lift
library(torch)       # tensors, modules
library(luz)         # high-level training loop 

knitr::opts_chunk$set(
  echo = TRUE,
  warning = FALSE,
  message = FALSE,
  results = 'hold'
)

# replace '1729' with the numeric part of your matric no
torch_manual_seed(1729) 
set.seed(1729)

X <- titanic_train %>%
  select(Survived, Pclass, Sex, Age, SibSp, Parch, Fare, Embarked) %>%
  drop_na()

head(X)

Question 1

Plot the bar chart that shows frequencies (not counts) of survival for men and women. Were men or women more likely to survive the sinking?

X %>% 
  group_by(Sex) %>%
  summarize(fraction_survived = mean(Survived), .groups = "drop") %>%
  ggplot(aes(x = Sex, y = fraction_survived)) +
  geom_col() +
  theme_minimal()

Here is how we can do it with tabyl from janitor package:

X %>% 
  tabyl(Sex, Survived) %>%
  adorn_percentages() %>%
  ggplot(aes(x = Sex, y = `1`)) +
  geom_bar(stat = "identity") +
  theme_bw() + ylab("Fraction survived")

Or post-processing in ggplot:

X %>% 
  ggplot(aes(x = Sex, y = Survived)) +
  stat_summary(fun = mean, geom = "bar") +
  scale_y_continuous(labels = scales::percent) +
  theme_light()

There is a way to include the death rate together with survival rate:

X %>% 
  mutate(Survived = ifelse(Survived, "Yes", "No")) %>%
  ggplot(aes(x = Sex, fill = Survived)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent) +
  labs(
    x = "Sex",
    y = "Percentage",
    fill = "Survival"
  ) +
  theme_linedraw()

In any case, women were more likely to survive the sinking.

Preparing the data (for torch)

We will:

  1. Split the response from predictors;
  2. One-hot encode categorical variables;
  3. Convert to matrices, then tensors;
  4. Split into training and test sets;
  5. Wrap tensors in dataloaders.
# Split indices (70/30)
ind <- which(runif(nrow(X)) < 0.7)

# One-hot encode predictors
x_all <- X %>%
  select(-Survived) %>%
  dummy_cols(remove_first_dummy = TRUE, remove_selected_columns = TRUE) %>%
  as.matrix()

y_all <- X$Survived %>% as.numeric()  # 0/1

# Train/test matrices
x_train <- x_all[ind, , drop = FALSE]
y_train <- y_all[ind]

x_test  <- x_all[-ind, , drop = FALSE]
y_test  <- y_all[-ind]

cat("Training data dimensions =", dim(x_train), "\n")
cat("Test data dimensions =", dim(x_test), "\n")

# 1) targets must be float and [N,1]
y_train_t <- torch_tensor(y_train, dtype = torch_float())$unsqueeze(2)  # [N,1]
y_test_t  <- torch_tensor(y_test,  dtype = torch_float())$unsqueeze(2)  # [N,1]

x_train_t <- torch_tensor(x_train, dtype = torch_float())
x_test_t  <- torch_tensor(x_test,  dtype = torch_float())

train_ds <- tensor_dataset(x_train_t, y_train_t)
test_ds  <- tensor_dataset(x_test_t,  y_test_t)

# 2) macOS stability: avoid extra workers and limit threads
train_dl <- dataloader(train_ds, batch_size = 64, shuffle = TRUE, num_workers = 0)
test_dl  <- dataloader(test_ds,  batch_size = 64, num_workers = 0)

torch_set_num_threads(2)  # optional but helps prevent UI freezes
# device <- torch_device("cpu")  # optional: force CPU if MPS causes issues
## Training data dimensions = 522 9 
## Test data dimensions = 192 9

Here is a helper function for you:

extract_y <- function(torch_data) {
  torch_data[1:length(torch_data)][[2]] %>% as.numeric()
}

extract_y(test_ds)
##   [1] 0 1 1 0 0 1 0 0 0 1 0 1 0 1 0 1 1 1 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 0 1 0
##  [38] 0 1 1 0 0 1 0 1 1 0 0 0 0 1 0 1 0 1 1 0 0 0 0 0 1 1 1 1 0 0 1 0 1 1 0 0 0
##  [75] 0 0 1 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 1 0 1 1 1 1 1 0 0 1 0 1 1 1 0 0 1 1 1
## [112] 1 0 1 1 0 0 1 1 0 0 1 1 1 1 1 0 1 1 1 0 1 1 0 0 1 1 0 1 0 0 1 1 0 0 0 0 0
## [149] 0 0 0 0 1 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 1 1 0 1 1 0 0 0
## [186] 0 0 1 1 1 1 1

Question 2

Here we define and train an artificial neural network that predicts the survival of a Titanic passenger. This is done incorrectly. Figure out what is wrong, correct the error, and train the model.

bad_net <- nn_module(
  initialize = function(input_dim) {
    self$fc1 <- nn_linear(input_dim, 1000)
    self$fc2 <- nn_linear(1000, 1000)
    self$fc3 <- nn_linear(1000, 1000)
    self$fc4 <- nn_linear(1000, 1000)
    self$fc5 <- nn_linear(1000, 1000)
    self$out <- nn_linear(1000, 1)        
  },
  forward = function(x) {
    x %>% self$fc1() %>% nnf_relu() %>%
      self$fc2() %>% nnf_relu() %>%
      self$fc3() %>% nnf_relu() %>%
      self$fc4() %>% nnf_relu() %>%
      self$fc5() %>% nnf_relu() %>%
      self$out()
  }
)

bad_fit <- bad_net %>%                            
  setup(
    loss = nn_l1_loss(),                          
    metrics = list(luz_metric_binary_accuracy()),
    optimizer = optim_rmsprop
  ) %>%
  set_hparams(input_dim = ncol(x_train)) %>%      
  set_opt_hparams(lr = 1e-3) %>%
  fit(train_dl, epochs = 20, valid_data = test_dl, verbose = TRUE)

bad_fit
## A `luz_module_fitted`
## ── Time ────────────────────────────────────────────────────────────────────────
## • Total time: 2.2s
## • Avg time per training epoch: 100ms
## 
## ── Results ─────────────────────────────────────────────────────────────────────
## Metrics observed in the last epoch.
## 
## ℹ Training:
## loss: 0.3399
## acc: 0.6667
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## An `nn_module` containing 4,015,001 parameters.
## 
## ── Modules ─────────────────────────────────────────────────────────────────────
## • fc1: <nn_linear> #10,000 parameters
## • fc2: <nn_linear> #1,001,000 parameters
## • fc3: <nn_linear> #1,001,000 parameters
## • fc4: <nn_linear> #1,001,000 parameters
## • fc5: <nn_linear> #1,001,000 parameters
## • out: <nn_linear> #1,001 parameters

Answer: For binary classification we should use the sigmoid activation in the output layer. There are two ways to do it:

  • Option A: use nnf_sigmoid activation in the last layer and nn_binary_cross_entropy() loss, or
  • Option B: (preferred) linear activation in the last layer and nn_bce_with_logits_loss(), which applies a numerically stable sigmoid internally.

We will implement Option B and we will train two versions of the model, with a high and a small learning rate.

good_net <- nn_module(
  initialize = function(input_dim) {
    self$fc1 <- nn_linear(input_dim, 1000)
    self$fc2 <- nn_linear(1000, 1000)
    self$fc3 <- nn_linear(1000, 1000)
    self$fc4 <- nn_linear(1000, 1000)
    self$fc5 <- nn_linear(1000, 1000)
    self$out <- nn_linear(1000, 1)        
  },
  forward = function(x) {
    x %>% self$fc1() %>% nnf_relu() %>%
      self$fc2() %>% nnf_relu() %>%
      self$fc3() %>% nnf_relu() %>%
      self$fc4() %>% nnf_relu() %>%
      self$fc5() %>% nnf_relu() %>%
      self$out()
  }
)

good_fit_1 <- good_net %>%                            
  setup(
    loss = nn_bce_with_logits_loss(),                          
    metrics = list(luz_metric_binary_accuracy()),
    optimizer = optim_rmsprop
  ) %>%
  set_hparams(input_dim = ncol(x_train)) %>%      
  set_opt_hparams(lr = 1e-3) %>%
  fit(train_dl, epochs = 100, valid_data = test_dl, verbose = FALSE)

good_fit_2 <- good_net %>%                            
  setup(
    loss = nn_bce_with_logits_loss(),                          
    metrics = list(luz_metric_binary_accuracy()),
    optimizer = optim_rmsprop
  ) %>%
  set_hparams(input_dim = ncol(x_train)) %>%      
  set_opt_hparams(lr = 1e-6) %>%
  fit(train_dl, epochs = 100, valid_data = test_dl, verbose = FALSE)

cat("Model with high learning rate:\n")
good_fit_1
cat("Model with low learning rate:\n")
good_fit_2
## Model with high learning rate:
## A `luz_module_fitted`
## ── Time ────────────────────────────────────────────────────────────────────────
## • Total time: 13.8s
## • Avg time per training epoch: 126ms
## 
## ── Results ─────────────────────────────────────────────────────────────────────
## Metrics observed in the last epoch.
## 
## ℹ Training:
## loss: 0.5897
## acc: 0.6762
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## An `nn_module` containing 4,015,001 parameters.
## 
## ── Modules ─────────────────────────────────────────────────────────────────────
## • fc1: <nn_linear> #10,000 parameters
## • fc2: <nn_linear> #1,001,000 parameters
## • fc3: <nn_linear> #1,001,000 parameters
## • fc4: <nn_linear> #1,001,000 parameters
## • fc5: <nn_linear> #1,001,000 parameters
## • out: <nn_linear> #1,001 parameters
## Model with low learning rate:
## A `luz_module_fitted`
## ── Time ────────────────────────────────────────────────────────────────────────
## • Total time: 12.8s
## • Avg time per training epoch: 116ms
## 
## ── Results ─────────────────────────────────────────────────────────────────────
## Metrics observed in the last epoch.
## 
## ℹ Training:
## loss: 0.6152
## acc: 0.659
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## An `nn_module` containing 4,015,001 parameters.
## 
## ── Modules ─────────────────────────────────────────────────────────────────────
## • fc1: <nn_linear> #10,000 parameters
## • fc2: <nn_linear> #1,001,000 parameters
## • fc3: <nn_linear> #1,001,000 parameters
## • fc4: <nn_linear> #1,001,000 parameters
## • fc5: <nn_linear> #1,001,000 parameters
## • out: <nn_linear> #1,001 parameters

Remark: In practice, such a small learning rate is unusual; we use it here just to provoke a visible effect.

Question 3

What is the issue with the model you trained in Question 2? Choose the correct option:

  • The learning rate is too high (gradient descent overshoots)
  • The learning rate is too low (convergence is very slow)
  • The model overfits
  • The model underfits

Or is the model just right? Explain your answer.

Answer: The default learning rate is too high - after 2 or 3 epochs, both training and test accuracy will fluctuate between 0.6 and 0.7, which is a sign that gradient descent overshoots. This is a sign that we should reduce it. However, setting the learning rate to \(10^{-6}\) will cause slow convergence, which we can see in these diagnostic plots:

history_of_training <- function(mlp_model) {
  # extract per-epoch losses
  train_loss <- map_dbl(mlp_model$records$metrics$train, ~ .x$loss)
  valid_loss <- map_dbl(mlp_model$records$metrics$valid, ~ .x$loss)
  
  # extract per-epoch accuracy
  train_acc <- map_dbl(mlp_model$records$metrics$train, ~ .x$acc)
  valid_acc <- map_dbl(mlp_model$records$metrics$valid, ~ .x$acc)
  
  tibble(
    epoch = seq_along(train_loss),
    train_loss = train_loss,
    valid_loss = valid_loss,
    train_acc = train_acc,
    valid_acc = valid_acc
  )
}

log_df_1 <- history_of_training(good_fit_1)
log_df_2 <- history_of_training(good_fit_2)

plot_acc_1 <- log_df_1 %>%
  pivot_longer(cols = ends_with("acc"), names_to = "Metric") %>%
  ggplot(aes(x = epoch, y = value, color = Metric)) +
  geom_line() + geom_point() + theme_minimal() + ggtitle("High Learning Rate")

plot_acc_2 <- log_df_2 %>%
  pivot_longer(cols = ends_with("acc"), names_to = "Metric") %>%
  ggplot(aes(x = epoch, y = value, color = Metric)) +
  geom_line() + geom_point() + theme_minimal() + ggtitle("Low Learning Rate")


plot_acc_1 / plot_acc_2

Here are the training and the test accuracy:

binary_prediction <- function(torch_model,
                              torch_data = test_ds) {
  pred_probs <- torch_model %>%
    predict(torch_data) %>%
    nnf_sigmoid() %>%
    as.numeric()
  
  0 + (pred_probs > 0.5)
}

torch_binary_acc <- function(torch_model, 
                      torch_data = test_ds) {
  preds <- binary_prediction(torch_model, torch_data)
  actual_y <- extract_y(torch_data)
  mean(preds == actual_y)
}


torch_conf_mat <- function(torch_model, 
                           torch_data = test_ds, 
                           predefined_classes = torch_data$classes) {
  preds <- binary_prediction(torch_model, torch_data) %>%
    ifelse("Died", "Survived") %>%
    as_factor()
  actual_y <- extract_y(torch_data) %>%
    ifelse("Died", "Survived") %>%
    as_factor()
  confusionMatrix(actual_y, preds)
}

cat("Training accuracy with high learning rate:",
    torch_binary_acc(good_fit_1, train_ds), "\n")

cat("Test accuracy with high learning rate:",
    torch_binary_acc(good_fit_1, test_ds), "\n")

cat("\nTraining accuracy with low learning rate:",
    torch_binary_acc(good_fit_2, train_ds),"\n")

cat("Test accuracy with low learning rate:",
    torch_binary_acc(good_fit_2, test_ds),"\n")
## Training accuracy with high learning rate: 0.6302682 
## Test accuracy with high learning rate: 0.6666667 
## 
## Training accuracy with low learning rate: 0.6992337 
## Test accuracy with low learning rate: 0.6666667

Depending on the initial random seed, we might observe what looks like overfitting - training accuracy might be much higher than the test accuracy. However, the true reason is lack of normalization. Note that independent variables have very different scale:

X %>%
  slice(ind) %>%
  select(where(is.numeric)) %>%
  summarise(across(everything(), 
                   list(min = ~ min(.x), max = ~ max(.x), sd  = ~ sd(.x)))) %>%
  pivot_longer(
    everything(),
    names_to = c("variable", ".value"),
    names_sep = "_"
  )

It causes gradients to be very imbalanced and the gradient descent hard to converge.

Question 4

What is your strategy to fix the issue you identified in Question 3? Explain your strategy and implement it by training a new model. Does your approach solve the issue?

Remark: if you know the correct strategy, but do not have time to tune your new model, just write the correct code in R and explain what you would have done if you had all the time in the world.

Answer: We will first standardize the data:

### First, we compute means and standard deviation of all the columns in x_train:
X_mean <- x_train %>% as_tibble() %>% map_dbl(mean)
X_sd   <- x_train %>% as_tibble() %>% map_dbl(sd)

### Then apply standardization: 
standardize_data <- function(x_matrix, x_mean = X_mean, x_sd = X_sd) {
  x_matrix %>% as_tibble() %>%
    mutate(across(everything(), 
                ~ (.x - x_mean[cur_column()]) / x_sd[cur_column()])) %>%
    as.matrix()
}

x_train_norm <- x_train %>% standardize_data()
x_test_norm <- x_test %>% standardize_data()

### Since we don't need the torch tensors with raw data anymore, we will just overwrite them
x_train_t <- torch_tensor(x_train_norm, dtype = torch_float())
x_test_t  <- torch_tensor(x_test_norm,  dtype = torch_float())

train_ds_norm <- tensor_dataset(x_train_t, y_train_t)
test_ds_norm  <- tensor_dataset(x_test_t,  y_test_t)

# 2) macOS stability: avoid extra workers and limit threads
train_dl <- dataloader(train_ds, batch_size = 64, shuffle = TRUE, num_workers = 0)
test_dl  <- dataloader(test_ds,  batch_size = 64, num_workers = 0)

torch_set_num_threads(2)  # optional but helps prevent UI freezes
# device <- torch_device("cpu")  # optional: force CPU if MPS causes issues

Then, since the model is huge, to prevent overfitting, we will reduce the number of parameters and introduce regularization. We will also change the optimizer - apparently, Adam works better than RMSProp for this problem.

better_net <- nn_module(
  initialize = function(input_dim, p_drop = 0.2) {
    self$fc1 <- nn_linear(input_dim, 200)
    self$fc2 <- nn_linear(200, 100)
    self$fc3 <- nn_linear(100, 50)
    self$drop <- nn_dropout(p = p_drop)
    self$out <- nn_linear(50, 1)
  },
  forward = function(x) {
    x %>% self$fc1() %>% nnf_relu() %>% self$drop() %>%
      self$fc2() %>% nnf_relu() %>% self$drop() %>%
      self$fc3() %>% nnf_relu() %>%
      self$out()
  }
)


better_fit <- better_net %>%
  setup(
    loss = nn_bce_with_logits_loss(),
    metrics = list(luz_metric_binary_accuracy()),
    optimizer = optim_adamw
  ) %>%
  set_hparams(input_dim = ncol(x_train)) %>%      
  set_opt_hparams(weight_decay = 1e-2) %>%
  fit(train_dl, epochs = 150, valid_data = test_dl, verbose = FALSE)

cat("Training accuracy (new model):",
    torch_binary_acc(better_fit, train_ds))

cat("\nTest accuracy (new model):",
    torch_binary_acc(better_fit, test_ds),"\n")
## Training accuracy (new model): 0.8180077
## Test accuracy (new model): 0.8072917

Here is the history of training (the learning rate is too high, but it does the job):

better_fit %>%
  history_of_training() %>%
  pivot_longer(cols = ends_with("acc"), names_to = "Metric") %>%
  ggplot(aes(x = epoch, y = value, color = Metric)) +
  geom_line() + geom_point() + theme_minimal() + ggtitle("Best Model")

REMARK Note that finding a good combination of layer dimensions, optimizer, learning rate, weight decay etc is a challenging task. In practice, it should be done by a systematic search.

For instance, a simple logistic regression is not going to be much less accurate and sometimes even more accurate than a very complex model:

### We will drop Embarked because it has 2 missing values, it doesn't make the model much worse

tit_train <- X %>% select(-Embarked) %>%
  mutate(Survived = ifelse(Survived, "yes", "no")) %>% slice(ind)

tit_test <- X %>% select(-Embarked) %>%
  mutate(Survived = ifelse(Survived, "yes", "no")) %>% slice(-ind)

logistic_model <- train(Survived ~ ., tit_train,
                        method = "glm", family = "binomial")
y_pred <- logistic_model %>% predict(tit_test) 
cat("Accuracy of logistic regression is", 
    round(100 * mean(y_pred == tit_test$Survived), 1), "%\n")
## Accuracy of logistic regression is 81.8 %

Even if we only use Sex as predictor, the logistic model (which just predicts “yes” for women and “no” for men) is better than the first attempt at a neural network with 4M parameters!!!

### We will drop Embarked because it has 2 missing values, it doesn't make the model much worse
sex_model <- train(Survived ~ Sex, tit_train,
                        method = "glm", family = "binomial")
y_pred <- sex_model %>% predict(tit_test) 
cat("Accuracy of logistic regression based on Sex is", 
    round(100 * mean(y_pred == tit_test$Survived), 1), "%\n")
## Accuracy of logistic regression based on Sex is 80.2 %

Declaration of Generative AI usage

I used ChatGPT 5.0 to translate the old handout based on keras library to this new version based on torch library.

Type your name to confirm: Fedor Duzhin