London Velib Notebook

##Clean DATA##

setwd("C:/Users/HSBC/OneDrive/桌面")

velib <- read.table("london_merged.txt", header = TRUE, sep =",",quote = "")

# Lignes restantes
n <- nrow(velib) - 72

# Les supprimes
velib <- velib[1:n, ]

# Suppr les guillemets des colonnes "timestamp" et "season".
velib$X.timestamp <- gsub('"', '', velib$X.timestamp)
velib$season. <- gsub('"', '', velib$season.)

# make sure it is in the date format
velib$X.timestamp <- as.POSIXct(velib$X.timestamp, format = "%Y-%m-%d %H:%M:%S")
# hours
velib$heure2 <- as.numeric(format(velib$X.timestamp, "%H"))
velib$season. <- as.numeric(velib$season.)
summary(velib)

  X.timestamp                          cnt             t1             t2       
 Min.   :2015-01-04 00:00:00.00   Min.   :   0   Min.   :-1.5   Min.   :-6.00  
 1st Qu.:2015-07-04 02:15:00.00   1st Qu.: 260   1st Qu.: 8.5   1st Qu.: 6.00  
 Median :2016-01-02 03:30:00.00   Median : 846   Median :12.5   Median :12.50  
 Mean   :2016-01-02 09:36:27.56   Mean   :1146   Mean   :12.5   Mean   :11.56  
 3rd Qu.:2016-07-02 09:45:00.00   3rd Qu.:1676   3rd Qu.:16.0   3rd Qu.:16.00  
 Max.   :2016-12-31 23:00:00.00   Max.   :7860   Max.   :34.0   Max.   :34.00  
      hum           wind_speed     weather_code      is_holiday     
 Min.   : 20.50   Min.   : 0.00   Min.   : 1.000   Min.   :0.00000  
 1st Qu.: 63.00   1st Qu.:10.00   1st Qu.: 1.000   1st Qu.:0.00000  
 Median : 74.50   Median :15.00   Median : 2.000   Median :0.00000  
 Mean   : 72.28   Mean   :15.92   Mean   : 2.721   Mean   :0.02076  
 3rd Qu.: 83.00   3rd Qu.:20.50   3rd Qu.: 3.000   3rd Qu.:0.00000  
 Max.   :100.00   Max.   :56.50   Max.   :26.000   Max.   :1.00000  
   is_weekend        season.          heure2     
 Min.   :0.0000   Min.   :0.000   Min.   : 0.00  
 1st Qu.:0.0000   1st Qu.:0.000   1st Qu.: 6.00  
 Median :0.0000   Median :1.000   Median :12.00  
 Mean   :0.2852   Mean   :1.486   Mean   :11.51  
 3rd Qu.:1.0000   3rd Qu.:2.000   3rd Qu.:18.00  
 Max.   :1.0000   Max.   :3.000   Max.   :23.00  

##ACP##

library(FactoMineR)

Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

velib2 <- velib[-1] #suppr timestamp
velib3 <- velib2[-3] #suppr t2
  
vel.pca <- PCA(velib2, graph = TRUE)

vel.pca <- PCA(velib2, quali.sup=6:10,graph=T)

pc_velib_norm <- prcomp(velib2, scale. = TRUE)
summary(pc_velib_norm)

Importance of components:
                          PC1    PC2    PC3    PC4     PC5     PC6     PC7
Standard deviation     1.7101 1.1078 1.0816 1.0426 0.97663 0.93496 0.88172
Proportion of Variance 0.2924 0.1227 0.1170 0.1087 0.09538 0.08741 0.07774
Cumulative Proportion  0.2924 0.4152 0.5322 0.6409 0.73623 0.82365 0.90139
                           PC8     PC9    PC10
Standard deviation     0.77013 0.61940 0.09674
Proportion of Variance 0.05931 0.03837 0.00094
Cumulative Proportion  0.96070 0.99906 1.00000

biplot(pc_velib_norm,col="navy", xlabs = rep("",17342),cex = 0.5)

##CAH##

library(ggplot2)
library(cluster)  # but


# normee, and remove time
velib_clust <- velib %>% select(cnt, t1, hum, wind_speed, is_holiday, is_weekend, season.)
velib_scaled <- scale(velib_clust)
# distance
d <- dist(velib_scaled, method = "euclidean")

if (!requireNamespace("dendextend", quietly = TRUE)) {
  install.packages("dendextend")
}
library(dendextend)

