1. Load the data [uber-data.csv]

uber<-fread("uber-data.csv",showProgress=FALSE)
#show the first 6 observations
knitr::kable(head(uber))
Date/Time Lat Lon Base
9/1/2014 0:01 40.2201 -74.0021 B02512
9/1/2014 0:01 40.7500 -74.0027 B02512
9/1/2014 0:03 40.7559 -73.9864 B02512
9/1/2014 0:06 40.7450 -73.9889 B02512
9/1/2014 0:11 40.8145 -73.9444 B02512
9/1/2014 0:12 40.6735 -73.9918 B02512 
#check the structure 
str(uber)
## Classes 'data.table' and 'data.frame':   1028136 obs. of  4 variables:
##  $ Date/Time: chr  "9/1/2014 0:01" "9/1/2014 0:01" "9/1/2014 0:03" "9/1/2014 0:06" ...
##  $ Lat      : num  40.2 40.8 40.8 40.7 40.8 ...
##  $ Lon      : num  -74 -74 -74 -74 -73.9 ...
##  $ Base     : chr  "B02512" "B02512" "B02512" "B02512" ...
##  - attr(*, ".internal.selfref")=<pointer: 0x000001a71b251230>
#summarize the data quality
initial_summary <- data.table(
  observations = nrow(uber),
  variables = ncol(uber),
  missing_values = sum(is.na(uber)),
  duplicate_rows = sum(duplicated(uber)),
  unique_bases = uniqueN(uber$Base)
)

knitr::kable(initial_summary)
observations variables missing_values duplicate_rows unique_bases
1028136 4 0 24037 5 
print(initial_summary)
##    observations variables missing_values duplicate_rows unique_bases
##           <int>     <int>          <int>          <int>        <int>
## 1:      1028136         4              0          24037            5

The “uber-data.csv” contains 1,028,136 observations and 4 variables (no missing values). With 24,037 duplicate records and 5 unique Uber base codes. I kept the duplicates because identical records may represent separate trips.

1.1. Check datetime and geographic coordinates

I am going to convert the date&time variable into a proper datetime format and check whether all timestamps are valid. Also check the latitude and longitude ranges to identify invalid geographical coordinates before spatial clustering.

#convert date and time column
uber[,datetime:=mdy_hm(`Date/Time`,tz="America/New_York")] #to use the NY local time zone 

#check datetime validity & coverage
datetime_check<-uber[,.(
  invalid_datetime = sum(is.na(datetime)),
  first_timestamp = min(datetime, na.rm = TRUE),
  last_timestamp = max(datetime, na.rm = TRUE)
)]

knitr::kable(datetime_check)
invalid_datetime first_timestamp last_timestamp
0 2014-09-01 2014-09-30 22:59:00 
#check if geographic coordinate is valid
coordinate_check <- uber[, .(
  invalid_coordinates = sum(
    !is.finite(Lat) | !is.finite(Lon) |
    Lat < -90 | Lat > 90 |
    Lon < -180 | Lon > 180
  ),
  min_latitude = min(Lat, na.rm = TRUE),
  max_latitude = max(Lat, na.rm = TRUE),
  min_longitude = min(Lon, na.rm = TRUE),
  max_longitude = max(Lon, na.rm = TRUE)
)]

knitr::kable(coordinate_check)
invalid_coordinates min_latitude max_latitude min_longitude max_longitude
0 39.9897 41.3476 -74.7736 -72.7163 
print(datetime_check)
##    invalid_datetime first_timestamp      last_timestamp
##               <int>          <POSc>              <POSc>
## 1:                0      2014-09-01 2014-09-30 22:59:00
print(coordinate_check)
##    invalid_coordinates min_latitude max_latitude min_longitude max_longitude
##                  <int>        <num>        <num>         <num>         <num>
## 1:                   0      39.9897      41.3476      -74.7736      -72.7163

All timestamps were converted into datetime format, with 0 invalid values. The geographic coordinates are also valid, with latitude from 39.9897 to 41.3476 and longitude from -74.7736 to -72.7163.

1.2. Create time variables and check temporal coverage

This step, we extract the date, hour, and day of the week from the datetime variable to support temporal clustering analysis and forecasting. Then, we also check the observation period and hourly coverage to see incomplete days.

