Import the Data

I begin by equipping myself with essential libraries like tidyverse, providing me with the ability to model, transform, and visualize data throughout the journey.

library(tidyverse)
library(skimr)
library(scales)

Then, I start importing all 12 monthly datasets and merging them into a single dataset named trip_data. To achieve this, I leverage two key functions:

1. list.files to generate a list of all file paths for efficient reading, and
2. read_csv to iteratively read each file into memory.

Once all files are read, I use the bind_rows function to seamlessly stitch them together into a unified dataset.

trip_data <- list.files(path = "./data/", pattern = "*-divvy-tripdata.csv", full.names = TRUE) %>% 
  lapply(read_csv) %>% 
  bind_rows %>% 
  arrange(started_at)

I would like to take a quick look with the dataset. Cyclists pedaled their way to more than 5 million trips on Cyclistic bikes in 2022, demonstrating the significant impact of bike-sharing programs.

dim(trip_data)

## [1] 5667717      13

head(trip_data)

Remove Duplicates

I use the duplicated and sum functions to meticulously check for duplicate trip data. The results reveal no duplicate values, indicating excellent data integrity and providing a reliable foundation for further analysis.

sum(duplicated(trip_data$ride_id))

## [1] 0

Correct Structural Errors

The next step is to address structural errors, those inconsistencies in data formatting that can lead to misinterpretations. This often takes the form of typos, inconsistent capitalization, or duplicate entries. With this in mind, I’ll primarily focus on cleaning up the member type, bike type, and start and end station columns, as they contain text data and are particularly susceptible to such issues.

member_type <- count(trip_data, member_casual, name = "count")
member_type

bike_type <- count(trip_data, rideable_type, name = "count")
bike_type

trip_data_v2 <- trip_data %>% 
  mutate(rideable_type = str_replace_all(rideable_type, "docked_bike", "classic_bike"))

bike_type_v2 <- count(trip_data_v2, rideable_type, name = "count")
bike_type_v2

start_station <- trip_data_v2 %>% 
  count(start_station_name, name = "count") %>% 
  arrange(start_station_name)

start_station

test_station_list <- c("Pawel Bialowas - Test- PBSC charging station", 
                       "Hastings WH 2", 
                       "DIVVY CASSETTE REPAIR MOBILE STATION", 
                       "Base - 2132 W Hubbard Warehouse", 
                       "Base - 2132 W Hubbard", 
                       "NewHastings", 
                       "WestChi", 
                       "WEST CHI-WATSON")

trip_data_v2 <- trip_data_v2 %>% 
  filter(!(trip_data_v2$start_station_name %in% test_station_list | 
           trip_data_v2$end_station_name %in% test_station_list))

words <- c("*", " - Charging", " (Temp)", "amp;", "Public Rack - ", 
           " - north corner", " - south corner", " - midblock south", " - midblock", 
           " - North", " - South", " - East", " - West", 
           " - NE", " - NW", " - SE", " - SW", 
           " - N", " - S", " - E", " - W", 
           " NE", " NW", " SE", " SW")

for (word in words) {
  trip_data_v2 <- trip_data_v2 %>% 
    mutate(start_station_name = str_replace_all(start_station_name, fixed(word, ignore_case = TRUE), "")) %>%
    mutate(end_station_name   = str_replace_all(end_station_name, fixed(word, ignore_case = TRUE), ""))
}

trip_data_v2 <- trip_data_v2 %>% 
  mutate(start_station_name = str_replace_all(start_station_name, regex(" (?<=\\s)[N|S|E|W]$", ignore_case = TRUE), "")) %>% 
  mutate(end_station_name   = str_replace_all(end_station_name, regex(" (?<=\\s)[N|S|E|W]$", ignore_case = TRUE), ""))

start_station_v2 <- trip_data_v2 %>% 
  count(start_station_name, name = "count") %>% 
  arrange(start_station_name)

start_station_v2

Handle Missing Data

I use the colSum and is.na functions to meticulously identify missing values within the dataset. The analysis reveals six columns with missing data: start_station_name (833,025 rows), start_station_id (833,025 rows), end_station_name (891,896 rows), end_station_id (891,896 rows), end_lat (5,858 rows), and end_lng (5,858 rows).

colSums(is.na(trip_data_v2))

##            ride_id      rideable_type         started_at           ended_at 
##                  0                  0                  0                  0 
## start_station_name   start_station_id   end_station_name     end_station_id 
##             833025             833025             891896             891896 
##          start_lat          start_lng            end_lat            end_lng 
##                  0                  0               5858               5858 
##      member_casual 
##                  0

To address these missing values, I prioritize start_station_name and end_station_name, as their missing entries can be imputed using geographic coordinates as a reference. I refrain from imputing the remaining four columns: start_station_id and end_station_id will be removed in the subsequent step due to their irrelevance to the analysis, while end_lat and end_lng lack suitable reference data for imputation and will also be discarded.

