---
title: "Unsupervised Learning with tidylearn"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Unsupervised Learning with tidylearn}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)
```

```{r setup}
library(tidylearn)
library(dplyr)
library(ggplot2)
```

## Two Ways In

tidylearn offers two entry points to the same algorithms, and which one you
want depends on what you are doing.

`tl_model(method = "kmeans")` puts clustering behind the same signature as
every supervised method, which is what you want when the clustering is a step
inside a larger workflow — see `vignette("integration-workflows")`.

The `tidy_*()` family is the fuller interface, and the subject of this
vignette. Each function returns a list of tibbles rather than a fitted object
you have to take apart, and each has a matching `augment_*()` that glues the
result back onto your data.

```{r}
# Same algorithm, two interfaces
model <- tl_model(iris[, 1:4], method = "kmeans", k = 3)
km <- tidy_kmeans(iris[, 1:4], k = 3)

names(km)
```

Both wrap `stats::kmeans()`; the algorithms are unchanged. Reach the raw
object through `model$fit$model` or `km$model` — an unsupervised `$fit`
is the list of tidied components, with the wrapped object among them.

**Wrapped packages:**

- stats (`prcomp()`, `kmeans()`, `hclust()`, `cmdscale()`)
- cluster (`pam()`, `clara()`)
- dbscan for density-based clustering
- MASS (`isoMDS()`, `sammon()`)
- smacof for MDS algorithms
- arules for association rules — see `vignette("market-basket")`

## Scale First

Every distance-based method here answers a question about distance, and
distance is measured in whatever units your columns happen to use. A variable
measured in thousands will dominate one measured in tenths, whatever its
actual relevance.

`standardize_data()` centres and scales:

```{r}
iris_scaled <- standardize_data(iris[, 1:4])

sapply(iris_scaled, function(x) round(c(mean = mean(x), sd = sd(x)), 3))
```

iris happens to have four variables on a similar scale, so the examples below
use the raw values and stay comparable to the species labels. On real data,
scale first.

## Principal Component Analysis

```{r}
pca <- tidy_pca(iris[, 1:4], scale = TRUE)
names(pca)
```

Accessors rather than list-digging:

```{r}
get_pca_variance(pca)
```

```{r}
get_pca_loadings(pca, n_components = 2)
```

Two components carry 96% of the variance. `plot_variance_explained()` marks
where a threshold is crossed:

```{r}
plot_variance_explained(get_pca_variance(pca), threshold = 0.9)
```

`tidy_pca_screeplot()` is the conventional scree plot:

```{r}
tidy_pca_screeplot(pca)
```

### Scores back on the data

`augment_pca()` returns the original data with component scores appended, so
the labels you already have stay attached:

```{r}
scored <- augment_pca(pca, iris, n_components = 2)
head(scored, 3)
```

```{r}
ggplot(scored, aes(x = PC1, y = PC2, color = Species)) +
  geom_point(size = 3, alpha = 0.7) +
  labs(
    title = "PCA of Iris",
    x = paste0("PC1 (",
               round(get_pca_variance(pca)$prop_variance[1] * 100, 1), "%)"),
    y = paste0("PC2 (",
               round(get_pca_variance(pca)$prop_variance[2] * 100, 1), "%)")
  ) +
  theme_minimal()
```

`tidy_pca_biplot()` overlays the loadings on the same scatter, which is how
you read what the components mean:

```{r}
tidy_pca_biplot(pca, color_by = iris$Species)
```

## How Many Clusters?

Guessing *k* and checking the answer against labels you happen to have is not
a method. `optimal_clusters()` runs three criteria at once.

```{r}
opt <- optimal_clusters(iris[, 1:4], max_k = 8)
names(opt)
```

```{r}
opt$wss
```

```{r}
opt$silhouette
```

```{r}
c(silhouette = attr(opt$silhouette, "optimal_k"),
  gap = opt$gap$recommended_k)
```

The three criteria disagree, which is normal. The elbow is a judgement call,
silhouette favours well-separated compact clusters, and the gap statistic
compares against a null of no structure. Silhouette says 2 here because
*versicolor* and *virginica* overlap; the botanical answer is 3.

```{r}
plot_elbow(opt$wss, suggested_k = 3)
```

```{r}
plot_gap_stat(opt$gap)
```

`calc_wss()` gives the within-cluster sums of squares on their own if that is
all you need:

```{r}
calc_wss(iris[, 1:4], max_k = 6)
```

## K-means

```{r}
km <- tidy_kmeans(iris[, 1:4], k = 3)
km$centers
```

```{r}
km$clusters
```

`augment_kmeans()` puts the assignment back on the data:

```{r}
iris_clustered <- augment_kmeans(km, iris)
table(Cluster = iris_clustered$cluster, Species = iris_clustered$Species)
```

```{r}
plot_cluster_sizes(km$clusters$cluster)
```

```{r}
plot_clusters(iris_clustered, cluster_col = "cluster",
              x_col = "Petal.Length", y_col = "Petal.Width")
