Tidy Tuesday: French Train Schedule Analysis with Tidyverts

Creating a Tsibble

Creating a tsibble object requires a key and index variable to be set. The key field contains a unique variable combination for each observation and the index field represents the time component of the series. Earlier I used the tsibble::yearmonth() function to specify that this is a monthly time series. Once in tsibble format, it’s easy to visualize, extract features, or forecast a time series data set with Tidyverts packages.

# tidyverts  --------------------------------------------------------------

library(tsibble)
library(fabletools)
library(feasts)

trains_ts <- trains %>%
  tsibble::tsibble(key = c(service, departure_station, arrival_station), index = yearmonth)

## Viz of all time series 
trains_ts %>%
autoplot(total_num_trips) +
  guides(col = FALSE) +
  labs(y = "Trips", x = NULL, title = "Monthly Total Trips by Train Leg",
       subtitle = "Each line represents a unique departure and arrival station combination")

Time Series Feature Extraction

The above graph shows us the time frame of our time series data set and the variation in trip volume among different train routes. Let’s extract some features from these time series in order to better understand differences among them. The fabletools::features() function does just this, transforming the tsibble object into a tibble with one row for each key. I used a set of features pre-defined by the feasts package, but a custom list can be created as well.

note: some features were not appropriate for this data set and had zero variance or NA results. I added an additional processing step to remove these

## Feature Extraction

train_features <- trains_ts %>%
  fabletools::features(total_num_trips, fabletools::feature_set(pkgs = "feasts")) 

## Identify columns with NA values and 0 variance
non_zero_var <- train_features %>%
  summarise(across(everything(), ~var(.x, na.rm = T))) %>%
  pivot_longer(cols = everything(), names_to = "column", values_to = "value") %>% 
  filter(value != 0 & !is.na(value)) %>%
  pull(column)

## Remove columns with NA values 0 variance
train_features <- train_features %>%
  select(service, departure_station, arrival_station, all_of(non_zero_var)) %>%
  mutate(across(everything(), ~replace_na(.x, 0)))

head(train_features)

## # A tibble: 6 x 42
##   service departure_stati~ arrival_station trend_strength spikiness linearity
##   <chr>   <chr>            <chr>                    <dbl>     <dbl>     <dbl>
## 1 Intern~ Barcelona        Paris Lyon               0.572   1686.       41.0 
## 2 Intern~ Francfort        Paris Est                0.663    438.      159.  
## 3 Intern~ Geneve           Paris Lyon               0.205   2137.      -41.3 
## 4 Intern~ Italie           Paris Lyon               0.307     50.7      -3.55
## 5 Intern~ Lausanne         Paris Lyon               0.570    147.       81.1 
## 6 Intern~ Madrid           Marseille St C~          0.611      1.46      9.32
## # ... with 36 more variables: curvature <dbl>, stl_e_acf1 <dbl>,
## #   stl_e_acf10 <dbl>, seasonal_strength_year <dbl>, seasonal_peak_year <dbl>,
## #   seasonal_trough_year <dbl>, acf1 <dbl>, acf10 <dbl>, diff1_acf1 <dbl>,
## #   diff1_acf10 <dbl>, diff2_acf1 <dbl>, diff2_acf10 <dbl>, season_acf1 <dbl>,
## #   pacf5 <dbl>, diff1_pacf5 <dbl>, diff2_pacf5 <dbl>, season_pacf <dbl>,
## #   nonzero_squared_cv <dbl>, lambda_guerrero <dbl>, nsdiffs <dbl>,
## #   bp_stat <dbl>, bp_pvalue <dbl>, lb_stat <dbl>, lb_pvalue <dbl>,
## #   var_tiled_var <dbl>, var_tiled_mean <dbl>, shift_level_max <dbl>,
## #   shift_level_index <dbl>, shift_var_max <dbl>, shift_var_index <dbl>,
## #   shift_kl_max <dbl>, shift_kl_index <dbl>, spectral_entropy <dbl>,
## #   n_crossing_points <dbl>, longest_flat_spot <dbl>, stat_arch_lm <dbl>