hc <- hclust(d, method = "complete")
# Convertir en objet dendrogramme
dend <- as.dendrogram(hc)

dend <- dend %>% 
  set("labels_cex", 0.7) %>%  # la taille de l'étiquette
  set("nodes_pch", 19) %>%    # la forme du nœud
  set("nodes_cex", 0.5)       # la size du nœud


par(mar = c(5, 5, 2, 2))  # marges graphiques
plot(dend, horiz = TRUE, main = "Cluster Dendrogram")

# Coefficients de contour pour différents nombres de grappes
sil_widths <- lapply(2:10, function(k) {
  clusters <- cutree(hc, k)
  sil <- silhouette(clusters, d)
  
  if (dim(sil)[1] == 0) {
    return(NA)  # # Si sil n'a pas de lignes, retourner NA
  } else {
    return(mean(sil[, 3]))  #  la largeur moyenne du profil
  }
})

# suppri NA
valid_sil_widths <- sil_widths[!is.na(sil_widths)]

if (length(valid_sil_widths) > 0) {
  best_k <- which.max(valid_sil_widths)
  
  plot(2:10, sil_widths, type = 'b', xlab = "Number of clusters", ylab = "Average Silhouette Width",
       main = "Silhouette analysis for choosing number of clusters")
  abline(v = best_k + 1, col = "red")  #le nmbr optimal de grappes.
} else {
  print("No valid silhouette widths were calculated. Please check data or distance measure.")
}

library(dplyr)

stats_cluster1 <- summary(velib[velib$cluster == 1, ])
stats_cluster2 <- summary(velib[velib$cluster == 2, ])

print("Descriptive statistics for Cluster 1:")

[1] "Descriptive statistics for Cluster 1:"

print(stats_cluster1)

  X.timestamp                         cnt             t1              t2       
 Min.   :2015-01-04 00:00:00.0   Min.   :   0   Min.   :-1.50   Min.   :-6.00  
 1st Qu.:2015-07-04 08:15:00.0   1st Qu.: 263   1st Qu.: 8.50   1st Qu.: 6.00  
 Median :2016-01-02 15:30:00.0   Median : 856   Median :12.50   Median :12.50  
 Mean   :2016-01-02 01:06:02.5   Mean   :1153   Mean   :12.52   Mean   :11.59  
 3rd Qu.:2016-07-03 03:45:00.0   3rd Qu.:1689   3rd Qu.:16.50   3rd Qu.:16.50  
 Max.   :2016-12-31 23:00:00.0   Max.   :7860   Max.   :34.00   Max.   :34.00  
      hum           wind_speed     weather_code      is_holiday   is_weekend    
 Min.   : 20.50   Min.   : 0.00   Min.   : 1.000   Min.   :0    Min.   :0.0000  
 1st Qu.: 63.00   1st Qu.:10.00   1st Qu.: 1.000   1st Qu.:0    1st Qu.:0.0000  
 Median : 74.50   Median :15.00   Median : 2.000   Median :0    Median :0.0000  
 Mean   : 72.22   Mean   :15.91   Mean   : 2.714   Mean   :0    Mean   :0.2912  
 3rd Qu.: 83.00   3rd Qu.:20.50   3rd Qu.: 3.000   3rd Qu.:0    3rd Qu.:1.0000  
 Max.   :100.00   Max.   :56.50   Max.   :26.000   Max.   :0    Max.   :1.0000  
    season.          heure2      cluster  
 Min.   :0.000   Min.   : 0.00   1:16982  
 1st Qu.:1.000   1st Qu.: 6.00   2:    0  
 Median :1.000   Median :12.00            
 Mean   :1.493   Mean   :11.51            
 3rd Qu.:2.000   3rd Qu.:18.00            
 Max.   :3.000   Max.   :23.00            

print("Descriptive statistics for Cluster 2:")

[1] "Descriptive statistics for Cluster 2:"

