---
title: "🗓️ Week 09: Anomaly Detection"
subtitle: "Theme: Unsupervised Learning"
author: "Dr. Ghita Berrada"
format: 
  html:
    from: markdown+emoji
    page-layout: full
    toc: true
    toc-depth: 2
self-contained: true
---

Welcome to Week 09! This week, we finish our exploration of unsupervised learning techniques with anomaly...

# :gear: Setup

Before you can load the libraries below, you might have to do the following package installs (feel free to skip if you already have done the installs before)

``` r
install.packages("fpc")
install.packages("factoextra")
install.packages("dbscan")
install.packages("tidytext")
install.packages("ClusterR")
install.packages("isotree")
install.packages("embed")
```

Then, let's load the libraries we'll need today.

```{r, echo=FALSE, message=FALSE, warning=FALSE}
# Tidyverse packages
library(tidyr)
library(dplyr)
library(ggplot2)
library(readr)
library(skimr)
library(tidytext)
library(viridis)

# Tidymodels packages, spelt out individually for clarity
library(recipes)

#new libraries for clustering
library(fpc) #library for Calinski-Harabasz index
library(factoextra) # library for elbow and average silhouette plots
library(dbscan) # library for DBSCAN
library(ClusterR)
library(cluster)

#library for umap
library(embed)

#library for isolation forest
library(isotree)

```

## Let's start with DBSCAN first!

Last week, we stopped by saying that k-means (and its variants) were not algorithms that were suitable for use in all (clustering) situations.

Let's see why again: let's check the pimpled smiley face example [here](https://www.naftaliharris.com/blog/visualizing-k-means-clustering/). Does K-Means capture the expected clusters?

Now let's check the same example with [DBSCAN](https://www.naftaliharris.com/blog/visualizing-dbscan-clustering/).

#### How does DBSCAN work?

![](./figures_W09/dbscan-principle.png)

Fig 1: Different types of points and parameters in DBSCAN

*A few definitions to start with*

-   Any point $x$ in the data set, with a neighbour count greater than or equal to `MinPts`, is a **core point**
-   $x$ is a **border point**, if the number of its neighbors is less than `MinPts`, but it belongs to the $\epsilon$-neighborhood of some core point $z$.
-   If a point is neither a core nor a border point, then it is called a **noise point** or an **outlier**.

E.g: The figure (Fig. 1) shows the different types of points (core, border and outlier points) using `MinPts = 6`:

-   $x$ is a core point because $neighbours_\epsilon(x)=6$

-   $y$ is a border point because $neighbours_\epsilon(y)<MinPts$, but it belongs to the $\epsilon$-neighborhood of the core point $x$.

-   Finally, $z$ is a noise point.

*A few more definitions:*

-   **Direct density reachable**: A point “A” is directly density reachable from another point “B” if: i) “A” is in the $\epsilon$-neighborhood of “B” and ii) “B” is a core point. E.g $y$ is directly density reachable
-   **Density reachable**: A point “A” is density reachable from “B” if there are a set of core points leading from “B” to “A.
-   **Density connected**: Two points “A” and “B” are density connected if there are a core point “C”, such that both “A” and “B” are density reachable from “C”.

In density-based clustering, a cluster is defined as a group of density connected points.

### The DBSCAN Algorithm

1.  For each point $x_i$, compute the distance between $x_i$ and the other points. Find all neighbor points within distance `eps` of the starting point ($x_i$). Each point, with a neighbor count greater than or equal to `MinPts`, is marked as core point or visited.

2.  For each core point, if it’s not already assigned to a cluster, create a new cluster. Find recursively all its density connected points and assign them to the same cluster as the core point.

3.  Iterate through the remaining unvisited points in the data set.

The points that do not belong to any cluster are treated as outliers or noise.