```

### Was it a good clustering?

Silhouette width scores each observation on how much better it fits its own
cluster than the next-nearest one. Values near 1 are comfortable, near 0
borderline, negative means the point is on the wrong side.

```{r}
dist_mat <- tidy_dist(iris[, 1:4])
sil <- tidy_silhouette(km$clusters$cluster, dist_mat)

sil$avg_width
```

```{r}
sil$cluster_avg
```

```{r}
plot_silhouette(sil)
```

Cluster 1 is clean; the other two are the *versicolor*/*virginica* boundary,
and their scores say so without needing the labels.

`calc_validation_metrics()` collects the summary numbers in one row:

```{r}
calc_validation_metrics(km$clusters$cluster, iris[, 1:4], dist_mat)
```

## PAM and CLARA

PAM picks actual observations as cluster centres, which makes it less
sensitive to outliers than k-means and gives you a representative row rather
than an average.

```{r}
pam_result <- tidy_pam(iris[, 1:4], k = 3)
pam_result$medoids
```

```{r}
pam_result$silhouette_avg
```

```{r}
table(Cluster = augment_pam(pam_result, iris)$cluster, Species = iris$Species)
```

CLARA samples rather than computing the full distance matrix, which is what
makes it usable when *n* is large enough that an *n × n* matrix is not:

```{r}
large_data <- iris[rep(seq_len(nrow(iris)), 10), 1:4]
clara_result <- tidy_clara(large_data, k = 3, samples = 5)

table(clara_result$clusters$cluster)
```

## Hierarchical Clustering

`tidy_hclust()` builds the tree; cutting it is a separate decision.

```{r}
hc <- tidy_hclust(iris[, 1:4], method = "average")
plot_dendrogram(hc, k = 3)
```

`optimal_hclust_k()` scores cut heights the way `optimal_clusters()` scores
*k*:

```{r}
optimal_hclust_k(hc, method = "silhouette", max_k = 8)$optimal_k
```

```{r}
cuts <- tidy_cutree(hc, k = 3)
head(cuts, 3)
```

```{r}
hc_data <- augment_hclust(hc, iris, k = 3)
table(Cluster = hc_data$cluster, Species = hc_data$Species)
```

Linkage matters more than most people expect. `"average"`, `"complete"`,
`"single"` and `"ward.D2"` can produce different trees from the same distances:

```{r}
linkages <- c("single", "average", "complete", "ward.D2")

sapply(linkages, function(m) {
  cl <- tidy_cutree(tidy_hclust(iris[, 1:4], method = m), k = 3)$cluster
  max(table(cl))
})
```

Single linkage chains, so it puts almost everything in one cluster. That is a
property of the linkage, not a finding about irises.

## DBSCAN

DBSCAN finds arbitrarily shaped clusters and labels sparse points as noise.
It needs `eps` (the neighbourhood radius) and `minPts`. Rather than guessing,
`suggest_eps()` reads it off the k-nearest-neighbour distance curve.

```{r}
eps_suggestion <- suggest_eps(iris[, 1:4], minPts = 5)
eps_suggestion$eps
```

```{r}
plot_knn_dist(iris[, 1:4], k = 5)
```

The elbow in that curve is where points stop having close neighbours, which
is the radius you want.

```{r}
db <- tidy_dbscan(iris[, 1:4], eps = eps_suggestion$eps, minPts = 5)