Now that we reduced the time series to a tidy data frame of features, we can work with our original data in a more traditional way. Let’s see which routes have a greatest overall trend strength (increase over the entire time period) and yearly seasonality.

train_features %>%
  mutate(train_leg = paste0(departure_station," to ", arrival_station)) %>%
  ggplot(aes(x = seasonal_strength_year, y = trend_strength)) +
  geom_point() +
  ggrepel::geom_text_repel(aes(label = if_else(seasonal_strength_year >= .87 | trend_strength >= .8, train_leg, NULL))) +
  labs(x = "Yearly Seasonal Strength", y = "Overall Trend Strength",
       title = "Trend vs Seasonality for All Train Legs",
       subtitle = "Higher values labelled for context")

Extracting features from the time series lets us easily answer certain questions compared to sifting through hundreds of time series graphs. The below graphs use simple filtering to answer questions I had about the data set.

## Most positive trend

trains_ts %>%
  inner_join(train_features %>% filter(trend_strength == max(trend_strength)), by = c("service","departure_station","arrival_station")) %>%
  ggplot(aes(x = yearmonth, y = total_num_trips)) + 
  geom_line() +
    labs(title = "Which Route Frequency has the Greatest Overall Trend Strength?", y = "Trips", x = NULL) +
  facet_wrap(vars(service, paste0("Departure: ",departure_station), paste0("Arrival: ",arrival_station)))

## Most Seasonal Trend

trains_ts %>%
  inner_join(train_features %>% filter(seasonal_strength_year == max(seasonal_strength_year, na.rm = T)), by = c("service","departure_station","arrival_station")) %>%
  ggplot(aes(x = yearmonth, y = total_num_trips)) + geom_line() +
  labs(title = "Which Route has the Most Seasonal Schedule?", y = "Trips", x = NULL) +
  facet_wrap(vars(service, paste0("Departure: ",departure_station), paste0("Arrival: ",arrival_station)))

## Highest spikiness 

trains_ts %>%
  inner_join(train_features %>% filter(spikiness == max(spikiness, na.rm = T)), by = c("service","departure_station","arrival_station")) %>%
  ggplot(aes(x = yearmonth, y = total_num_trips)) + geom_line() +
  labs(title = "Which Route has the Most Drastic Changes in Volumes?", y = "Trips", x = NULL) +
  facet_wrap(vars(service, paste0("Departure: ",departure_station), paste0("Arrival: ",arrival_station)))

PCA Analysis

Feature extraction lets us expand our exploratory analysis toolkit for time series data. I decided to perform principal component analysis on the train features to better understand the main sources of variation among monthly train volumes.

# Principal Component Analysis on Time Series Features --------------------

features_pcs <- train_features %>%
  select(-service, -departure_station, -arrival_station) %>%
  prcomp(scale = TRUE)


features_pcs$rotation %>% 
  as.data.frame() %>%
  rownames_to_column() %>% 
  pivot_longer(cols = -rowname, names_to = "PC") %>%
  group_by(PC) %>%
  #top_n(30, abs(value)) %>%
  filter(PC %in% c("PC1","PC2","PC3")) %>%
  mutate(rowname_ro = tidytext::reorder_within(rowname, value, PC)) %>%
  ggplot(aes(x = value, y = rowname_ro, fill = rowname)) +
  geom_col() +
  tidytext::scale_y_reordered() +
  facet_wrap(~PC, scales = "free_y") +
  guides(fill = FALSE) +
  labs(title = "What Features Account for the Top 3 Principal Components?",
        y = NULL)

Clustering Analysis

Let’s do the same thing with k-means clustering to group our different train routes and see what characterizes each group. Below I perform a traditional k-means clustering analysis, look at the main differences between each cluster, and join the cluster number back to the train features and our original data set.

# K means clustering ------------------------------------------------------

set.seed(456)

## Scale train features
train_features_scaled <- train_features %>%
  mutate(across(where(is.numeric), ~scale(.x)))

