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

13 Oct 2025, midnight

Libraries

We will work with a dataset of information on HDB blocks that comes together with this lab. The original dataset is available here:

https://data.gov.sg/dataset/hdb-property-information

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

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE,
                      out.width = "90%", out.height = "40%",
                      fig.align='center')

library(tidyverse)    # for manipulation with data

## ── 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(GGally)       # Plot of the whole dataset
library(ggfortify)    # Principal components

## Warning: package 'ggfortify' was built under R version 4.6.1

library(ggdendro)     # Hierarchical clustering
library(janitor)      # Helper functions

## 
## Attaching package: 'janitor'
## 
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test

library(broom)        # For tidying

X <- read_csv("hdb_info.csv") %>%
  select(blk_no, street, max_floor_lvl, year_completed, 
         residential, commercial, market_hawker, miscellaneous, multistorey_carpark,
         bldg_contract_town, total_dwelling_units, 
         `1room_sold`, `2room_sold`, `3room_sold`, `4room_sold`, 
         `exec_sold`, multigen_sold, studio_apartment_sold)

## Rows: 12268 Columns: 24
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (9): blk_no, street, residential, commercial, market_hawker, miscellane...
## dbl (15): max_floor_lvl, year_completed, total_dwelling_units, 1room_sold, 2...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(X)

set.seed(1729) # replace the number '1729' with your matric no

Question 1

Change names of numeric variables to shorter ones and produce the pairs plots of all numeric variables in the dataset. Which variables have the highest correlation?

X <- X %>%
  rename(lvl = max_floor_lvl,
         year = year_completed,
         num_flats = total_dwelling_units,
         R1 = `1room_sold`,
         R2 = `2room_sold`,
         R3 = `3room_sold`,
         R4 = `4room_sold`,
         exec = exec_sold,
         MG = multigen_sold,
         studio = studio_apartment_sold)

X %>%
  select_if(is.numeric) %>%
  ggpairs(progress = F) + theme_minimal()

Let’s find out which variables have the highest absolute correlation:

X %>%
  select(where(is.numeric)) %>%
  cor %>%
  as_tibble(rownames = "variable_1") %>%
  pivot_longer(cols = -variable_1, 
               names_to = "variable_2", 
               values_to = "correlation") %>%
  filter(correlation < 1) %>%
  slice(which.max(abs(correlation)))

Question 2

Run principal component analysis on this dataset and print the matrix of loadings. What are interpretations of the first three principal components? Produce the scree plot that shows variance and cumulative variance explained by the principal components.

The matrix of loadings is

hdb_pca <- X %>%
  select_if(is.numeric) %>%
  prcomp(scale = TRUE)

hdb_pca$rotation %>% 
  as.data.frame() %>%
  select(1:3) %>%
  round(3)

The first principal component is the size of a building (it is mostly generated by the number of floors and the number of flats). The second principal component is mostly explained by the year of construction minus the number of 3-room apartments. It probably has to do with the type of architecture (older buildings have more 3-room apartments, which can be seen from the pairs plot). The third principal component shows if the building is made out of small apartments (2-room or studios), i.e., another aspect of the type of architecture.

Below is the scree plot of variance explained by principal components:

hdb_pca %>%
  tidy(matrix = "d") %>%
  pivot_longer(cols = percent:cumulative) %>%
  ggplot(aes(x = PC, y = value)) +
  geom_point() + geom_line() +
  facet_grid(.~ name) + theme_minimal() +
  labs(y = "Variance",
       title = "Variance explained by each principal component")

Question 3

Plot the data as a scatterplot in axes (PC1, PC2) coloured according to the town.

hdb_pca %>%
  augment(X) %>%
  ggplot(aes(x = .fittedPC1, y = .fittedPC2, color = bldg_contract_town)) +
  geom_point() + theme_minimal()

Question 4

Use \(K\)-means clustering with \(K=7\) to produce clusters for these HDB blocks. Print the contingency table of towns vs clusters. What is the most common HDB type (cluster) in Jurong West and what is the most common type in Queenstown? More generally, are Quenstown and Jurong West similar or different in terms of types of HDB blocks?

Here is the table of towns vs clusters:

hdb_clust <- X %>%
  select(where(is.numeric)) %>%
  scale %>%
  kmeans(7)

hdb_clust %>%
  augment(X) %>%
  tabyl(bldg_contract_town, .cluster)

To answer questions about Jurong West and Queenstown, it is more convenient to flip the table. Distributions of block clusters (aka types of blocks) in Jurong West and Queenstown are

hdb_clust_table <- hdb_clust %>%
  augment(X) %>%
  tabyl(.cluster, bldg_contract_town)

hdb_clust_table %>%
  select(.cluster, JW, QT)

The most common clusters are below:

hdb_clust_table %>%
  slice(which.max(JW)) %>%
  pull(.cluster) %>% 
  cat("The most common cluster in Jurong West is", ., "\n")

## The most common cluster in Jurong West is 2

hdb_clust_table %>%
  slice(which.max(QT)) %>%
  pull(.cluster) %>% 
  cat("The most common cluster in Queenstown is", .)

## The most common cluster in Queenstown is 5

Let’s convert raw numbers of clusters in Jurong West and Queenstown to frequencies so that we can compare them directly:

hdb_clust_table %>%
  select(.cluster, JW, QT) %>%
  mutate(across(c(JW, QT), ~. / sum(.)))

Generally, Jurong West and Queenstown are quite different in the composition of HDB blocks. To make sure that we are not mistaken, we can even run the \(\chi^2\)-test:

hdb_clust %>%
  augment(X) %>%
  filter(bldg_contract_town %in% c("JW", "QT")) %>%
  tabyl(.cluster, bldg_contract_town) %>%
  chisq.test()

## 
##  Pearson's Chi-squared test
## 
## data:  .
## X-squared = 183.89, df = 6, p-value < 2.2e-16

The \(p\)-value is nearly 0, which means that it is highly unlikely that the difference in block type decomposition is due to random chance.

Note that we could get fractions directly with tabyl():

hdb_clust %>%
  augment(X) %>%
  tabyl(.cluster, bldg_contract_town) %>%
  adorn_percentages("col") %>%
  select(.cluster, JW, QT)

Lab 8 — graded report.

Fedor Duzhin