c(clusters = db$n_clusters, noise = db$n_noise)
```

```{r}
db_data <- augment_dbscan(db, iris)
table(Cluster = db_data$cluster, Species = db_data$Species)
```

Cluster 0 is noise, not a cluster.

`explore_dbscan_params()` sweeps the two parameters together, which is more
informative than tuning either alone:

```{r}
explore_dbscan_params(
  iris[, 1:4],
  eps_values = c(0.4, 0.6, 0.8, 1.0),
  minPts_values = c(4, 5, 10)
)
```

Read `prop_noise` alongside `n_clusters`: a setting that finds many clusters
by discarding a third of the data has not found structure.

## Multidimensional Scaling

MDS places observations so that their plotted distances reproduce their
distances in the original space. Unlike PCA it can work from any distance
matrix, including non-Euclidean ones.

```{r}
mds <- tidy_mds(iris[, 1:4], method = "classical", ndim = 2)
head(mds$config, 3)
```

```{r}
plot_mds(mds, color_by = iris$Species, label_points = FALSE)
```

`method` also takes `"metric"` and `"nonmetric"` (both via smacof),
`"sammon"` and `"kruskal"`. The last two minimise a stress function by
dividing through the observed distances, so they need every pairwise distance
to be strictly positive. iris contains one duplicated row, which is enough to
stop them:

```{r, error = TRUE}
tidy_mds(iris[, 1:4], method = "sammon", ndim = 2)
```

Drop the duplicates and they run:

```{r}
distinct_iris <- iris[!duplicated(iris[, 1:4]), 1:4]
sammon <- tidy_mds(distinct_iris, method = "sammon", ndim = 2)
sammon$stress
```

Stress is the number to check before reading anything into a non-metric
layout: below about 0.05 the picture is a faithful rendering of the
distances, and above about 0.2 it is decoration.

## Comparing Clusterings

`compare_clusterings()` scores several partitions of the same data side by
side.

```{r}
comparison <- compare_clusterings(
  list(
    kmeans = km$clusters$cluster,
    pam = pam_result$clusters$cluster,
    hclust = cuts$cluster,
    dbscan = db$clusters$cluster
  ),
  iris[, 1:4],
  dist_mat
)

comparison
```

```{r}
plot_cluster_comparison(
  iris[, 1:4] %>%
    mutate(kmeans = km$clusters$cluster, hclust = cuts$cluster),
  cluster_cols = c("kmeans", "hclust"),
  x_col = "Petal.Length",
  y_col = "Petal.Width"
)
```

The distance metric is a choice too. `compare_distances()` computes several
so you can see whether your conclusion depends on it:

```{r}
names(compare_distances(iris[, 1:4]))
```

```{r}
plot_distance_heatmap(dist_mat)
```

Look for block structure on the diagonal — contiguous runs of small
distances are what a cluster looks like here.

## A Worked Sequence

Putting the pieces in the order they belong:

```{r}
data_matrix <- standardize_data(iris[, 1:4])

# 1. How many clusters does the data support?
choice <- optimal_clusters(data_matrix, max_k = 8)
k <- attr(choice$silhouette, "optimal_k")
k
```

```{r}
# 2. Cluster at that k
final_km <- tidy_kmeans(data_matrix, k = k)

# 3. Score the result before believing it
final_sil <- tidy_silhouette(final_km$clusters$cluster, tidy_dist(data_matrix))
final_sil$avg_width
```

```{r}
# 4. Attach the assignment and look at it
final_data <- augment_kmeans(final_km, iris)
table(Cluster = final_data$cluster, Species = final_data$Species)
```

```{r}
plot_clusters(final_data, cluster_col = "cluster",
              x_col = "Petal.Length", y_col = "Petal.Width")
```

Silhouette chose 2, and the table shows what that means: *setosa* separated,
the other two species merged. That is the honest reading of this data at this
metric — the third cluster exists botanically but is not well separated in
these four measurements.

## Function Reference

| Task | Function |
|---|---|
| Prepare | `standardize_data()`, `tidy_dist()`, `compare_distances()` |
| Choose *k* | `optimal_clusters()`, `calc_wss()`, `tidy_gap_stat()`, `optimal_hclust_k()` |
| Reduce | `tidy_pca()`, `tidy_mds()` |
| Cluster | `tidy_kmeans()`, `tidy_pam()`, `tidy_clara()`, `tidy_hclust()`, `tidy_dbscan()` |
| Attach results | `augment_pca()`, `augment_kmeans()`, `augment_pam()`, `augment_hclust()`, `augment_dbscan()` |
| Validate | `tidy_silhouette()`, `calc_validation_metrics()`, `compare_clusterings()` |
| Tune DBSCAN | `suggest_eps()`, `plot_knn_dist()`, `explore_dbscan_params()` |
| Plot | `plot_elbow()`, `plot_silhouette()`, `plot_gap_stat()`, `plot_clusters()`, `plot_cluster_sizes()`, `plot_dendrogram()`, `plot_mds()`, `plot_distance_heatmap()`, `plot_variance_explained()` |
| Accessors | `get_pca_variance()`, `get_pca_loadings()`, `tidy_cutree()` |

Association rules have their own vignette: `vignette("market-basket")`.