## Cluster and store results for K 1 through 9
train_cluster <- tibble(k = 1:9) %>%
  mutate(
    kclust = map(k, ~kmeans(train_features_scaled %>% select(where(is.numeric)), .x)),
    tidied = map(kclust, tidy),
    glanced = map(kclust, glance),
    augmented = map(kclust, augment, train_features_scaled %>% select(where(is.numeric)))
    )

## Visualize results to detemine optimal clusters to use
train_cluster %>%
  unnest(cols = c(glanced)) %>%
  ggplot(aes(k, tot.withinss)) +
  geom_line() +
  geom_point() +
  scale_x_continuous(breaks = seq(1,9,1)) +
  labs(title = "Elbow Analysis of Train Route Feature Clusters",
       subtitle = "Suggests 4 is the optimal number of clusters")

## Add cluster column to train features
train_features_clusters <- train_features %>%
  mutate(cluster = train_cluster %>%
           filter(k == 4) %>%
           unnest(augmented) %>%
           pull(.cluster))

If we facet our original time series graph by cluster number we can begin to see the main differences between our groupings. There looks to be a set of new train lines that began operating in 2018 that make up an entie cluster.

trains_ts %>%
  left_join(train_features_clusters) %>%
  ggplot(aes(x = yearmonth, y = total_num_trips, col = cluster, group = interaction(service, departure_station, arrival_station))) +
  geom_line() +
  facet_wrap(~cluster) +
  labs(title = "Train Route Volumes Clustered", x = NULL, y = "Total Trips")

Below are those feature differences between clusters quantified. By graphing a smaller sample of routes from each cluster it’s easier to see the main features represented in the data.

Here are a few of my initial observations:

• Cluster 1 - Most seasonal French train routes. Possibly tied to vacation travel patterns.
• Cluster 2 - Has the highest overall trends. Route volumes tend to steadily increase or decrease over time and with little seasonal pattern.
• Cluster 3 - New train lines that began operating in 2018. Many appear to be international routes.
• Cluster 4 - Flat overall trend with drastic variations indicated by downward spikes. Google searches for some of these spikes coincide with train worker strikes and other service outage events.

## Bar chart of cluster differences
train_cluster %>%
  filter(k == 4) %>%
  unnest(cols = c(tidied)) %>%
  select(-kclust,-glanced, -augmented, -k) %>% 
  pivot_longer(cols = trend_strength:stat_arch_lm) %>%
  group_by(cluster) %>%
  top_n(20, abs(value)) %>%
  #mutate(name_ro = tidytext::reorder_within(name, value, cluster)) %>%
  ggplot(aes(x = value, y = name, fill = name)) +
  geom_col() +
  #tidytext::scale_y_reordered() +
  facet_wrap(~cluster, nrow = 1) +
  guides(fill = FALSE) +
  labs(x = "Scaled Value", y = NULL, title = "What Features Define each Cluster?")

## Sample routes from each cluster 
set.seed(1234)

trains_ts %>%
  inner_join(train_features_clusters %>%
               group_by(cluster) %>%
               slice_sample(n = 6)) %>%
  ggplot(aes(x = yearmonth, y = total_num_trips, col = cluster, group = interaction(service, departure_station, arrival_station))) +
  geom_line(size = 1) +
  facet_wrap(~cluster) +
  labs(title = "Train Route Samples by Cluster ", x = NULL, y = "Total Trips")

Combining PCA and Clustering

I thought it’d be interesting to compare the PCA and clustering results to see if there are any discrepancies. After combining the results, I created a scatter plot of routes in PC space and colored them by cluster. The results appear very consistent, with the first 2 principal components graph clearly grouping clusters. This also shows that cluster 3 differs most from all others and clusters 4 and 1 appear most alike with overlapping points.

### Combines Graph of PC and Clustering


features_pcs %>% 
  broom::augment(train_features) %>%
  left_join(train_features_clusters %>% select(service, departure_station, arrival_station, cluster)) %>%
  ggplot(aes(x=.fittedPC1, y=.fittedPC2, col = cluster)) +
  geom_point(size = 2) + 
  theme(aspect.ratio=1) +
  labs(title = " Train Routes in PC Space")