4-ClusterAnalysis.R

#' 
#' #### Example exercise: Airports
#' 
#' **Your task**: Create and evaluate the many types of clustering methods. 
#' 
#' #### Variables: 
#' 
#' * `Code`: Code of the airport;  
#' * `Airport`: Name of the airport;  
#' * `Ordem`: ID of the observations;  
#' * `Passengers`: Number of passengers;  
#' * `Movements`: Number of flights;  
#' * `Numberofairlines`: Number of airlines in each airport;  
#' * `Mainairlineflightspercentage`: Percentage of flights of the main airline of each airport;  
#' * `Maximumpercentageoftrafficpercountry`: Maximum percentage of flights per country;  
#' * `NumberofLCCflightsweekly`: Number of weekly low cost flights`;   
#' * `NumberofLowCostAirlines`: Number of low cost airlines of each airport;  
#' * `LowCostAirlinespercentage`: Percentage of the number of low cost airlines in each airport;  
#' * `Destinations`: Number of flights arriving at each airport;  
#' * `Average_route_Distance`: Average route distance in km;  
#' * `DistancetoclosestAirport`: Distance to closest airport in km
#' * `DistancetoclosestSimilarAirport`: Distance to closest similar airport in km;  
#' * `AirportRegionalRelevance`: Relevance of the airport in a regional scale (0 - 1);  
#' * `Distancetocitykm`: Distance between the airport and the city in km;  
#' * `Inhabitantscorrected`: Population of the city;  
#' * `numberofvisitorscorrected`: Number of vistors that arrived in the airport;  
#' * `GDP corrected`: Corrected value of the Gross Domestic Product;  
#' * `Cargoton`: Cargo ton. The total number of cargo transported in a certain period multiplied by the number o flights.
#' 
#' ## Startup
#' ##### Import Libraries
library(readxl) # Reading excel files
library(skimr) # Summary statistics
library(tidyverse) # Pack of useful tools
library(mclust) # Model based clustering
library(cluster) # Cluster analysis
library(factoextra) # Visualizing distances

#' 
#' ##### Import dataset as a dataframe
dataset <- read_excel("Data/Data_Aeroports_Clustersv1.xlsX")
df <- data.frame(dataset)

#' 
#' ## Get to know your data
#' ##### Summary statistics
skim(df)

#' ##### Sneak the plot
#' Now let us plot an example and take a look
plot(Numberofairlines ~ Destinations, df, xlim = c(min(df$Destinations), 1.4*max(df$Destinations)), ylim = c(min(df$Numberofairlines), 1.2*max(df$Numberofairlines))) #plot
text(Numberofairlines ~ Destinations, df, label = Airport, pos = 4, cex = 0.6) #labels over the previous plot

#' 
#' By looking at the plot, you may already have a clue on the number of clusters with this two variabls. However, this is not clear and it does not consider the other variables in the analysis. 
#' 
#' ### Prepare data before performing a cluster analysis
#' ##### Deal with missing data
table(is.na(df))

#' In this example we do not have missing values. In case you do have in the future, you can take out the missing values with _listwise deletion_ (`df <- na.omit(df)`) or use other ways of treating missing values.
#' 
#' ##### Continuous variables
#' Leave only continuous variables, and make `Ordem` as the row ID variable
df_reduced = df[,!(names(df) %in% c("Code","Airport"))]
df_reduced = data.frame(df_reduced, row.names = 1) #Ordem is the 1st variable in the df

#' 
#' Take a look at the scale of the variables. See how they are different!
head(df_reduced)

#' 
#' ##### Standardize variables
#' Z-score standardization:  xi - xmean / sd
mean <- apply(df_reduced, 2, mean) # The "2" in the function is used to select the columns. MARGIN: c(1,2)
sd <- apply(df_reduced, 2, sd)
df_scaled <- scale(df_reduced, mean, sd)

#'   
#' ## Hierarchical Clustering
#' 
#' ##### Measuring Similarity through Euclidean distances
#' 
distance <- dist(df_scaled, method = "euclidean")

#' 
#' > **Note:** There are other forms of distance measures that can be used such as:  
#'   i) Minkowski distance; 
#'   ii) Manhattan distance; 
#'   iii) Mahanalobis distance.
#'   
#' ##### Visualize distances in heatmap
#'
fviz_dist(distance, gradient = list(low = "#00AFBB", mid = "white", high = "#FC4E07"), order = FALSE)

#' 
#' #### Types of hierarchical clustering  
#' There are many types of hierarchical clustering. Let's explore some.
#' 
#' **1. Single linkage (nearest neighbor) clustering algorithm**
#' 
#' Based on a bottom-up approach, by linking two clusters that have the closest distance between each other. 
models <- hclust(distance, "single")
plot(models, labels = df$Airport, xlab = "Distance - Single linkage", cex=0.6, hang = -1)

rect.hclust(models, 4, border = "purple") # Visualize the cut on the dendogram, with 4 clusters

#' 
#' 
#' **2. Complete linkage (Farthest neighbor) clustering algorithm**
#' 
#' Complete linkage is based on the maximum distance between observations in each cluster.
modelc <- hclust(distance, "complete")
plot(modelc, labels = df$Airport, xlab = "Distance - Complete linkage", cex=0.6, hang = -1)
rect.hclust(modelc, 4, border = "blue") 