#create temporal variables
uber[,`:=`(
  date = as.Date(datetime,tz="America/New_York"),
  hour = lubridate::hour(datetime),
  day_of_week = lubridate::wday(datetime, week_start = 1)
)]

#check first observations
knitr::kable(head(uber[, .(datetime, date, hour, day_of_week)]))
datetime date hour day_of_week
2014-09-01 00:01:00 2014-09-01 0 1
2014-09-01 00:01:00 2014-09-01 0 1
2014-09-01 00:03:00 2014-09-01 0 1
2014-09-01 00:06:00 2014-09-01 0 1
2014-09-01 00:11:00 2014-09-01 0 1
2014-09-01 00:12:00 2014-09-01 0 1 
#summarize the observation period
time_summary<-uber[, .(
  first_timestamp = min(datetime),
  last_timestamp = max(datetime),
  total_dates = uniqueN(date),
  total_hours = uniqueN(hour)
)]

print(time_summary)
##    first_timestamp      last_timestamp total_dates total_hours
##             <POSc>              <POSc>       <int>       <int>
## 1:      2014-09-01 2014-09-30 22:59:00          30          24
#check hourly coverage for each date
hour_coverage <- uber[, .(
  trips = .N,
  observed_hours = uniqueN(hour)
), by = date][order(date)]

#display dates with incomplete hourly coverage
print(hour_coverage[observed_hours < 24])
##          date trips observed_hours
##        <Date> <int>          <int>
## 1: 2014-09-30 33431             23
#check original and converted datetime
print(head(uber[, .(original = `Date/Time`, parsed = datetime)]))
##         original              parsed
##           <char>              <POSc>
## 1: 9/1/2014 0:01 2014-09-01 00:01:00
## 2: 9/1/2014 0:01 2014-09-01 00:01:00
## 3: 9/1/2014 0:03 2014-09-01 00:03:00
## 4: 9/1/2014 0:06 2014-09-01 00:06:00
## 5: 9/1/2014 0:11 2014-09-01 00:11:00
## 6: 9/1/2014 0:12 2014-09-01 00:12:00
#check datetime validity
print(sum(is.na(uber$datetime)))
## [1] 0

As we can see, the dataset covers 30 days in September 2014, with observations across all 24hrs of the day. September 30 contains only 23 observed hours (33,431 records).

2. Spatial clustering

2.1. Select algorithm and optimal clusters

We use K-means to identify geographical groups of uber trips based on lat and lon. I believe K-means is suitable for this analysis because it provides cluster centers that can be used in the subsequent temporal analysis and forecasting.

To specify number of clusters, we compare different values of k using the Elbow Method, Silhouette Score, and Calinski–Harabasz Index. Since the dataset is large, model selection is performed on a reproducible random sample of training observations.

#select training observations
train_rows<-which(uber$date < as.Date("2014-09-24"))

#create a reproducible sample for model selection
set.seed(42)
sample_rows <- sample(train_rows, size = min(20000, length(train_rows)))
sample_data <- uber[sample_rows, .(Lon, Lat)]

#convert coordinates to a projected system in meters
sample_sf <- st_as_sf(
  sample_data,
  coords = c("Lon", "Lat"),
  crs = 4326
)

sample_sf <- st_transform(sample_sf, crs = 32618)
X <- st_coordinates(sample_sf)

#select observations for Silhouette evaluation
set.seed(42)
sil_rows <- sample(seq_len(nrow(X)), size = 1000)
sil_dist <- dist(X[sil_rows, , drop = FALSE])

#test different numbers of clusters
k_values <- 2:12

k_results <- data.table::rbindlist(lapply(k_values, function(k) {

  #fit K-means for the current number of clusters
  model <- kmeans(X, centers = k, nstart = 10, iter.max = 500, algorithm = "Lloyd")

  #calculate average Silhouette Score
  sil <- cluster::silhouette(model$cluster[sil_rows], sil_dist)
  avg_sil <- mean(sil[, "sil_width"])

  #calculate Calinski-Harabasz Index
  ch <- (model$betweenss / (k - 1)) /
        (model$tot.withinss / (nrow(X) - k))

  data.table::data.table(
    k = k,
    inertia = model$tot.withinss,
    silhouette = avg_sil,
    calinski_harabasz = ch
  )
}))