start_station_location <- trip_data_v2 %>% 
  count(start_lat, start_lng , start_station_name, name = "count") %>% 
  arrange(start_lat, start_lng)

start_station_location

To effectively handle missing station names, I create four temporary columns representing increasingly granular levels of start and end coordinates (ranging from 5 digits to 1 digit). These columns facilitate accurate matching and assignment of station names to missing entries, prioritizing those with more precise location data. This iterative approach ensures that the most reliable matches are filled in first.

digit <- 5

while (digit > 1) {
  trip_data_v2 <- trip_data_v2 %>% 
    mutate(start_lat_round = round(start_lat, digits = digit), 
           start_lng_round = round(start_lng, digits = digit), 
           end_lat_round   = round(end_lat,   digits = digit), 
           end_lng_round   = round(end_lng,   digits = digit))
  
  trip_data_v2 <- trip_data_v2 %>% 
    group_by(start_lat_round, start_lng_round) %>% 
    fill(start_station_name, .direction = "downup") %>% 
    fill(start_station_id,   .direction = "downup") %>% 
    ungroup()
  
  trip_data_v2 <- trip_data_v2 %>% 
    group_by(end_lat_round, end_lng_round) %>%  
    fill(end_station_name, .direction = "downup") %>% 
    fill(end_station_id,   .direction = "downup") %>% 
    ungroup()
  
  digit <- digit - 1
}

trip_data_v2 <- trip_data_v2 %>% 
  select(!c(start_lat_round, start_lng_round, end_lat_round, end_lng_round))

The imputation process yields a significant reduction in missing values: start_station_name dropped by 11,821 records each, while end_station_name saw a decrease of 41,386 records each.

colSums(is.na(trip_data_v2))

##            ride_id      rideable_type         started_at           ended_at 
##                  0                  0                  0                  0 
## start_station_name   start_station_id   end_station_name     end_station_id 
##              11821              11821              41386              41386 
##          start_lat          start_lng            end_lat            end_lng 
##                  0                  0               5858               5858 
##      member_casual 
##                  0

start_station_location_v2 <- trip_data_v2 %>% 
  count(start_lat, start_lng, start_station_name, name = "count") %>% 
  arrange(start_lat, start_lng)

start_station_location_v2

To ensure a focus on complete trips, I subsequently remove all remaining rows with missing values using the drop_na function. This results in 5,614,669 remaining trips, representing a relatively minor reduction from the original dataset.

trip_data_v2 <- drop_na(trip_data_v2)
dim(trip_data_v2)

## [1] 5614669      13

Add Relevant Columns

To facilitate time-series analysis, I create three new columns: ride_length_min, day_of_week, and month.

trip_data_v2$ride_length_min <- as.double(difftime(trip_data_v2$ended_at, trip_data_v2$started_at, units = "mins"))
trip_data_v2$day_of_week <- wday(trip_data_v2$started_at, label = TRUE)
trip_data_v2$month <- format(trip_data_v2$started_at, "%b")

However, both day_of_week and month are initially stored as unsorted text strings and not categorized in their natural order. Therefore, I utilize the ordered function to transform them based on their natural order (Monday to Sunday for days and January to December for months), ensuring their suitability for further analysis.

trip_data_v2 <- within(trip_data_v2, {
  day_of_week <- ordered(day_of_week, levels = c("Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"))
  month <- ordered(month, levels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))
})

Drop Unnecessary Data

The final step involves identifying and addressing outliers in the data. This is crucial to ensure the validity of the numbers and prevent distortion of the analysis. Outliers are identified using interquartile range (IQR) method, and their relevance is assessed based on their impact on the analysis and potential for being errors. If an outlier is deemed irrelevant or a mistake, it’s removed from the dataset. This process may involve dropping irrelevant columns, correcting error inputs, and filtering out potential outliers.

trip_data_v2 <- trip_data_v2 %>% 
  select( !c(start_station_id, end_station_id) )

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -10353.35      5.80     10.25     16.29     18.38  32035.45

trip_data_v2 <- trip_data_v2 %>% 
  filter( !(trip_data_v2$start_lat == 0 | trip_data_v2$start_lng == 0 | trip_data_v2$end_lat == 0 | trip_data_v2$end_lng == 0 | 
            trip_data_v2$ride_length_min < 1 | trip_data_v2$ride_length_min > 1440) )

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    6.05   10.48   16.55   18.67 1439.37

ggplot(data = trip_data_v2, aes(x = member_casual, y = ride_length_min, fill = member_casual)) + 
  geom_boxplot() +
  coord_flip() + 
  theme(legend.position = "none") + 
  labs(x = "Member type", 
       y = "Ride length (in minutes)", 
       title = "Box plot showing 'ride_length_min' before removing outliers")