#' 
#' **3. Average linkage between groups** 
#' 
#' The average linkage considers the distance between clusters to be the average of the distances between observations in one cluster to all the members in the other cluster. 
modela <- hclust(distance, "average")
plot(modela, labels = df$Airport, xlab = "Distance - Average linkage", cex=0.6, hang = -1)
rect.hclust(modela, 4, border = "red")

#' 
#' **4. Ward`s method**
#'   
#' The Ward`s method considers the measures of similarity as the sum of squares within the cluster summed over all variables. 
modelw <- hclust(distance, "ward.D2")
plot(modelw, labels = df$Airport, xlab = "Distance - Ward method", cex=0.6, hang = -1)
rect.hclust(modelw, 4, border = "orange")

#' 
#' **5. Centroid method**
#' 
#' The centroid method considers the similarity between two clusters as the distance between its centroids.
modelcen <- hclust(distance, "centroid")
plot(modelcen, labels = df$Airport, xlab = "Distance - Centroid method", cex=0.6, hang = -1)
rect.hclust(modelcen, 4, border = "darkgreen")

#' 
#' ##### Comparing results from different hierarchical methods
#' Now lets evaluate the **membership** of each observation with the `cutree` function for each method.
#' 
member_single <- cutree(models, 4)
member_com <- cutree(modelc, 4)
member_av <- cutree(modela, 4)
member_ward <- cutree(modelw, 4)
member_cen <- cutree(modelcen, 4)

#'   
#' 
#' Compare how common each method is to each other.
table(member_com, member_av) # compare the complete linkage with the average linkage

#' 
#' > **Note:** Try comparing other methods, and evaluate how common they are.   
#'   
#' ##### Silhouette Plots: evaluate which method is more appropriate 
#' The silhouette plot evaluates how similar an observation is to its own cluster compared to other clusters. The clustering configuration is appropriate when most objects have high values. Low or negative values indicate that the clustering method is not appropriate or the number of clusters is not ideal. 
#' 
plot(silhouette(member_single, distance))
plot(silhouette(member_com, distance))
plot(silhouette(member_av, distance))
plot(silhouette(member_ward, distance))
plot(silhouette(member_cen, distance))

#' 
#' 
#' ## Non-Hirarchical Clustering 
#' 
#' ##### K-means clustering
#' * k-means with n=3 clusters
km_clust <- kmeans(df_scaled, 3)
km_clust #print the results

#' 
#' * Other ways of setting the number of clusters  
#' 
#' This algorithm will detect how many clusters from 1 to 10 explains more variance
#' 
  k <- list()
  for(i in 1:10){
    k[[i]] <- kmeans(df_scaled, i)
  }

#'   
#' > **Note**: Try printing the k value and take a look at the ratio `between_SS` / `total_SS`. Evaluate how it varies when you add clusters. 
#' 
#' Now, lets  plot `between_SS` / `total_SS` into a scree plot
#' 
betSS_totSS <- list()
for(i in 1:10){
betSS_totSS[[i]] <- k[[i]]$betweenss/k[[i]]$totss
}
plot(1:10, betSS_totSS, type = "b", ylab = "Between SS / Total SS", xlab = "Number of clusters")

#' 
#' 
#' Let take out the outliers and see the difference in the k-means clustering:
#'  
#' * **Examine the boxplots**
#'   
par(cex.axis=0.6, mar=c(11,2,1,1))# Make labels fit in the boxplot
boxplot(df_scaled, las = 2) #labels rotated to vertical

#' 
#' * **Detect the outliers** 
#' 
outliers <- boxplot.stats(df_scaled)$out
outliers

#' 
#' * **Remove the biggest outlier** 
#' 
nrow(df_scaled) #32
out_ind <- which(df_scaled %in% c(outliers)) #the row.names that contain outliers
df_big_outlier = df_scaled[-c(out_ind),] #remove those rows from the df_scaled
nrow(df_big_outlier) #31

#' > **Note:** There are many methods to treat outliers. This is just one of them. Note that it is not very appropriate, since it removes many observations that are relevant for the analysis. Try using other methods and evaluate the difference.  
#' 
#' 
#' Execute a k-means clustering with the dataset without the outliers and see the difference. 
#' 
km_big_outlier <- kmeans(df_big_outlier, 3)
km_big_outlier

#' 
#' #### Ploting the clusters
#' Finally, plotting the clusters results to check if they make sense.  
#' Let us go back to first example and take a look. 
#' 
#' * **K-means with outliers**
#' 
plot(Numberofairlines ~ Destinations, df, col = km_clust$cluster, xlim = c(1.1*min(df$Destinations), 1.2*max(df$Destinations)), ylim = c(min(df$Numberofairlines), 1.2*max(df$Numberofairlines)))
with(df, text(Numberofairlines ~ Destinations, label = Airport, pos = 1, cex = 0.4))

#' 
#' * **K-means without outliers**
#'
plot(Numberofairlines ~ Destinations, df, col = km_big_outlier$cluster, xlim = c(1.1*min(df$Destinations), 1.4*max(df$Destinations)), ylim = c(1.2*min(df$Numberofairlines), 1.2*max(df$Numberofairlines)))
with(df, text(Numberofairlines ~ Destinations, label = Airport, pos = 1, cex = 0.6))