#display quality measures
print(round(k_results, 3))
##         k      inertia silhouette calinski_harabasz
##     <num>        <num>      <num>             <num>
##  1:     2 626060409323      0.652          8664.566
##  2:     3 455702670263      0.387          9689.338
##  3:     4 369815930522      0.413          9507.310
##  4:     5 305865314267      0.439          9666.043
##  5:     6 244766108085      0.433         10660.835
##  6:     7 201839938490      0.435         11481.561
##  7:     8 170296582423      0.449         12192.632
##  8:     9 149005208970      0.407         12549.438
##  9:    10 136125041255      0.434         12420.100
## 10:    11 113977206844      0.434         13737.950
## 11:    12 104958181362      0.437         13717.689
#plot Elbow Method
ggplot(k_results, aes(x = k, y = inertia)) +
  geom_line() +
  geom_point() +
  scale_x_continuous(breaks = k_values) +
  labs(title = "Elbow Method", x = "Number of clusters",
       y = "Within-Cluster Sum of Squares") +
  theme_minimal()

#plot Silhouette Scores
ggplot(k_results, aes(x = k, y = silhouette)) +
  geom_line() +
  geom_point() +
  scale_x_continuous(breaks = k_values) +
  labs(title = "Silhouette Score by Number ofclusters",
       x = "Number of Clusters", y = "Average Silhouette") +
  theme_minimal()

#compare cluster sizes for k=5 & k=8
set.seed(42)
model_5 <- kmeans(X, centers = 5, nstart = 10,iter.max = 500, algorithm = "Lloyd")

model_8 <- kmeans(X, centers = 8, nstart = 10,iter.max = 500, algorithm = "Lloyd")

#calculate cluster sizes,  percentages and display
print(size_5 <- data.table(
  cluster = 1:5,
  count = model_5$size,
  percentage = round(100 * model_5$size / nrow(X), 2)
))
##    cluster count percentage
##      <int> <int>      <num>
## 1:       1   660       3.30
## 2:       2  1368       6.84
## 3:       3   215       1.07
## 4:       4  7463      37.31
## 5:       5 10294      51.47
print(size_8 <- data.table(
  cluster = 1:8,
  count = model_8$size,
  percentage = round(100 * model_8$size / nrow(X), 2)
))
##    cluster count percentage
##      <int> <int>      <num>
## 1:       1   528       2.64
## 2:       2   215       1.07
## 3:       3   472       2.36
## 4:       4  2088      10.44
## 5:       5  6958      34.79
## 6:       6    73       0.36
## 7:       7   899       4.50
## 8:       8  8767      43.84

The Elbow Method doesn’t show a clear optimal number of clusters. Although k=2 gets the highest Silhouette score =0.652 (limited geographical detail).

Among the larger solutions, k=8 gets a slightly higher Silhouette Score =0.449 than k=5 0.439. However, k=8 also produces a very small cluster (only 0.36% of the sampled observations).

Both solutions show unequal cluster sizes. Therefore, for now we will keep k=5 and k=8 for further geographical evaluation.

2.2. Check clustering results and quality

We compare the clustering results for k=5 and k=8 using the Silhouette Score, Calinski–Harabasz Index, and within-cluster sum of squares.

#compare clustering quality
models<-list(model_5, model_8)

quality_results <- rbindlist(lapply(seq_along(models), function(i) {
  model <- models[[i]]
  k <- c(5, 8)[i]

  #Silhouette Score
  sil <- silhouette(model$cluster[sil_rows], sil_dist)

  #Calinski-Harabasz Index
  ch <- (model$betweenss / (k - 1)) /
    (model$tot.withinss / (nrow(X) - k))

  data.table(
    k = k,
    inertia = round(model$tot.withinss, 2),
    silhouette = round(mean(sil[, "sil_width"]), 3),
    calinski_harabasz = round(ch, 2)
  )
}))

