Taxi Trip ETA
2023 — A model that predicts how many minutes a taxi ride will take
01 Problem
A rider who knows when the taxi will arrive waits less, a dispatcher who knows how long each trip takes assigns cars better, and a fleet that stops guessing burns less fuel and idle time.
The goal is a model that predicts the duration of a New York green taxi trip at the moment it is booked, so the arrival time can be promised before the ride starts.
02 Cleaning
The raw data holds 83,691 trips from the NYC Taxi and Limousine Commission, with no duplicated rows. Columns with missing values and columns that only break down the fare, such as tax, tip, and tolls, are dropped, since none of them is known when a trip is booked.
The summary exposes three kinds of impossible records, which are trips that travelled zero distance, one trip of 260,517 miles, and negative fares. Trips shorter than 0.1 mile, fares under $1, and distance outliers are removed.
trip <- trip_raw %>%
select_if(~ all(!is.na(.))) %>%
select(-extra,
-mta_tax,
-tip_amount,
-improvement_surcharge,
-total_amount,
-store_and_fwd_flag,
-tolls_amount) %>%
mutate_if(is.character, ymd_hms) %>%
rename(PU_datetime = lpep_pickup_datetime,
DO_datetime = lpep_dropoff_datetime) %>%
mutate_at(vars(PULocationID, DOLocationID), as.factor)
trip <- trip %>%
filter(trip_distance >= 0.1,
fare_amount >= 1) %>%
killOutliers(trip_distance)03 Features
Three features come out of the timestamps. They are the trip duration in minutes, which becomes the target, the pickup hour, and whether the trip starts on a Saturday or Sunday.
Fare amount is dropped. It is highly correlated with distance, and a fare is only settled after the ride ends, so keeping it would hand the model the answer. Durations under one minute and duration outliers, the longest at nearly 24 hours, are removed as well.
trip <- trip %>%
mutate(trip_duration = as.numeric(round(difftime(DO_datetime, PU_datetime, units = "mins"), 2)),
is_weekend = as.factor(ifelse(wday(PU_datetime, week_start = 1) > 5, 1, 0)),
hour = as.factor(hour(PU_datetime))) %>%
select(-PU_datetime,
-DO_datetime,
-fare_amount) %>%
killOutliers(trip_duration) %>%
filter(trip_duration >= 1)04 Duration
After cleaning, 69,715 trips remain. Most rides are short, and the most common trip lasts about 8 minutes.
ggplot(trip, aes(x = trip_duration)) +
geom_density(fill = "#56b4fc", col = "white", alpha = 0.5) +
labs(x = "Trip Duration",
y = "Density",
title = "Distribution of Trip Duration") +
theme_minimal()05 Demand
Demand climbs from 6 in the morning, stays high from 8 to 18, and fades through the night, which is why the pickup hour is kept as a feature.
plot_data <- trip %>%
group_by(hour, is_weekend) %>%
summarise(count = n()) %>%
ungroup()
ggplot(plot_data, aes(x = hour, y = count)) +
geom_bar(stat = "identity", aes(fill = count), show.legend = F) +
facet_grid(. ~ is_weekend,
labeller = labeller(is_weekend = c("0" = "Weekday", "1" = "Weekend"))) +
theme_minimal()06 Split
The trips are split 80/20 into 55,774 for training and 13,941 for testing. The test set is never used for training, validation, or stacking, so every score below is measured on trips the models have never seen.
set.seed(1)
indices <- createDataPartition(trip$trip_duration,
p = 0.8,
list = FALSE)
train_data <- trip[indices, ] %>% select(-PULocationID, -DOLocationID)
test_data <- trip[-indices, ] %>% select(-PULocationID, -DOLocationID)07 Baselines
A linear regression sets the floor at 4.11 minutes of mean absolute error. A random forest of 55 trees, which can bend around non-linear traffic patterns, brings the percentage error down from 36.8% to 28.9%.
model_lm <- lm(formula = trip_duration ~ .,
data = train_data)
model_rf <- randomForest(x = X_train,
y = y_train,
ntree = 55,
mtry = ncol(X_train),
importance = T)08 Neural network
The pickup and dropoff zones are one-hot encoded into 535 features, which lets a feed-forward network learn how long it takes to cross between specific parts of the city.
Three hidden layers of 256, 128, and 64 ReLU units feed a single linear output, trained with Adam on mean squared error. A fifth of the training data is held back for validation, and the network cuts the error to 3.21 minutes.
The validation error keeps falling alongside the training error through all 18 epochs, so the network is still learning rather than memorizing when training stops.
model_fnn <- keras_model_sequential(name = "ETA-Predictor_NYC-Taxi_1.0.0") %>%
layer_dense(units = 256, activation = "relu", input_shape = ncol(X_train), name = "hidden_1") %>%
layer_dense(units = 128, activation = "relu", name = "hidden_2") %>%
layer_dense(units = 64, activation = "relu", name = "hidden_3") %>%
layer_dense(units = 1, activation = "linear", name = "output")
model_fnn %>% compile(
loss = "mean_squared_error",
optimizer = optimizer_adam(learning_rate = 0.01),
metrics = "mean_absolute_percentage_error"
)
tensorflow::set_random_seed(1)
history_fnn <- model_fnn %>% fit(
X_train, y_train,
epochs = 18,
batch_size = 8192,
validation_split = 0.2
)09 Stacking
A k-nearest-neighbour regressor averages the durations of the five most similar trips and scores 3.40 minutes on its own. Its prediction becomes an extra input for a second network.
Each training row gets its KNN prediction from a model fitted on the other four folds, so a trip is never its own neighbour and the network cannot learn to trust a hint that already contains the answer.
set.seed(1)
folds <- createFolds(y_train, k = 5)
pred_knn_train <- numeric(length(y_train))
for (fold in folds) {
pred_knn_train[fold] <- knn.reg(train = X_train[-fold, ],
test = X_train[fold, ],
y = y_train[-fold],
k = 5)$pred
}
pred_knn_test <- knn.reg(train = X_train,
test = X_test,
y = y_train,
k = 5)$pred
X_train <- cbind(X_train, pred_knn = pred_knn_train)
X_test <- cbind(X_test, pred_knn = pred_knn_test)10 Selection
The stacked network wins on every metric, with 3.09 minutes of mean absolute error, 22.4% mean absolute percentage error, and 4.50 minutes of root mean squared error.
11 Assumptions
Predictions follow the actual durations along the diagonal, and the residuals bounce randomly around zero with no visible pattern, so linearity and homoscedasticity hold. The residual density peaks at zero.
Every feature has a variance inflation factor close to 1, far below the threshold of 10, so no feature is a disguised copy of another.
residuals_fnn_knn <- y_test - pred_fnn_knn
ggplot(data.frame(y_test, pred_fnn_knn), aes(x = y_test, y = pred_fnn_knn)) +
geom_point() +
geom_abline(intercept = 0, slope = 1, color = "red") +
theme_minimal()
vif(model_lm)Taxi Trip ETA
Lands within 3.1 minutes of the real duration on average