If this explanation confuses you, you can have a look at [this Youtube video](https://www.youtube.com/watch?app=desktop&v=RDZUdRSDOok)

### How would we apply DBSCAN to the quake example?

```{r}
quakes_features<- recipe(~ ., data = quakes) %>%
  step_rm(stations) %>% 
  step_normalize(all_predictors()) %>% 
  prep() %>% 
  bake(new_data = NULL)

normalized_quakes_df <- 
  quakes_features %>% 
  select(depth,mag) 

df<- normalized_quakes_df %>% select(where(is.numeric))
db <- dbscan::dbscan(df, eps = 0.26, minPts = 4)

# Plot DBSCAN results
fviz_cluster(db, data = df, stand = FALSE,
             ellipse = FALSE, show.clust.cent = FALSE,
             geom = "point", ggtheme = theme_bw(base_size=20))
```

What do you make of this plot? What do you think the black dots in it represent?

For a quick reminder, the result with k-means clustering was as follows:

```{r}
normalized_fitted_kmeans <- kmeans(normalized_quakes_df, centers = 4, iter.max = 10, nstart = 5)
fviz_cluster(normalized_fitted_kmeans, data = normalized_quakes_df, stand = FALSE,
             ellipse = FALSE, show.clust.cent = FALSE,
             geom = "point", ggtheme = theme_bw(base_size=20))
```

### How do we choose the `eps` parameter in DBSCAN

Now that we now how DBSCAN works, the million dollar question is how do we set its parameters? Well, once again, it's rather tricky...

The value of `minPts` should be at least one greater than the number of dimensions of the dataset, i.e., **`minPts>=Dimensions+1`**. It does not make sense to take `minPts` as 1 because the result will be each point being a separate cluster. Therefore, it must be at least 3. Generally, we go for at least twice the dimensions, e.g if we only only select two variables in our data (i.e `depth` and `mag` ), we'd select `minPts` as 4.

Now that `minPts` is set, we can select `eps` using the k-distance graph. In practice, this means that:

-   you compute the average distance (by default Euclidean distance) between each point and its k nearest neighbors, where k = the `minPts` value you selected.

-   the average k-distances are then plotted in ascending order on a k-distance graph. You’ll find the optimal value for `eps` at the point of maximum curvature (i.e. where the graph has the greatest slope).

    The code below gives an example of the k-distance graph for our ongoing quakes clustering example...

```{r}
dbscan::kNNdistplot(df, k = 4)+
abline(h = 0.26, lty = 2) +

scale_y_continuous(n.breaks = 40)     
```

As we can see, the point of maximum curvature is around 0.26, which sets our `eps` value...

### Quick note: Similarity metrics

Underlying the discussion about clustering, are similarity metrics. We mostly assumed until now that the similarity metric that was used was the Euclidean distance (and for k-means, this is a fair assumption to make). But there are cases where the Euclidean distance would not be an appropriate choice (e.g when comparing strings or sets or for computer vision or recommendation engines applications).

Here is a quick run-through of some of the most common similarity metrics:

*Real-valued data:*

-   Euclidian distance: $d(A,B)=\sqrt{\sum_{i=1}^n |A_i-B_i|^2}$

-   Cosine distance: $S_C(A,B)=\frac{\sum_{i=1}^n A_iB_i}{\sqrt{\sum_{i=1}^n A_i^2}\sqrt{\sum_{i=1}^n B_i^2}}$

where **A** and **B** are two n-dimensional vectors of attributes, and $A_i$ and $B_i$ the $i$th components of vectors **A** and **B** respectively

(metric frequently used in NLP and computer vision - see \[this page\](<https://www.datastax.com/guides/real-world-applications-of-cosine-similarity>) for more information)

*String similarity:*

-   Hamming distance: given two strings of equal length, it is the number of positions at which the corresponding symbols (i.e either letters, bits, or decimal digits) are different e.g:

    -   [0000]{style="color:#BA160C"} and [1111]{style="color:#BA160C"} is 4
    -   ca[t]{style="color:#BA160C"} and ca[r]{style="color:#BA160C"} is 1
    -   c[a]{style="color:#BA160C"}l[l]{style="color:#BA160C"} and c[o]{style="color:#BA160C"}l[d]{style="color:#BA160C"} is 2

-   Levenstein distance: given two strings, it's the minimal number of edits (substitutions, deletions, insertions) needed to turn one string into another e.g:

    The Levenshtein distance between *kitten* and *sitting* is 3, since the following 3 edits change one into the other, and there is no way to do it with fewer than 3 edits:

    1.  [k]{style="color:#BA160C"}itten → [s]{style="color:#BA160C"}itten (substitution of "s" for "k"),
    2.  sitt[e]{style="color:#BA160C"}n → sitt[i]{style="color:#BA160C"}n (substitution of "i" for "e"),
    3.  sittin → sittin[g]{style="color:#BA160C"} (insertion of "g" at the end).

*Set similarity:*

A set is a collection of items with no order or repetition. Sets are typically used to represent relationships or associations between objects or even people.

-   Jaccard index: it is often used in recommendation systems and social media analysis and measures the similarity between two sets based on the number of items present in both sets relative to the total number of items:

    $J(A,B)=\frac{|A\cap B|}{|A\cup B|}$

    where A and B are two sets, $|A\cap B|$ the number of items included in both sets and $|A\cup B|$ the total number of items

Another useful resources: \[this page\](<https://scikit-learn.org/stable/modules/clustering.html>) compares various clustering algorithms (among which k-means and DBSCAN) and shows their parameters, use cases,...

## Part 2: Anomaly detection

### What are anomalies?

Anomalies are data points that deviate significantly from the expected or “normal behavior.”

We generally distinguish between univariate and multivariate outliers/anomalies:

-   **univariate outliers** exist in a single feature in isolation: they are extreme or abnormal values that deviate from the typical range of values for that specific feature. Say that the temperatures in London in December are typically between 2 and 7 degrees. A temperature of -10 degree in December would be considered a univariate outlier/anomaly.

-   **Multivariate outliers** are found by combining the values of multiple variables at the same time: taken separately, the values for each variable appear normal but in conjunction, they are not. Take, for example, a month of December in London where the temperatures are persistently on the higher end of the usual range (i.e close to 7 degrees) and where there were fewer than 5 days of rain: though the temperatures were in the normal range and it rained as usual, it is unusual for the temperatures to consistently stay on the higher end of the scale for such a long period and for it to rain for a short period, and therefore, taken cumulatively, the pattern for December is clearly an outlier!

### Running example for this part: Indicators of Heart disease

Let's introduce our main example for this part.

This dataset is a dataset sourced from [Kaggle](https://www.kaggle.com/datasets/kamilpytlak/personal-key-indicators-of-heart-disease), originates from an 2022 annual CDC survey data of 400k+ adults related to their health status but more specifically focuses on indicators of heart disease.

```{r}
heart_df <- read_csv("data/heart_2020_cleaned.csv")
skim(heart_df)
```

No missing values! 4 numeric columns and 14 categorical columns.

To make full use of this data for clustering/anomaly detection, we would so far either discard the categorical features or use a `step_dummy` to convert categorical features to numerical ones.

Let's convert the categorical to numeric features.

```{r}

heart_slice <- heart_df %>% slice_sample(prop = 0.02) %>%
               mutate(across(where(is.character), as.factor))# slicing the data to allow computations to fit within memory!

heart_transformed<- recipe(~ ., data = heart_slice) %>%
  step_dummy(all_nominal_predictors(), one_hot = TRUE) %>% 
  step_normalize(all_predictors()) %>% 
  prep() %>% 
  bake(new_data = NULL)

head(heart_transformed)
skim(heart_transformed)
```

I now have 52 features!

### Anomaly detection with DBSCAN

We've already seen earlier (in the quakes example) that DBSCAN simply treats outliers or anomalies as a "separate cluster", that points that are not core or border points are labeled as outliers/anomalies and that therefore the number/extent of anomalies directly relates to the DSBCAN parameters chosen (i.e `eps` and `minPts`).

If you want to find out more about what to do when you get new data and about other clustering-based anomaly detection techniques, head [this way](https://towardsdatascience.com/best-clustering-algorithms-for-anomaly-detection-d5b7412537c8).

Let's go back to our indicators of heart disease example.

First, let's see if we can determine optimal parameters for DBSCAN.

We have 52 features in our dataset so we need to choose `minPts` at least equal to 104.

Then, let's determine `eps` :

```{r}
dbscan::kNNdistplot(heart_transformed, k = 104)+
  abline(h = 10.2, lty = 2) +
scale_y_continuous(n.breaks = 40)   
```

The point of maximum curvature is around 10.2, which sets our `eps` value...

Let's now run DBSCAN.

```{r}
heart_dbscan<-dbscan::dbscan(heart_transformed,eps=10.2,minPts = 104)

```

Let's visualise the results:

```{r}
heart_transformed <- heart_transformed %>%
  mutate(cluster = as.factor(heart_dbscan$cluster)) # I associate each point with their cluster label

# Step 2: Create a recipe with UMAP
umap_recipe <- recipe(cluster ~ ., data = heart_transformed) %>%
  step_umap(all_numeric_predictors(), num_comp = 2) %>% # Reduce to 2 dimensions for visualization with UMPAP
  prep() %>%
  bake(new_data=NULL)

# Step 4: Visualize with ggplot2
ggplot(umap_recipe, aes(UMAP1, UMAP2, color = cluster)) +
  geom_point(alpha = 0.7, size = 1.5) +
  theme_minimal() +
  labs(title = "UMAP Visualization of DBSCAN Clusters",
       x = "UMAP Dimension 1",
       y = "UMAP Dimension 2",
       color = "Cluster")
```

### Anomaly detection with Isolation Forest

#### How Isolation Forest works:

1.  One feature within the dataset is selected

2.  A split value within the range of the selected feature is randomly selected for the instances currently in the node.

3.  The data is into two child nodes based on the split value selected in step 2.

4.  The process from steps 1 to 3 is repeated recursively for each child node until one of the following conditions is met:

-   Each leaf node has only one instance.

-   A predefined maximum depth is reached.

So, basically one tree (called isolation tree) is created per feature.

How do we then calculate anomaly scores and determine the anomalies?

The basic idea is that we account for the length of the path for each observation within a tree (we call that $h(x)$ ).

The shorter that path is, the more likely that observation is to be anomalous for that particular feature.

We could be tempted to simply take the average of path lengths (called $E[h(x)]$) over all isolation trees to get our anomaly score: however, not all features and paths are directly comparable.

So, we need a normalization term: it's based on sample size $n$ (and Binary Search Tree theory)

$$
c(n)=\begin{cases}
     2H(n-1)-2\frac{n-1}{n} \textrm{ if } n>2\\
     1 \textrm{ if } n=2\\
     0 \textrm{ otherwise}
     \end{cases}
$$

where $H(i)$ is the harmonic number that can be estimated by $ln(i) + 0.5772156649$ (Euler’s constant)

and the anomaly score is:

$$
s(x,n)=2^{-\frac{E(h(x))}{c(n)}}
$$

​The score ranges between 0 and 1, and when interpreting it, there are three scenarios:

-   If a point is an anomaly then $E[h(x)] < c(n)$ and $s → 1$

-   If a point is not an anomaly then $E[h(x)] > c(n)$ and $s → 0$ 

-   If $E[h(x)] = c(n)$ then $s = 0.5$

In other words, when the average path length of an instance is shorter than the expected average path length, the anomaly score will be closer to 1.

Any score lower than 0.5 will be identified as normal instances.

To read more about Isolation Forest, you can refer to the original paper: Fei Tony Liu, Kai Ming Ting and Zhi-Hua Zhou, ["Isolation Forest,"](https://ieeexplore.ieee.org/document/4781136) *2008 Eighth IEEE International Conference on Data Mining*, Pisa, Italy, 2008, pp. 413-422, doi: 10.1109/ICDM.2008.17

How do we implement Isolation Forest on the earthquakes data?

```{r}
iso_model <- isolation.forest(
  normalized_quakes_df,
  ndim = 1,          # Number of random projections per tree
  ntrees = 100,      # Number of trees
  sample_size = 256  # Subsample size for building trees
)
```

How do you set the parameters `ntrees`,`ndim` and `sample_size`?

-   `ntrees` is number of Isolation Trees used in the ensemble. A value of `100` is used in the Isolation Forest paper. Through experimentation, the researchers found this to produce good results over a variety of datasets. This would generally be the **default** value and would be particularly recommended for datasets with a number of observations lower than 10000. For larger datasets, you can increase `ntrees` to 200-300. The idea is that more trees improve the stability of anomaly scores but increase computation time. You can try and tune this parameter to check if more trees significantly improve results.

-   `sample_size` is the number of observations used to build each tree. The idea is that smaller samples isolate anomalies more effectively. The **default** value is 256 (or the total number of observations in the dataset if the dataset has fewer than 256 observations).

-   `ndim` is the number of features randomly selected for splitting at each node. By default, it is set to either 1 (e.g you only have 2 features) or 2: you should use $\sqrt{\text{number of features}}$ as a rule of thumb for `ndim` (`ndim` can be increased to 2 or 3 for faster splits while maintaining randomness of splits in higher-dimensional data). You can keep `ndim` (e.g 1 or 2) low if the features are independent to allow for simpler splits and increase it.

```{r}
anomaly_scores <- predict(iso_model, normalized_quakes_df, type = "score")

# Add the scores to the dataset
normalized_quakes_df <- normalized_quakes_df %>%
mutate(anomaly_score = anomaly_scores)
```

```{r}
threshold <- 0.5 #anomaly score threshold

# Add an anomaly flag to the dataset
normalized_quakes_df <- normalized_quakes_df %>%
  mutate(is_anomaly = anomaly_score > threshold)
```

We plot our data with observations colored by anomaly score.

```{r}
ggplot(normalized_quakes_df, aes(x = depth, y = mag, color = anomaly_score)) +
  geom_point(size = 2, alpha = 0.7) +
  scale_color_viridis_c() +  # Viridis color gradient
  theme_minimal() +
  labs(
    title = "Isolation Forest Anomaly Detection",
    x = "Normalized Depth",
    y = "Normalized Magnitude",
    color = "Anomaly Score"
  )
```

We also plot our data, highlighting anomalies (based on the threshold we set!).

```{r}
ggplot(normalized_quakes_df, aes(x = depth, y = mag, color = is_anomaly)) +
  geom_point(size = 2, alpha = 0.7) +
  scale_color_manual(values = c("FALSE" = "cornflowerblue", "TRUE" = "indianred")) +
  theme_minimal() +
  labs(
    title = "Anomaly Detection Results (Isolation Forest)",
    x = "Normalized Depth",
    y = "Normalized Magnitude",
    color = "Anomaly (Threshold=0.5)"
  )

```

#### **How do you determine the anomaly threshold?**

We previously set the anomaly threshold to a 0.5 threshold, meaning that any point with a score higher than 0.5 would be flagged as an anomaly.

More generally, the anomaly threshold determines whether a point is flagged as anomalous. How do we set it? There are several approaches to set it:

**1. Default Heuristic**

-   In isolation forests, an anomaly score \>0.5 generally indicates some level of anomaly.

-   You can start with a threshold of **0.6–0.7** for the anomaly score and adjust the threshold based on the fraction of data flagged as anomalies.

**2. Use Percentiles of Anomaly Scores**

Anomaly scores fall within the range \[0, 1\] where scores closer to 1 indicate a higher likelihood of being an anomaly.

Here is how you set the anomaly threshold using percentiles:

1.  Calculate the anomaly scores for the dataset.

2.  Use a high percentile (e.g., 95th or 99th percentile) as the threshold (points with scores above this threshold are flagged as anomalies).

    ```{r}
    threshold <- quantile(anomaly_scores, 0.95)  # 95th percentile
    threshold
    ```

Based on this method, you'd be setting the threshold to 0.6 (you could obviously choose a lower percentile if you want to capture more anomalies).

**3. Visual Inspection**:

-   Plot a histogram of the anomaly scores:

    ```{r}
    hist(anomaly_scores, breaks = 30, col = adjustcolor("cornflowerblue", alpha.f = 0.6), main = "Anomaly Scores")
    ```

-   Identify natural separations or clusters in the scores and set the threshold accordingly

Looking at the histogram above for our data, most scores are clustered between 0.4 and 0.5 but we also have a second smaller cluster between around 0.52 to 0.59. It looks relatively safe to say that anomalies are points whose score is above 0.6.

#### **4. Based on Domain Knowledge**

-   Use domain expertise to decide the acceptable anomaly level.

    -   For example, if anomalies occur in 1% of observations, set the threshold such that \~1% of the points are flagged as anomalies.

:house_with_garden: **TAKE-HOME EXERCISE**: Can you do the same with the Heart Disease indicators data?

### Anomaly detection with LOF (Local Outlier Factor)

Both Isolation Forest and DBSCAN are methods that detect global outliers. Local Outlier Factor, in contrast, identifies local outliers.

Let's start with a few definition again:

-   K-distance is the distance between a point and it's K$^{th}$ nearest neighbor.
-   K-neighbors, denoted by $N_k(A)$ is the set of points that lie in or on the circle of radius K-distance. K-neighbors can be greater or equal to the value of K.

Let's see an example:

![](./figures_W09/k_distance.png){width="60%"}

Let's say we have four points: A, B, C, and D.

If K=2, K-neighbors of A include C, B and D and therefore $||N_2(A)||$ = 3.

-   Reachability distance:

    -   Given two points $x_j$ and $x_i$, it is defined as the maximum of K-distance of $x_j$ and the distance between $x_i$ and $x_j$. The distance measure is problem-specific (e.g Euclidean, Manhattan).

    Written more formally, this would look like this: $RD(x_i,x_j)=\max(\textrm{K-distance}(x_j),\textrm{distance}(x_i,x_j))$

    Said differently, if a point $x_i$ belongs to the K-neighbors of $x_j$ , the reachability distance will be K-distance of $x_j$ (blue line), otherwise the reachability distance will be the distance between $x_i$ and $x_j$ (orange line)

![](./figures_W09/reachability_LOF.png){width="60%"}

-   Local reachability density: It is inverse of the average reachability distance of a point A from its neighbors. The exact formula is:

    $LRD_k(A)=\frac{1}{\sum_{x_j\in N_k(A)} \frac{RD(A,x_j)}{||N_k(A)||}}$

    Intuitively according to the LRD formula, the higher the average reachability distance (i.e., neighbors are far from the point), the smaller the density of points are present around a particular point. This tells us how far a point is from the nearest cluster of points. Low LRD value imply that the closest cluster is far from the point.

    In other words, the LRD relies on the intuition that the longer the distance to the next neighbors, the sparser the area the respective point is located in. It tells us how far we have to travel from our point to reach the next point or cluster of points. So the lower it is, the less dense it is, the longer we have to travel.

-   The LRD of each point is compared with the average LRD of its K-neighbors. We compute the local outlier factor or LOF which is the ratio of between the average LRD of the K-neighbors of point A and the LRD of A.

Or more specifically: $LOF_k(A)=\frac{\sum_{x_j\in N_k(A)} LRD_k(x_j)}{||N_k(A)||}\times\frac{1}{LRD_k(A)}$

-   Intuitively, if the point is not an outlier (i.e is an inlier), the ratio of average LRD of neighbors is approximately equal to the LRD of the point (because the density of a point and its neighbors are roughly equal). In that case, LOF is nearly equal to 1. On the other hand, if the point is an outlier, the LRD of a point is less than the average LRD of neighbors. Then LOF value will be high.

-   Generally, if LOF\> 1, the point is considered an outlier, but that is not always true. Let’s say we know that we only have one outlier in the data, then we take the maximum LOF value among all the LOF values, and the point corresponding to the maximum LOF value will be considered as an outlier.

![](./figures_W09/LOF_example.png){width="65%"}

-   4 points: A(0,0), B(1,0), C(1,1) and D(0,3) and K=2.
-   We will use LOF to detect one outlier (D) among these 4 points.
-   First, calculate the K-distance, distance between each pair of points, and K-neighborhood of all the points with K=2
-   We will be using the Manhattan distance as a measure of distance i.e $d(x,y)=\sum_{i=1}^n |x_i-y_i|$

| Distance | Value |
|----------|-------|
| d(A,B)   | 1     |
| d(A,C)   | 2     |
| d(A,D)   | 3     |
| d(B,C)   | 1     |
| d(B,D)   | 4     |
| d(C,D)   | 3     |

K-distance(A): since C is the 2$^{nd}$ nearest neighbor of A, K-distance(A)=d(A,C) =2

K-distance(B): since A, C are the 2$^{nd}$ nearest neighbor of B, K-distance(B)=distance(B,C) OR K-distance(B)=distance(B,A) = 1

K-distance(C): since A is the 2$^{nd}$ nearest neighbor of C –\> K-distance(C)=distance(C,A) =2

K-distance(D): since A,C are the 2$^{nd}$ nearest neighbor of D –\> K-distance(D)=distance(D,A) or K-distance(D)=distance(D,C) =3

-   K-neighborhood (A) = {B,C} , $||N_2(A)|| =2$
-   K-neighborhood (B) = {A,C}, $||N_2(B)|| =2$
-   K-neighborhood (C)= {B,A}, $||N_2(C)|| =2$
-   K-neighborhood (D) = {A,C}, $||N_2(D)|| =2$

![](./figures_W09/LOF_example.png){width="65%"}

We use the K-distance and K-neighbourhood to compute the LRD for each point

$LRD_2(A)=\frac{1}{\frac{RD(A,B)+RD(A,C)}{||N_2(A)||}}=\frac{\frac{1}{1+2}}{2}=0.667$ $LRD_2(B)=\frac{1}{\frac{RD(B,A)+RD(B,C)}{||N_2(B)||}}=\frac{\frac{1}{2+2}}{2}=0.50$ $LRD_2(C)=\frac{1}{\frac{RD(C,A)+RD(C,B)}{||N_2(C)||}}=\frac{\frac{1}{1+2}}{2}=0.667$ $LRD_2(D)=\frac{1}{\frac{RD(D,A)+RD(D,C)}{||N_2(D)||}}=\frac{\frac{1}{3+3}}{2}=0.337$

$LOF_2(A)=\frac{LRD_2(B)+LRD_2(C)}{||N_2(A)||}\times\frac{1}{LRD_2(A)}=\frac{0.5+0.667}{2}\times\frac{1}{0.667}=0.87$ $LOF_2(B)=\frac{LRD_2(A)+LRD_2(C)}{||N_2(B)||}\times\frac{1}{LRD_2(B)}=\frac{0.667+0.667}{2}\times\frac{1}{0.5}=1.334$ $LOF_2(C)=\frac{LRD_2(A)+LRD_2(B)}{||N_2(C)||}\times\frac{1}{LRD_2(C)}=\frac{0.667+0.5}{2}\times\frac{1}{0.667}=0.87$ $LOF_2(D)=\frac{LRD_2(A)+LRD_2(C)}{||N_2(D)||}\times\frac{1}{LRD_2(D)}=\frac{0.667+0.667}{2}\times\frac{1}{0.337}=2$

The highest LOF among the four points is LOF(D). Therefore, D is an outlier/anomaly.

### Local outlier factor (LOF): Pros and Cons

*Pros*

-   LOF identifies outliers based on the local neighborhood
-   A point will be considered as an anomaly/outlier if it is at a short distance from an extremely dense cluster.
-   A global approach may not consider that point as an outlier. But LOF can effectively identify local outliers. It gives better results than the global approach to find outliers

*Cons*

-   Since LOF is a ratio, it is difficult to interpret
-   There is no threshold value to LOF and the selection of a point as an anomaly/outlier is user-dependent (and also problem-dependent)

Let's start with the earthquake data:

```{r}
set.seed(665544)

### calculate LOF score with a neighborhood of 4 points
lof_results <- lof(df, minPts = 4) # we compute the LOF score with a number of neighbours of 4 (the same number that we chose for DBSCAN clustering!)
```

Let's check our distribution of Local Outlier Factors!

```{r}
### distribution of outlier factors
summary(lof_results)
hist(lof_results, breaks = 10, main = "LOF (minPts = 4)")
```

```{r}
### plot sorted lof. Looks like outliers start around a LOF of 2.4.
plot(sort(lof_results), type = "l",  main = "LOF (minPts = 4)",
  xlab = "Points sorted by LOF", ylab = "LOF")
abline(h=2.4,lty=2)
```

And finally...

```{r}
anomaly_data <- df %>% add_column(score = lof_results)
ggplot(anomaly_data, aes(x = depth, y = mag, color = score, size = score)) +
  geom_point(alpha = 0.5) +  # Add transparency for better visualization
  scale_color_viridis(option = "viridis", direction = 1) +  # Viridis color scale
  scale_size_continuous(range = c(1, 10)) +  # Size range for LOF scores
  theme_minimal() +
  labs(
    title = "Scatter Plot with LOF Scores (Viridis)",
    x = "X Coordinate",
    y = "Y Coordinate",
    color = "LOF Score",
    size = "LOF Score"
  )
```

For more details on LOF, see Breunig, Markus M., Hans-Peter Kriegel, Raymond T. Ng, and Jörg Sander. "LOF: identifying density-based local outliers." In Proceedings of the 2000 ACM SIGMOD international conference on Management of data, pp. 93-104. 2000.

:house_with_garden: **TAKE-HOME EXERCISE**: Can you do the same with the Heart Disease indicators data?

For a more involved example of how to apply LOF (with PCA included!) in a fraud detection context, see [https://rpubs.com/muh_risman/LOF](here)