print(stats_cluster2)

  X.timestamp                       cnt               t1              t2       
 Min.   :2015-04-03 00:00:00   Min.   :  14.0   Min.   : 2.00   Min.   :-2.00  
 1st Qu.:2015-05-25 17:45:00   1st Qu.: 179.2   1st Qu.: 8.00   1st Qu.: 5.50  
 Median :2016-01-01 11:30:00   Median : 450.0   Median :11.50   Median :11.50  
 Mean   :2016-01-19 02:54:00   Mean   : 788.0   Mean   :11.38   Mean   :10.34  
 3rd Qu.:2016-05-30 05:15:00   3rd Qu.:1087.5   3rd Qu.:15.00   3rd Qu.:15.00  
 Max.   :2016-12-27 23:00:00   Max.   :3100.0   Max.   :23.00   Max.   :23.00  
      hum           wind_speed     weather_code      is_holiday   is_weekend
 Min.   : 40.50   Min.   : 1.00   Min.   : 1.000   Min.   :1    Min.   :0   
 1st Qu.: 66.00   1st Qu.:10.00   1st Qu.: 1.000   1st Qu.:1    1st Qu.:0   
 Median : 77.00   Median :15.00   Median : 2.000   Median :1    Median :0   
 Mean   : 74.86   Mean   :16.08   Mean   : 3.039   Mean   :1    Mean   :0   
 3rd Qu.: 85.00   3rd Qu.:20.00   3rd Qu.: 4.000   3rd Qu.:1    3rd Qu.:0   
 Max.   :100.00   Max.   :54.00   Max.   :26.000   Max.   :1    Max.   :0   
    season.          heure2      cluster
 Min.   :0.000   Min.   : 0.00   1:  0  
 1st Qu.:0.000   1st Qu.: 5.75   2:360  
 Median :0.000   Median :11.50          
 Mean   :1.133   Mean   :11.50          
 3rd Qu.:3.000   3rd Qu.:17.25          
 Max.   :3.000   Max.   :23.00          

# Température et cnt par groupe
library(ggplot2)
ggplot(velib, aes(x = t1, y = cnt, color = factor(cluster))) + 
  geom_point(alpha = 0.5) +
  facet_wrap(~cluster) +  # afficher chaque groupe de manière facettée
  labs(title = "Temperature vs Bike Count by Cluster", x = "Temperature (°C)", y = "Bike Count") +
  scale_color_manual(values = c("red", "blue"), labels = c("Cluster 1", "Cluster 2")) +
  theme_minimal()

clusters <- cutree(hc, k = 2)  
velib$cluster <- as.factor(clusters)  

summary_stats <- aggregate(velib[, c("cnt", "t1", "hum", "wind_speed")], 
                           by = list(velib$cluster), 
                           FUN = function(x) c(mean = mean(x), sd = sd(x)))

print(summary_stats)

p <- ggplot(velib, aes(x = t1, y = cnt, color = cluster)) +
  geom_point(alpha = 0.5) +
  labs(title = "Temperature vs Bike Count by Cluster",
       x = "Temperature (°C)",
       y = "Bike Count") +
  scale_color_manual(values = c("red", "blue"), labels = c("Cluster 1", "Cluster 2")) +
  theme_minimal()
print(p)

p2 <- ggplot(velib, aes(x = wind_speed, y = hum, color = cluster)) +
  geom_point(alpha = 0.5) +
  labs(title = "Wind Speed vs Humidity by Cluster",
       x = "Wind Speed (km/h)",
       y = "Humidity (%)") +
  scale_color_manual(values = c("red", "blue"), labels = c("Cluster 1", "Cluster 2")) +
  theme_minimal()
print(p2)

# Matrice de corrélation pour chaque groupe
cor_cluster1 <- cor(velib[velib$cluster == '1', c("t1", "hum", "wind_speed", "cnt")])
cor_cluster2 <- cor(velib[velib$cluster == '2', c("t1", "hum", "wind_speed", "cnt")])

# print corr
print("Correlation Matrix for Cluster 1:")

[1] "Correlation Matrix for Cluster 1:"

print(cor_cluster1)

                   t1        hum wind_speed        cnt
t1          1.0000000 -0.4458184  0.1476201  0.3848495
hum        -0.4458184  1.0000000 -0.2903531 -0.4584796
wind_speed  0.1476201 -0.2903531  1.0000000  0.1179926
cnt         0.3848495 -0.4584796  0.1179926  1.0000000

print("Correlation Matrix for Cluster 2:")

[1] "Correlation Matrix for Cluster 2:"

print(cor_cluster2)

                    t1        hum   wind_speed          cnt
t1          1.00000000 -0.4106834  0.024545124  0.556330336
hum        -0.41068340  1.0000000 -0.168119965 -0.659225420
wind_speed  0.02454512 -0.1681200  1.000000000  0.008825807
cnt         0.55633034 -0.6592254  0.008825807  1.000000000

# 
library(corrplot)
corrplot(cor_cluster1, method = "circle", title = "Correlation Matrix for Cluster 1")

corrplot(cor_cluster2, method = "circle", title = "Correlation Matrix for Cluster 2")