print(quality_results)
##        k      inertia silhouette calinski_harabasz
##    <num>        <num>      <num>             <num>
## 1:     5 300313928341      0.429           9937.13
## 2:     8 170296582423      0.449          12192.63
#extract cluster centers and sizes
centers <- rbindlist(list(
  data.table(
    solution = "k = 5",
    cluster = 1:5,
    x = model_5$centers[, 1],
    y = model_5$centers[, 2],
    count = model_5$size
  ),
  data.table(
    solution = "k = 8",
    cluster = 1:8,
    x = model_8$centers[, 1],
    y = model_8$centers[, 2],
    count = model_8$size
  )
))

#convert cluster centers back to longitude and latitude
centers_sf <- st_as_sf(
  centers, coords = c("x", "y"), crs = 32618
)

centers_sf <- st_transform(centers_sf, 4326)
center_coords <- st_coordinates(centers_sf)

centers[, `:=`(
  Lon = center_coords[, 1],
  Lat = center_coords[, 2],
  percentage = round(100 * count / nrow(X), 2)
)]

#display cluster centers
print(centers[, .(
  solution, cluster, Lon, Lat, count, percentage
)])
##     solution cluster       Lon      Lat count percentage
##       <char>   <int>     <num>    <num> <int>      <num>
##  1:    k = 5       1 -73.76770 40.67029   660       3.30
##  2:    k = 5       2 -73.88954 40.80268  1368       6.84
##  3:    k = 5       3 -74.20406 40.69833   215       1.07
##  4:    k = 5       4 -73.98668 40.70974  7463      37.31
##  5:    k = 5       5 -73.98058 40.75778 10294      51.47
##  6:    k = 8       1 -73.77873 40.65348   528       2.64
##  7:    k = 8       2 -74.20406 40.69833   215       1.07
##  8:    k = 8       3 -73.92343 40.84774   472       2.36
##  9:    k = 8       4 -73.96666 40.67811  2088      10.44
## 10:    k = 8       5 -73.97356 40.76582  6958      34.79
## 11:    k = 8       6 -73.58872 40.90679    73       0.36
## 12:    k = 8       7 -73.86748 40.76017   899       4.50
## 13:    k = 8       8 -73.99452 40.72985  8767      43.84
#sample observations for geographic visualization
set.seed(42)
plot_rows <- sample(seq_len(nrow(X)), 7000)

plot_data <- rbindlist(list(
  data.table(
    Lon = sample_data$Lon[plot_rows],
    Lat = sample_data$Lat[plot_rows],
    cluster = factor(model_5$cluster[plot_rows]),
    solution = "k = 5"
  ),
  data.table(
    Lon = sample_data$Lon[plot_rows],
    Lat = sample_data$Lat[plot_rows],
    cluster = factor(model_8$cluster[plot_rows]),
    solution = "k = 8"
  )
))

#visualize
ggplot(plot_data, aes(x = Lon, y = Lat, color = cluster)) +
  geom_point(alpha = 0.4, size = 0.5) +
  geom_point(
    data = centers,
    aes(x = Lon, y = Lat),
    inherit.aes = FALSE,
    shape = 4,
    color = "black",
    size = 3,
    stroke = 1.2
  ) +
  facet_wrap(~solution) +
  coord_quickmap() +
  labs(
    title = "Geographical Clustering",
    subtitle = "Black crosses indicate cluster centers",
    x = "Longitude",
    y = "Latitude",
    color = "Cluster"
  ) +
  theme_minimal()

k=8 solution gets a higher Silhouette Score 0.449 and Calinski–Harabasz Index 12,192.63 than k=5 (0.429 and 9,937.13,)

The geographical visualization shows k=8 provides a more detailed spatial separation of uber trip locations.

Based on the results, I’ll choose k=8 for further analysis.

2.3. Cluster visualization

Use k=8 clustering model to the dataset using the established cluster centers then examine the number of trips in each cluster and visualize their geographical distribution to check the main areas of uber activity.

#project geographic coordinates
xy<-sf::sf_project("EPSG:4326", "EPSG:32618", as.matrix(uber[, .(Lon, Lat)]))

#assign observations to nearest cluster center
nearest_dist <- rep(Inf, nrow(uber))
uber[, cluster := NA_integer_]

