Visualising High-dimensional Data with R

Dianne Cook
Monash University

Session 2

Outline

Show code

plan <- tribble(~time, ~topic,
                "30", "Understanding clusters in data using visualisation",
                "30", "Building better classification models with visual input")
knitr::kable(plan)

time	topic
30	Understanding clusters in data using visualisation
30	Building better classification models with visual input

Clustering

Method 1: Spin-and-brush

Show code

library(detourr)
set.seed(645)
detour(p_tidy_std[,2:5], 
       tour_aes(projection = bl:bm)) |>
       tour_path(grand_tour(2), fps = 60, 
                 max_bases=40) |>
       show_scatter(alpha = 0.7, 
                    axes = FALSE)

DEMO

What are clusters?

Ideal thinking of neatly separated clusters, but it is rarely encountered in data

Objective is to organize the cases into groups that are similar in some way. You need a measure of similarity (or distance).

Why visualise? Which is the better?

Nuisance cases

Show code

# Nuisance observations
set.seed(20190514)
x <- (runif(20)-0.5)*4
y <- x
d1 <- data.frame(x1 = c(rnorm(50, -3), 
                            rnorm(50, 3), x),
                 x2 = c(rnorm(50, -3), 
                            rnorm(50, 3), y),
                 cl = factor(c(rep("A", 50), 
                             rep("B", 70))))
d1 <- d1 |> 
  mutate_if(is.numeric, function(x) (x-mean(x))/sd(x))

# Nuisance variables
set.seed(20190512)
d2 <- data.frame(x1=c(rnorm(50, -4), 
                            rnorm(50, 4)),
                 x2=c(rnorm(100)),
                 cl = factor(c(rep("A", 50), 
                             rep("B", 50))))
d2 <- d2 |> 
  mutate_if(is.numeric, function(x) (x-mean(x))/sd(x))

d1_hs <- hclust(dist(d1[,1:2]),
                method="single")
d1 <- d1 |>
  mutate(cls = factor(cutree(d1_hs, 2)))
pc_d1s <- ggplot(d1) +
  geom_point(aes(x=x1, y=x2, colour=cls), 
             size=2, alpha=0.8) +
  scale_colour_discrete_divergingx(palette = "Zissou 1",
                                   nmax=4, rev=TRUE) +
  ggtitle("Result 1") +
  theme(legend.position="none")

d1_hw <- hclust(dist(d1[,1:2]),
                method="ward.D2")
d1 <- d1 |>
  mutate(clw = factor(cutree(d1_hw, 2)))
pc_d1w <- ggplot(d1) +
  geom_point(aes(x=x1, y=x2, colour=clw), 
             size=2, alpha=0.8) +
  scale_colour_discrete_divergingx(palette = "Zissou 1",
                                   nmax=4, rev=TRUE) +
  ggtitle("Result 2") +
  theme(legend.position="none")

pc_d1s + pc_d1w +
  plot_layout(ncol=2)

Nuisance variables

Show code

d2_hc <- hclust(dist(d2[,1:2]),
                method="complete")
d2 <- d2 |>
  mutate(clc = factor(cutree(d2_hc, 2)))
pc_d2c <- ggplot(d2) +
  geom_point(aes(x=x1, y=x2, colour=clc), 
             size=2, alpha=0.8) +
  scale_colour_discrete_divergingx(palette = "Zissou 1",
                                   nmax=4, rev=TRUE) +
  ggtitle("Result 1") + 
  theme(legend.position="none")

d2_hw <- hclust(dist(d2[,1:2]),
                method="ward.D2")
d2 <- d2 |>
  mutate(clw = factor(cutree(d2_hw, 2)))
pc_d2w <- ggplot(d2) +
  geom_point(aes(x=x1, y=x2, colour=clw), 
             size=2, alpha=0.8) +
  scale_colour_discrete_divergingx(palette = "Zissou 1",
                                   nmax=4, rev=TRUE) +
  ggtitle("Result 2") + 
  theme(legend.position="none")
  
pc_d2c + pc_d2w +
  plot_layout(ncol=2)

To decide on a best result, you need to see how it divides the data into clusters. The cluster statistics, like dendrogram, or cluster summaries, or gap statistics might all look good but the result is bad. You need to see the model in the data space!

Model-based clustering (1/3)

Model-based clustering fits a multivariate normal mixture model to the data.

\[ \Sigma_k = \lambda_kD_kA_kD_k^\top, ~~~k=1, \dots, g \]

where

\(\Sigma_k\) is the variance-covariance of cluster \(k\),

\(g=\)number of clusters,

\(D_k\) describes the orientation of a cluster,

\(A_k\) describes the variance in different variables,

\(\lambda_k\) is an overall size.

Model-based clustering (2/3)

Clustering this data. What do you expect?

Can we assume the shape of the clusters is elliptical?

Show code

p_BIC <- mclustBIC(p_tidy[,2:5])
ggmc <- ggmcbic(p_BIC, cl=2:9, top=7) + 
  scale_color_discrete_divergingx(palette = "Roma") 