quantiles <- as.numeric(quantile(trip_data_v2$ride_length_min, probs = c(0.25, 0.50, 0.75), na.rm = FALSE))
iqr_value <- IQR(trip_data_v2$ride_length_min)

## [1] "Q1: 6.05 minutes"

## [1] "Q3: 18.66667 minutes"

## [1] "IQR: 12.61667 minutes"

lower_fence <- quantiles[1] - ( 1.5 * iqr_value )
upper_fence <- quantiles[3] + ( 1.5 * iqr_value )

## [1] "Lower Fence: -12.875 minutes"

## [1] "Upper Fence: 37.59167 minutes"

trip_data_v2 <- trip_data_v2 %>%
  filter(!(trip_data_v2$ride_length_min < lower_fence | trip_data_v2$ride_length_min > upper_fence))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   5.783   9.717  11.958  16.150  37.583

ggplot(data = trip_data_v2, aes(x = member_casual, y = ride_length_min, fill = member_casual)) + 
  geom_boxplot() +
  coord_flip() + 
  theme(legend.position = "none") + 
  labs(x = "Member type", 
       y = "Ride length (in minutes)", 
       title = "Box plot showing 'ride_length_min' after removing outliers")

Validate the Data

After data cleaning, I re-examined the dataset, finding it now contains 5,091,142 trip records, a reduction of 576,575 (10%) compared to the original 5,667,717. There are 14 columns in the data, all of which have no missing values, ensuring a complete dataset for thorough analysis.

skim_without_charts(trip_data_v2)

Variable type: character

Variable type: factor

Variable type: numeric

Variable type: POSIXct

To verify if the data follows a normal distribution, I use the Empirical Rule, which states that 99.7% of data points in a normal distribution cluster within three standard deviations of the mean (68% within 1 SD, 95% within 2 SD). I create the calculate_percentage function to measure the data percentage within specific standard deviation ranges. Its output reveals that 72.10%, 94.08%, and 99.38% of the data fall within one, two, and three standard deviations of the mean, respectively, confirming a close resemblance to a normal distribution. This suggests the data can be reliably analyzed with various statistical methods.

summary_stats <- summarise(trip_data_v2,
                           sd = sd(ride_length_min),
                           mean = mean(ride_length_min),
                           count = n())

calculate_percentage <- function(n_sd) {
  filtered_count <- trip_data_v2 %>%
  filter(between(ride_length_min, 
                 summary_stats$mean - n_sd * summary_stats$sd, 
                 summary_stats$mean + n_sd * summary_stats$sd)) %>%
  summarise(count = n())
  round((filtered_count$count / summary_stats$count) * 100, 2)
}

percentage_sd1 <- calculate_percentage(1)
percentage_sd2 <- calculate_percentage(2)
percentage_sd3 <- calculate_percentage(3)

## [1] "One standard deviation: 72.10%"

## [1] "Two standard deviations: 94.08%"

## [1] "Three standard deviations: 99.38%"

Overall Summary of Cyclistic Rides in 2022

As shown in the below table, it presents a statistical summary of time duration for Cyclistic’s users, including casual riders and annual members.

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   5.783   9.717  11.958  16.150  37.583

The time duration of Cyclistic’s rides can also be visualized by a histogram. The histogram shows that the distribution of ride times is positively skewed, with a high cluster of lower values and a spread-out tail on the right. This means that there are a few rides that are much longer than the majority of rides.

Because of the positive skew, the mean ride time is not a very accurate measure of the typical ride time. The median ride time is a more accurate measure, as it is not as affected by the outliers. So, the median ride time will be chosen to conduct further descriptive analysis.

Key Findings

Based on the data collected, it can be concluded that

• Casual riders: They are more likely to use Cyclistic bikes for leisure purpose, including enjoyment, fitness, exploring, and much more. They spend longer periods of time and have a higher number of rides on weekends. They take a ride throughout the day, but have higher rides at the evening.
• Annual members: They tend to use Cyclistic bikes for specific purpose, such as commute to work. They spend shorter periods of time and have a higher number of rides on weekdays. They take a ride during rush hour , especially at 8 AM, 12 PM, and 5 PM.
• Warmer months: Both casual riders and annual members take more rides during Spring (in the beginning of March) and peak usage at the end of Summer (July and August).

Recommendations

The report makes the following recommendations:

1. Price Incentives: Offer incentives for casual riders to become annual members, such as discounts, free ride credits, or gift cards. This would make it a more attractive option for casual riders who are not sure if they will use Cyclistic bikes enough to justify the full annual membership price.
2. Personalized Marketing: Personalize the marketing messages that we send to casual riders. This will help to ensure that we are reaching them with the right message at the right time. For example, we could send a message to casual riders who have taken a certain number of rides in a month, or who have used Cyclistic bikes during rush hour.
3. New Member Plans: Consider to add a monthly or quarterly membership option. It would also be a good option for casual riders who only bike during certain times of the year, such as springtime and summertime.