for (i in 1:8) {
  d <- (xy[, 1] - model_8$centers[i, 1])^2 + (xy[, 2] - model_8$centers[i, 2])^2
  nearest <- d < nearest_dist
  uber[nearest, cluster := i]
  nearest_dist[nearest] <- d[nearest]
}

#remove temporary objects
rm(xy, d, nearest, nearest_dist)

#summarize cluster sizes and centers
final_summary <- uber[, .(trips = .N), by = cluster][order(cluster)]
final_summary[, percentage := round(100 * trips / sum(trips), 2)]
final_summary <- merge(final_summary, centers[solution == "k = 8", .(cluster, Lon, Lat)], by = "cluster")

print(final_summary)
## Key: <cluster>
##    cluster  trips percentage       Lon      Lat
##      <int>  <int>      <num>     <num>    <num>
## 1:       1  29147       2.83 -73.77873 40.65348
## 2:       2  10203       0.99 -74.20406 40.69833
## 3:       3  24930       2.42 -73.92343 40.84774
## 4:       4 108247      10.53 -73.96666 40.67811
## 5:       5 359839      35.00 -73.97356 40.76582
## 6:       6   3452       0.34 -73.58872 40.90679
## 7:       7  43268       4.21 -73.86748 40.76017
## 8:       8 449050      43.68 -73.99452 40.72985
knitr::kable(final_summary, digits = 4)
cluster trips percentage Lon Lat
1 29147 2.83 -73.7787 40.6535
2 10203 0.99 -74.2041 40.6983
3 24930 2.42 -73.9234 40.8477
4 108247 10.53 -73.9667 40.6781
5 359839 35.00 -73.9736 40.7658
6 3452 0.34 -73.5887 40.9068
7 43268 4.21 -73.8675 40.7602
8 449050 43.68 -73.9945 40.7298 
#sample observations for visualization
set.seed(42)
plot_final <- uber[sample.int(nrow(uber), 15000), .(Lon, Lat, cluster)]

#visualize final clusters and their centers
ggplot(plot_final, aes(x = Lon, y = Lat, color = factor(cluster))) +
  geom_point(alpha = 0.3, size = 0.5) +
  geom_point(data = final_summary, aes(x = Lon, y = Lat), inherit.aes = FALSE, shape = 4, color = "black", size = 3) +
  geom_text(data = final_summary, aes(x = Lon, y = Lat, label = cluster), inherit.aes = FALSE, nudge_y = 0.025, size = 3) +
  coord_quickmap() +
  labs(title = "Clustering k=8", x = "Longitude", y = "Latitude", color = "Cluster") +
  theme_minimal()

k=8 model was applied to all 1,028,136 uber observations. Clusters 5 & 8 contain the majority of trips 78.68%, showing a strong concentration of uber activity in central New York City. 6 is the smallest 0.34% and represents a more distant geographical area. Overall, the results reveal an uneven spatial distribution of uber trips across the eight clusters.

3. Analyze the cluster centers by time and date

3.1. By time

We analyze uber trip activity and the geographical centers of the 8 clusters across different hrs of the day. For each cluster and hr, we calculate the number of trips and the average geographical location and also measure how far these hourly centers are from the original cluster centers to identify temporal changes in spatial activity.

#convert coordinates to meters
xy<-sf::sf_project("EPSG:4326", "EPSG:32618", as.matrix(uber[, .(Lon, Lat)]))
uber[, `:=`(x = xy[, 1], y = xy[, 2])]
rm(xy)

#hourly cluster centers
hourly_centers <- uber[, .(trips = .N, center_x = mean(x), center_y = mean(y)), by = .(cluster, hour)][order(cluster, hour)]

#distance from the original cluster center
hourly_centers[, shift_km := sqrt((center_x - model_8$centers[cluster, 1])^2 + (center_y - model_8$centers[cluster, 2])^2) / 1000]