ggmc

Volume, Shape, Orientation

Model-based clustering (3/3)

Show code

p_mc <- Mclust(p_tidy[,2:5], 
                      G=4, 
                      modelNames = "VEE")
p_mce <- mc_ellipse(ps_mc)
p_cl <- p_tidy
p_cl$cl <- factor(p_mc$classification)

p_mc_data <- p_cl |>
  select(bl:bm, cl) |>
  mutate(type = "data") |>
  bind_rows(bind_cols(p_mce$ell,
                      type=rep("ellipse",
                               nrow(p_mce$ell)))) |>
  mutate(type = factor(type))

Show code

animate_xy(p_mc_data[,1:4],
           col=p_mc_data$cl,
           pch=c(4, 20 )[as.numeric(p_mc_data$type)], 
           axes="off")

load("data/penguins_tour_path.rda")
render_gif(p_mc_data[,1:4], 
           planned_tour(pt1), 
           display_xy(col=p_mc_data$cl,
               pch=c(4, 20)[
                 as.numeric(p_mc_data$type)], 
                      axes="off",
               half_range = 0.7),
           gif_file="gifs/penguins_best_mc.gif",
           frames=500,
           loop=FALSE)

p_mc <- Mclust(p_tidy[,2:5], 
                      G=3, 
                      modelNames = "EEE")
p_mce <- mc_ellipse(p_mc)
p_cl <- p_tidy
p_cl$cl <- factor(p_mc$classification)

p_mc_data <- p_cl |>
  select(bl:bm, cl) |>
  mutate(type = "data") |>
  bind_rows(bind_cols(p_mce$ell,
                      type=rep("ellipse",
                               nrow(p_mce$ell)))) |>
  mutate(type = factor(type))

animate_xy(p_mc_data[,1:4],
           col=p_mc_data$cl,
           pch=c(4, 20)[as.numeric(p_mc_data$type)], 
           axes="off")

# Save the animated gif
load("data/penguins_tour_path.rda")
render_gif(p_mc_data[,1:4], 
           planned_tour(pt1), 
           display_xy(col=p_mc_data$cl,
               pch=c(4, 20)[
                 as.numeric(p_mc_data$type)], 
                      axes="off",
               half_range = 0.7),
           gif_file="gifs/penguins_simpler_mc.gif",
           frames=500,
           loop=FALSE)

Four-cluster VEE

Three-cluster EEE

Models (ellipses) are overlaid on the data. Which is the best fit?

How do you draw ellipses in high-d?

Extract the estimated model parameters

Show code

p_mc <- Mclust(
  p_tidy[,2:5], 
  G=3, 
  modelNames = "EEE")
p_mc$parameters$mean

   [,1] [,2] [,3]
bl   39   48   49
bd   18   15   18
fl  190  217  196
bm 3693 5076 3754

Show code

p_mc$parameters$variance$sigma[,,1]

      bl    bd     fl     bm
bl   8.4   1.6    8.3    755
bd   1.6   1.2    3.5    318
fl   8.3   3.5   42.5   1751
bm 754.7 318.4 1751.2 211467

Generate data that represents the ellipse(s) to overlay on the data.

Show code

p_mce <- mc_ellipse(p_mc)

Sample points uniformly on a pD sphere
Transform into an ellipse using the inverse variance-covariance matrix

Your turn

Use the spin-and-brush approach to extract the clusters from the c1 data.

Classification

What should you visualise?

Understand any clustering related to the known classes.
Obtain a sense of where boundaries might be placed.
Examine where the fitted model fits the data well, and where poorly.
Understand the misclassifications, whether they are reasonable given uncertainty in the data, or due to an ill-fitting or poorly specified model.
Understand what can happen with model fitting and pattern recognition with sparse data.

Example: Linear DA

Linear discriminant analysis is the ideal classifier for this data.

Show code

p_lda <- lda(species ~ ., p_tidy[,1:5], prior = c(1/3, 1/3, 1/3))
p_lda_pred <- predict(p_lda, p_tidy)
p_lda_pred_x1 <- data.frame(p_lda_pred$x)
p_lda_pred_x1$species <- p_tidy$species
ggplot(p_lda_pred_x1, 
       aes(x=LD1, y=LD2, 
           colour=species)) +
  geom_point(alpha=0.5) +  
  scale_color_discrete_divergingx("Zissou 1") +
#  scale_shape_manual(values=c(1, 2, 3)) +
  theme(legend.position = "bottom",
        legend.title = element_blank())

Random forests (1/2)

A random forest is the simplest classifier to fit for complicated boundaries. It is built from multiple trees generated by randomly sampling the cases and the variables. The random sampling (with replacement) of cases has the fortunate effect of creating a training (“in-bag”) and a test (“out-of-bag”) sample for each tree computed. The most beautiful results are obtaining diagnostics that help us to assess the model are the votes, the measure of variable importance, and the proximity matrix.

Show code

data(bushfires)