#identify the busiest hour for each cluster
peak_hours <- hourly_centers[, .SD[which.max(trips)], by = cluster][, .(cluster, peak_hour = hour, peak_trips = trips)]
print(peak_hours)
##    cluster peak_hour peak_trips
##      <int>     <int>      <int>
## 1:       1        14       2468
## 2:       2        15        845
## 3:       3        21       1663
## 4:       4        21       6628
## 5:       5        17      29238
## 6:       6        17        238
## 7:       7        20       3177
## 8:       8        18      32768
#summarize hourly activity and center shifts
hourly_summary <- hourly_centers[, .(min_trips = min(trips), max_trips = max(trips), mean_shift_km = round(mean(shift_km), 2)), by = cluster]
print(hourly_summary)
##    cluster min_trips max_trips mean_shift_km
##      <int>     <int>     <int>         <num>
## 1:       1        90      2468          0.88
## 2:       2        37       845          1.16
## 3:       3       279      1663          0.48
## 4:       4      1453      6628          0.36
## 5:       5      2266     29238          0.27
## 6:       6        24       238          2.95
## 7:       7       280      3177          0.49
## 8:       8      5072     32768          0.11
#hourly trip activity
ggplot(hourly_centers, aes(x = hour, y = factor(cluster), fill = trips)) +
  geom_tile() +
  scale_x_continuous(breaks = seq(0, 23, by = 2)) +
  labs(title = "Hourly uber trip activity by cluster", x = "hour", y = "cluster", fill = "trips") +
  theme_minimal()

#hourly center shifts
ggplot(hourly_centers, aes(x = hour, y = shift_km, color = factor(cluster))) +
  geom_line() + geom_point(size = 0.8) +
  facet_wrap(~cluster, scales = "free_y") +
  scale_x_continuous(breaks = seq(0, 23, by = 6)) +
  labs(title = "Hourly changes in cluster centers", x = "hour", y = "distance from original center (km)") +
  theme_minimal() + theme(legend.position = "none")

It shows uber trip activity varies across clusters throughout the day. Clusters 5 & 8 show the highest activity (peak at 17:00 and 18:00), while clusters 3 and 4 peak at 21:00.

The geographical centers also show different levels of hourly variation. Cluster 8 keeps stable (average 0.11km), while cluster 6 has the largest variation (2.95 km), probably due to its smaller number of observations.

3.2. By date

Similar to section (a), we analyze daily uber trip activity and changes in cluster centers from September 1 to 30, 2014.

#daily cluster centers
daily_centers <- uber[, .(trips = .N, center_x = mean(x), center_y = mean(y)), by = .(cluster, date)][order(cluster, date)]

#distance from the original cluster center
daily_centers[, shift_km := sqrt((center_x - model_8$centers[cluster, 1])^2 + (center_y - model_8$centers[cluster, 2])^2) / 1000]

#identify the busiest date for each cluster
peak_dates <- daily_centers[, .SD[which.max(trips)], by = cluster][, .(cluster, peak_date = date, peak_trips = trips)]
print(peak_dates)
##    cluster  peak_date peak_trips
##      <int>     <Date>      <int>
## 1:       1 2014-09-01       1530
## 2:       2 2014-09-01        445
## 3:       3 2014-09-27       1133
## 4:       4 2014-09-13       5686
## 5:       5 2014-09-16      15027
## 6:       6 2014-09-27        176
## 7:       7 2014-09-21       1929
## 8:       8 2014-09-13      20725
#summarize daily activity and center shifts
daily_summary <- daily_centers[, .(total_trips = sum(trips), mean_daily_trips = round(mean(trips), 1), mean_shift_km = round(mean(shift_km), 2)), by = cluster]
print(daily_summary)
##    cluster total_trips mean_daily_trips mean_shift_km
##      <int>       <int>            <num>         <num>
## 1:       1       29147            971.6          0.39
## 2:       2       10203            340.1          0.93
## 3:       3       24930            831.0          0.41
## 4:       4      108247           3608.2          0.15
## 5:       5      359839          11994.6          0.11
## 6:       6        3452            115.1          2.89
## 7:       7       43268           1442.3          0.31
## 8:       8      449050          14968.3          0.06
#daily trip activity
ggplot(daily_centers, aes(x = date, y = trips, color = factor(cluster))) +
  geom_line() + geom_point(size = 0.6) +
  facet_wrap(~cluster, scales = "free_y") +
  scale_x_date(date_breaks = "7 days", date_labels = "%d %b") +
  labs(title = "Daily uber trip activity by cluster", x = "date", y = "number of trips") +
  theme_minimal() + theme(legend.position = "none")

#daily center shifts
ggplot(daily_centers, aes(x = date, y = shift_km, color = factor(cluster))) +
  geom_line() + geom_point(size = 0.6) +
  facet_wrap(~cluster, scales = "free_y") +
  scale_x_date(date_breaks = "7 days", date_labels = "%d %b") +
  labs(title = "Daily changes in cluster centers", x = "date", y = "Distance from Original Center (km)") +
  theme_minimal() + theme(legend.position = "none")

It shows recurring fluctuations, with clusters 5 & 8 accounting for most observations. 8 has the most stable geographical center (0.06 km average shift), while 6 shows the largest variation (2.89 km), likely due to fewer observations.

4. Apply the dataset for forecasting

4.1. Data preparation

For this data preparation step, we prepare hourly uber trip counts for each of the 8 geographical clusters and split the data into training (September 1–23) and testing (September 24–30) sets. Missing cluster-hour combinations are filled with 0, while hrs not observed in the original dataset are excluded to avoid misleading forecast evaluations.

#count hourly trips for each cluster
trip_counts<-uber[, .(trips = .N), by = .(cluster, date, hour)]

#create a complete hourly grid for all clusters
observed_hours <- unique(uber[, .(date, hour)])
forecast_data <- CJ(cluster = sort(unique(uber$cluster)), date = sort(unique(uber$date)), hour = 0:23)

#exclude hours missing from the original dataset
forecast_data <- merge(forecast_data, observed_hours, by = c("date", "hour"))

#merge trip counts and fill missing cluster-hour counts with 0
forecast_data <- merge(forecast_data, trip_counts, by = c("cluster", "date", "hour"), all.x = TRUE)
forecast_data[is.na(trips), trips := 0L]
setorder(forecast_data, cluster, date, hour)

#split data into training and testing
train_data <- forecast_data[date < as.Date("2014-09-24")]
test_data <- forecast_data[date >= as.Date("2014-09-24")]

#summarize training and testing
split_summary <- rbindlist(list(
  train_data[, .(period = "Train", first_date = min(date), last_date = max(date), hourly_rows = .N, trips = sum(trips))],
  test_data[, .(period = "Test", first_date = min(date), last_date = max(date), hourly_rows = .N, trips = sum(trips))]
))
print(split_summary)
##    period first_date  last_date hourly_rows  trips
##    <char>     <Date>     <Date>       <int>  <int>
## 1:  Train 2014-09-01 2014-09-23        4416 789372
## 2:   Test 2014-09-24 2014-09-30        1336 238764
#check hourly coverage
print(forecast_data[, .(observed_hours = uniqueN(hour)), by = date][observed_hours < 24])
##          date observed_hours
##        <Date>          <int>
## 1: 2014-09-30             23
#check aggregated trip counts
stopifnot(sum(forecast_data$trips) == nrow(uber), !anyNA(forecast_data$trips))

Training (789,372 trips) and testing (238,764 trips). Hourly observations were prepared for all 8 clusters, with zero-trip hours included where appropriate. The incomplete hr on September 30 was excluded to ensure valid forecast evaluation.

4.2. Cluster-based forecasting model

We forecast hourly uber trip counts for each cluster using 2 methods (Seasonal Naive and Historical Hourly Profile). The first one uses observations from the same hr of the last week, while the second uses historical averages by cluster, weekday and hr (both use only training data to predict trips in the testing period).

#seasonal naive forecast using the previous week
naive_lookup <- train_data[, .(cluster, date = date + 7L, hour, pred_naive = trips)]

#add weekday information
train_data[, weekday := lubridate::wday(date, week_start = 1)]
test_data[, weekday := lubridate::wday(date, week_start = 1)]

#historical hourly averages
hourly_profile <- train_data[, .(pred_profile = mean(trips)), by = .(cluster, weekday, hour)]

#generate forecasts for the testing period
forecast_results <- merge(test_data, naive_lookup, by = c("cluster","date","hour"), all.x = TRUE)
forecast_results <- merge(forecast_results, hourly_profile, by = c("cluster","weekday","hour"), all.x = TRUE)
forecast_results[, pred_profile := round(pred_profile)]
setorder(forecast_results, cluster, date, hour)