bushfires_sub <- bushfires[,c(5, 8:45, 48:55, 57:60)] |>
  mutate(cause = factor(cause))

bushfires_rf <- randomForest(cause~.,
                             data=bushfires_sub,
                             importance=TRUE)
bushfires_rf


Call:
 randomForest(formula = cause ~ ., data = bushfires_sub, importance = TRUE) 
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 7

        OOB estimate of  error rate: 11%
Confusion matrix:
            accident arson burning_off lightning
accident          73     3           0        62
arson             11     8           1        17
burning_off        3     0           3         3
lightning         14     0           0       823
            class.error
accident          0.471
arson             0.784
burning_off       0.667
lightning         0.017

Random forests (2/2)

The votes matrix can be considered to be predictive probabilities, where the values for each observation sum to 1. With 3 classes it is a 2D triangle. For 4 or more classes it is a simplex and can be examined in a tour.

Votes matrix for the random forest fit on penguins

Votes matrix for bushfire model fit

Show code

# Create votes matrix data
bushfires_rf_votes <- bushfires_rf$votes |>
  as_tibble() |>
  mutate(cause = bushfires_sub$cause)

# Project 4D into 3D
proj <- t(geozoo::f_helmert(4)[-1,])
b_rf_v_p <- as.matrix(bushfires_rf_votes[,1:4]) %*% proj
colnames(b_rf_v_p) <- c("x1", "x2", "x3")
b_rf_v_p <- b_rf_v_p |>
  as.data.frame() |>
  mutate(cause = bushfires_sub$cause)
  
# Add simplex
simp <- simplex(p=3)
sp <- data.frame(simp$points)
colnames(sp) <- c("x1", "x2", "x3")
sp$cause = ""
b_rf_v_p_s <- bind_rows(sp, b_rf_v_p) |>
  mutate(cause = factor(cause))
labels <- c("accident" , "arson", 
                "burning_off", "lightning", 
                rep("", nrow(b_rf_v_p)))

animate_xy(b_rf_v_p_s[,1:3], col = b_rf_v_p_s$cause, 
           axes = "off", half_range = 1.3,
           edges = as.matrix(simp$edges),
           obs_labels = labels)

Exploring misclassifications

Show code

library(crosstalk)
library(plotly)
library(viridis)
p_cl <- p_tidy_std |>
  mutate(pspecies = p_lda_pred$class) |>
  dplyr::select(bl:bm, species, pspecies) |>
  mutate(sp_jit = jitter(as.numeric(species), 0.5),
         psp_jit = jitter(as.numeric(pspecies), 0.5))
p_cl_shared <- SharedData$new(p_cl)

detour_plot <- detour(p_cl_shared, tour_aes(
  projection = bl:bm,
  colour = species)) |>
    tour_path(grand_tour(2), 
                    max_bases=50, fps = 60) |>
       show_scatter(alpha = 0.9, axes = FALSE,
                    width = "100%", height = "450px")

conf_mat <- plot_ly(p_cl_shared, 
                    x = ~psp_jit,
                    y = ~sp_jit,
                    color = ~species,
                    colors = viridis_pal(option = "D")(3),
                    height = 450) |>
  highlight(on = "plotly_selected", 
              off = "plotly_doubleclick") |>
    add_trace(type = "scatter", 
              mode = "markers")
  
bscols(
     detour_plot, conf_mat,
     widths = c(5, 6)
 )

DEMO

Your turn

Explore the misclassifications in the random forest fit of the penguins data, using the code provided in the slides2.R file.

Cautions about high-dimensions

Space is big.

What might appear to be structure is only sampling variability.

Show code

w <- matrix(runif(48*40), ncol=40) |>
  as.data.frame() |>
  mutate(cl = factor(rep(c("A", "B", "C", "D"), rep(12, 4))))
w_lda <- lda(cl~., data=w)
w_pred <- predict(w_lda, w, dimen=2)$x
w_p <- w |>
  mutate(LD1 = w_pred[,1],
         LD2 = w_pred[,2])
ggplot(w_p, aes(x=LD1, y=LD2, colour=cl)) + 
  geom_point() +
  scale_colour_discrete_divergingx(palette = "Zissou 1",
                                   nmax=4, rev=TRUE)

\(n=48, p=40\)

Show code

animate_xy(w[,1:40], guided_tour(lda_pp(w$cl)), 
           col=w$cl,
           sphere=TRUE, 
           axes="off",
           half_range=4)

Permutation is your friend, for high-dimensional data analysis.

Permute the class labels.

Show code

set.seed(951)
ws <- w |>
  mutate(cl = sample(cl))

Show code

ws_lda <- lda(cl~., data=ws)
ws_pred <- predict(ws_lda, ws, dimen=2)$x
ws_p <- ws |>
  mutate(LD1 = ws_pred[,1],
         LD2 = ws_pred[,2])
ggplot(ws_p, aes(x=LD1, y=LD2, colour=cl)) + 
  geom_point() +
  scale_colour_discrete_divergingx(palette = "Zissou 1",
                                   nmax=4, rev=TRUE)