#check results
stopifnot(nrow(forecast_results) == nrow(test_data), !anyNA(forecast_results[, .(pred_naive, pred_profile)]))

#show actual vs predicted trips
print(head(forecast_results[, .(cluster, date, hour, trips, pred_naive, pred_profile)], 12))
##     cluster       date  hour trips pred_naive pred_profile
##       <int>     <Date> <int> <int>      <int>        <num>
##  1:       1 2014-09-24     0     8          6           10
##  2:       1 2014-09-24     1     4          3            2
##  3:       1 2014-09-24     2     2          0            1
##  4:       1 2014-09-24     3     2          3            4
##  5:       1 2014-09-24     4    26         13           16
##  6:       1 2014-09-24     5    36         31           42
##  7:       1 2014-09-24     6    39         32           44
##  8:       1 2014-09-24     7    23         26           21
##  9:       1 2014-09-24     8    19         15           17
## 10:       1 2014-09-24     9    21         18           20
## 11:       1 2014-09-24    10    21         20           24
## 12:       1 2014-09-24    11    23         20           23

Both methods generated hourly trip forecasts for all 8 clusters during the testing. All predictions are available, and they use only historical training data.

4.3. Visualization

We’re going to use MAE, RMSE and WAPE to evaluate the forecasting accuracy of Seasonal Naive and Historical Hourly Profile. Then, compare actual and predicted hourly trip counts to see which method performs better during the testing.

#forecasting error metrics
evaluate<-function(pred, label) {
  err <- forecast_results$trips - pred
  data.table(method = label, MAE = round(mean(abs(err)), 2),
             RMSE = round(sqrt(mean(err^2)), 2),
             WAPE = round(100 * sum(abs(err)) / sum(forecast_results$trips), 2))
}

#compare forecasting methods
evaluation <- rbindlist(list(
  evaluate(forecast_results$pred_naive, "Seasonal Naive"),
  evaluate(forecast_results$pred_profile, "Historical hourly profile")
))
print(evaluation)
##                       method   MAE  RMSE  WAPE
##                       <char> <num> <num> <num>
## 1:            Seasonal Naive 28.46 66.95 15.93
## 2: Historical hourly profile 25.30 55.02 14.16
knitr::kable(evaluation)
method MAE RMSE WAPE
Seasonal Naive 28.46 66.95 15.93
Historical hourly profile 25.30 55.02 14.16 
#aggregate actual and predicted trips across clusters
plot_forecast <- forecast_results[, .(Actual = sum(trips),
  Seasonal_Naive = sum(pred_naive),
  Historical_Profile = sum(pred_profile)), by = .(date, hour)]

#create hourly timestamps
plot_forecast[, datetime := as.POSIXct(paste(date, sprintf("%02d:00:00", hour)),
  format = "%Y-%m-%d %H:%M:%S", tz = "America/New_York")]

#prepare data
plot_long <- melt(plot_forecast, id.vars = "datetime",
  measure.vars = c("Actual", "Seasonal_Naive", "Historical_Profile"),
  variable.name = "method", value.name = "trips")
## Warning in melt.data.table(plot_forecast, id.vars = "datetime", measure.vars =
## c("Actual", : 'measure.vars' [Actual, Seasonal_Naive, Historical_Profile, ...]
## are not all of the same type. By order of hierarchy, the molten data value
## column will be of type 'double'. All measure variables not of type 'double'
## will be coerced too. Check DETAILS in ?melt.data.table for more on coercion.
#visualize actual vs predicted trips
ggplot(plot_long, aes(x = datetime, y = trips, color = method)) +
  geom_line(linewidth = 0.6) +
  scale_x_datetime(date_breaks = "1 day", date_labels = "%d %b") +
  labs(title = "Actual vs Predicted hourly uber trips",
       x = "Date", y = "Number of trips", color = "Method") +
  theme_minimal()

The Historical Hourly Profile performs better than Seasonal Naive, with lower MAE (25.30), RMSE (55.02) and WAPE (14.16%). Both methods capture hourly trip patterns, but forecasting accuracy is limited by the short period.