report.Rmd

---
title: "Flight Delays: An Analysis of Patterns and Predictive Modeling"
output: 
  html_document:
    toc: true
    code_folding: hide
    toc_float: true
---


## Introduction

New York City, a global fusion of diverse culture and busy travel, is serviced by three major airports: John F. Kennedy International Airport (JFK), LaGuardia Airport (LGA), and Newark Liberty International Airport (EWR). These airports are not only pivotal in the United States' air traffic network but also play a significant role in international air travel. This report incorporates an in-depth predictive analysis of flight delays at these three airports, which reveal patterns and probabilities of such delays.This study seeks to provide a robust framework for predicting flight delays, and thus providing better planning and management guideline for passengers.

<div style="text-align: center;">
  <img src="image/nyc-airport-transportation.jpg" style="width: 100%;">
</div>

## Motivation

The airline industry is highly dependent on timely flight operations, and flight delays are a common concern for both passengers and airlines. Understanding and predicting these delays can lead to improving airlines efficiency and customer satisfaction. In addition, the New York City area's unique geographical and diverse cultural characteristics present an intriguing case for study. The region's congested airspace, fluctuated weather conditions, and intense flight schedules create an environment where delays are frequent. Analyzing this environment offers valuable insights into urban airport operations and delay management, potentially serving as a model for other major cities that share similar backgrounds as well. 

## Inspiration and Related Work
In 2022, during the Christmas as well as many year-end holidays, the Southwest Airlines had experienced an unprecedented meltdown. According to the article "Southwest meltdown may cost the airline up to $825 million" published by Chris Isidore(2023) on [CNN](https://www.cnn.com/2023/02/08/business/southwest-meltdown-hearings/index.html), shares of Southwest had already lost 8% since December 21, after the start of the incident. <br>
One crucial factor that caused this massive meltdown is due to the fact that the Southwest airline failed to prepare for the winter storm and accumulated uncontrollable cancellations. Therefore, we were surprised by the amount of money lost during the meltdown and wished to build a model to predict the effects of factors that may cause flights delayed, including weather conditions. 


## Initial Questions

For the purpose of the project, we will be answering the following questions: <br>
- Which factors may build a model for predicting flight delays? <br>
- Does our model based on data from 2013 could be potentially used to predict flight delays in 2017? <br>
- Can our model be generalized based on the reults of the previous question. <br>


## Data Source

In this project, we located several publicly available data sources shown as following: 

* `nycflights13` from R library, which includes flight information and weather conditions in 2013
*  The datasets concerning with flight and weather in 2017 can be accessed through this [link](https://github.com/jayleetx/nycflights/tree/master). The data in 2017 would be used for prediction purpose. 

Note: all relevant datasets for this study were contained in the folder `data`. 

## Data Cleaning 

Files and libraries required at this stage:

* `tidyverse`
* `dplyr`
* `library ("nycflights13")`
* `flights_2017.rda`
* `weather_2017.rda`

The weather condition and flight information would be merged based on the `origin`, `year`, `month`, `day` and `hour`. As we perceived that large portion of values (>0.9) in some columns  were missing in `2017 weather dataset`, we removed these columns on both 2013 and 2017 datasets. Lastly, we removed the records containing any NA values. . 

```{r, message=FALSE, warning=FALSE}
library(tidyverse)
library(dplyr)
library("nycflights13")
knitr::opts_chunk$set(echo = TRUE,
                      warning = FALSE,
                      fig.align = 'center',
                      message = FALSE)
flights_2013 = flights |> 
  janitor::clean_names()
weather_2013 = weather |> 
  janitor::clean_names()

load("data/flights_2017.rda")
load("data/weather_2017.rda")

flights_2017 = flights |> 
  janitor::clean_names()
weather_2017 = weather |> 
  janitor::clean_names()


flights_2013_clean = 
  flights_2013 |> 
  drop_na() |> 
  unique()  |> 
  select(
    -dep_time,
    -arr_time,
    -dep_delay,
    -sched_arr_time,
    -sched_dep_time,
    -time_hour
  )

weather_2013_clean = 
  weather_2013 |> 
  select(-c("temp", "dewp", "humid")) |> 
  drop_na() |> 
  select(-time_hour)

merge_data_2013 = merge(flights_2013_clean, weather_2013_clean, by =c("origin", 
                                                                      "year", 
                                                                      "month",
                                                                      "day",
                                                                    "hour"))


flights_2017_clean = 
  flights_2017 |> 
  drop_na() |> 
  unique()  |> 
  select(
    -dep_time,
    -arr_time,
    -dep_delay,
    -sched_arr_time,
    -sched_dep_time,
    -time_hour
  )

weather_2017_clean = 
  weather_2017 |> 
  select(-c("temp", "dewp", "humid")) |> 
  drop_na() |> 
  select(-time_hour)

merge_data_2017 = merge(flights_2017_clean, weather_2017_clean, by =c("origin", 
                                                                      "year", 
                                                                      "month",
                                                                      "day",
                                                                    "hour"))


merge_data_2013 = merge_data_2013 |> 
  select(-flight, -tailnum)

merge_data_2017 = merge_data_2017 |>
  select(-flight, -tailnum)

```


The final clean column names and descriptions:

* `origin`: NYC flight origin (LGA, JFK or EWR)
* `year`:  Year of departure
* `month`: Month of departure
* `day`: Date of departure
* `hour` and `minute`: Time of departure, recorded in 24-hour clock, breaking into hour and minutes
* `arr_delay`: Arrival delays in minutes
* `carrier`: Two letter flight carrier abbreviation
* `dest`: destination
* `air_time`: Flight duration in minutes
* `distance`: Distance between origin and destination in miles
* `wind_dir`: Wind direction in degree
* `wind_speed`: Wind speed in mph
* `wind_gust`: Wind gust speed in mph
* `precip`: Precipitation in inches
* `pressure`: Sea level pressure in millibars
* `visib`: Visibility in miles



## Exploratory Analysis

To provide an overall view of the distribution in arrival delay (in minutes), we created a histogram for the arrival delay in 2013, a barplot showing the arrival status categorized by three NYC departure airports (`LGA`, `JFK`, `EWR`) and a boxplot below displaying the top 6 average arrival delay grouped by destinations. Additionally, we then visualized the correlation between arrival delay and weather features (e.g.,`pressure`, `wind`, `visibility`, etc.) and some other external factors including `flight carriers`, `month` and specific datetime features. 

Note: All visualization parts were performed on data consisting of flight and weather information in 2013. Data in 2017 was reserved for testing purpose only. 

```{r load_library, warning = FALSE, message = FALSE}
library(tidyr)
library(tidyverse)
library(rvest)
library(dplyr)
library(cowplot)
library(gridExtra)
library(RColorBrewer)
library(plotly)
library(corrplot)
library(ggpubr)
library(viridisLite)
library(patchwork)

knitr::opts_chunk$set(
    echo = TRUE,
    warning = FALSE,
    fig.width = 8, 
  fig.height = 6,
  out.width = "90%"
)

theme_set(theme_minimal() + 
          theme(legend.position = "bottom",
                plot.title= element_text(hjust = 0.5))
          )

options(
  ggplot2.continuous.colour = "viridis",
  ggplot2.continuous.fill = "viridis"
)

scale_colour_discrete = scale_colour_viridis_d
scale_fill_discrete = scale_fill_viridis_d

## Load dataset
df_2013 = 
  read_csv("data/merge_data_2013.csv", show_col_types = FALSE)


# Convert month to factor with levels in ascending order and labels as month abbreviations
df_2013$month <- factor(df_2013$month, levels = 1:12, labels = month.abb[1:12])

df_2013 <- df_2013 %>%
  mutate(date = paste(month, day, sep = "_"))


average_delay_by_date <- df_2013 %>%
  # group_by(date) %>%
  group_by(month, date) %>%
  summarise(
    avg_arr_delay = mean(arr_delay, na.rm = TRUE),
    avg_precip = mean(precip, na.rm = TRUE),
    avg_wind_dir = mean(wind_dir, na.rm = TRUE),
    avg_wind_speed = mean(wind_speed, na.rm = TRUE),
    avg_wind_gust = mean(wind_gust, na.rm = TRUE),
    avg_pressure = mean(pressure, na.rm = TRUE),
    avg_visib = mean(visib, na.rm = TRUE)
    )


# Aggregate the data by month and calculate averages
df_2013_avg <- df_2013 %>%
  group_by(month) %>%
  summarize(
    avg_arr_delay = mean(arr_delay, na.rm = TRUE),
    avg_precip = mean(precip, na.rm = TRUE),
    avg_wind_dir = mean(wind_dir, na.rm = TRUE),
    avg_wind_speed = mean(wind_speed, na.rm = TRUE),
    avg_wind_gust = mean(wind_gust, na.rm = TRUE),
    avg_pressure = mean(pressure, na.rm = TRUE),
    avg_visib = mean(visib, na.rm = TRUE)
  )
```


### Arrival Delay Summary

```{r delay_summary, warning = FALSE, message = FALSE,fig.width=10}
dd_hist = 
  ggplot(data = df_2013, aes(x = arr_delay)) +
  geom_histogram(fill = 'skyblue', color = 'black')+
  labs(title = "Distribution of Arrival Delay",
       x = "Arrival Delay (minutes)",
       y = "Frequency")

origin_delay = df_2013 |> 
  mutate(
    arrival_type = ifelse(arr_delay <= 0, "on time", "delayed")) |> 
  ggplot(aes(x = origin, fill = arrival_type)) +
    geom_bar()+
    scale_fill_brewer(palette = "Set3") + 
    labs(title = "Arrival Status Depatured from NYC Airports",
         x = "Origin",
         y = "Counts")

# Calculate average arrival delay by destination
avg_arr_delay_by_dest <- df_2013 %>%
  group_by(dest) %>%
  summarize(avg_arr_delay_dest = mean(arr_delay, na.rm = TRUE)) %>%
  arrange(desc(avg_arr_delay_dest)) %>%
  head(6)

# Select only the top 6 destinations
top_destinations <- avg_arr_delay_by_dest$dest

# Filter the data for the top 6 destinations
filtered_data <- df_2013 %>% filter(dest %in% top_destinations)

# Create a boxplot for the top 6 destinations
top6_dest = 
  ggplot(filtered_data, aes(x = reorder(dest, -arr_delay), y = arr_delay, fill = dest)) +
  geom_boxplot() +
  scale_fill_brewer(palette = "Set3") +
  labs(title = "Arrival Delay for Top 6 Destinations",
       x = "Destination",
       y = "Arrival Delay (minutes)")


dd_hist | origin_delay / top6_dest

```


```{r}
df_2013 |> 
  summarise(min_delay = min(arr_delay),
            avg_delay = mean(arr_delay),
            median_delay = median(arr_delay),
            max_delay = max(arr_delay)
            ) |>
  knitr::kable(
    digits = 3,
    col.names = c("Min Delay Time (Minutes)","Average Delay Time (Minutes)", "Median Delay Time (Minutes)", 
                  "Max Delay Time (Minutes)")
    )
```


```{r overview statistical table}
df_2013 |> 
  mutate(
    arrival_type = ifelse(arr_delay <= 0, "on time", "delayed")) |> 
  group_by(arrival_type) |> 
  summarise(count = n(),
            average_delay_time = mean(arr_delay)) |> 
  knitr::kable(
    digits = 3,
    col.names = c("Arrival Type", "Count", "Average Delay Time (minutes)")
    )
```

The overall distribution of arrival delay was highly skewed to the right, with the average delay time 7.28 minutes. Based on the barplot of arrival delay grouped by origins, we found that EWR had the highest portion of delays in respect to its total flights, and JFK had the the lowest number of delays.Based on the boxplot displaying average time of arrival delay, we did not perceive any major differences regarding arrival delay time under the top 6 destinations.Approximately 58% of records (n=42,150) in the dataset were reported as on-time, and among them, the average time in earlier arrival was 15 minutes. On the contrary, in the delay group, the delay time was calculated to be 38 minutes on average. 


### Correlation Plot
```{r}
df_2013 = read.csv("data/merge_data_2013.csv") |> 
  mutate(month = factor(month, levels = 1:12, labels = month.abb[1:12]))
df_corr = df_2013 |> 
  select(-year,-flight, -day, -minute, -hour) |> 
  select_if( is.numeric)
corrplot(cor(df_corr), type="upper", order="hclust",
         col=brewer.pal(n=8, name="RdYlBu"))
```

The correlation heatmap plot depicted some highly correlated pair of features such as `wind_speed` and `precip`, and `air_time` and `distance`. Thus, in the modeling process, we will only keep one variable for each pair to avoid the multicollinearity issue. In this case, we removed features `wind speed` and `distance` when constructing the model.


### Arrival Delay & Weather Factors
#### Pressure

```{r pressure_date, warning = FALSE, message = FALSE}
# Average pressure against average arr_delay by month
pressure_delay_date = 
  ggplot(average_delay_by_date, aes(x = avg_pressure, y = avg_arr_delay)) +
  geom_point(size = 3, aes(x = avg_pressure, y = avg_arr_delay, color = month)) +
  geom_smooth(method = "lm", se = FALSE, color = "coral", size = 0.6) +
  stat_cor(method = "pearson", label.x = 1020, label.y = -30, size = 5) +  
  labs(title = "Scatter Plot of Average Pressure against Arrival Delay by Date",
       x = "Avg Pressure by Date",
       y = "Avg Arrival Delay (minutes)")

pressure_delay_date

```

From the plot we could observe a trend that as pressure increased, the arrival delay time decreased. Since the trend seemed not very clear, we then calculated the correlation coefficient between the two variables. The correlation coefficient was -0.26, suggesting a weak negative correlation between pressure (mmhg) and arrival delay in minutes. While the p-value was significantly small, it might be driven by a large sample size (n=72,734). 



#### Visibility

```{r visib}

df_2013 |> 
  group_by(visib) |> 
  mutate(count_delay = if_else(arr_delay>0 , 1, 0)) |> 
  summarise(avg_delay = mean(arr_delay),
            number_of_delays = sum(count_delay)
            ) |> 
  ggplot(aes(x=visib,
             y=avg_delay, color=number_of_delays))+
   geom_point() +
  labs(title = "Arrival Delay Vs. Visibility",
      x = "Visibility",
      y = "Average of Arrival Delay Time (Minutes)") 

```

Based on the graph, it appeared that the incremental of visibility did not linearly correlate with the reduction in the number of delays. In the highest visibility 10, we perceived the highest number of delays. There was no major difference in the number of delays in visibility ranging from 0 to 9. However, we observed that visibility negatively correlated with average delay time in minutes. 

#### Wind
```{r plot_weather_wind, warning = FALSE, message = FALSE,  fig.height=7}
# Average Wind Direction against average arr_delay by date
wind_dir_delay_date = 
  ggplot(average_delay_by_date, aes(x = avg_wind_dir, y = avg_arr_delay)) +
  geom_point(size = 3, aes(x = avg_wind_dir, y = avg_arr_delay, color = month)) +
  geom_smooth(method = "lm", se = FALSE, color = "coral", size = 0.6) +
  stat_cor(method = "pearson", size = 5) + 
  guides(color = 'none') +  
  labs(title = "Avg Wind Direction VS Arrival Delay",
       x = "Avg Daily Wind Direction",
       y = "Avg arr_delay")

# Average wind_gust against average arr_delay by date
wind_gust_delay_date = 
  ggplot(average_delay_by_date, aes(x = avg_wind_gust, y = avg_arr_delay)) +
  geom_point(size = 3, aes(x = avg_wind_dir, y = avg_arr_delay, color = month)) +
  geom_smooth(method = "lm", se = FALSE, color = "coral", size = 0.6) +
  stat_cor(method = "pearson", size = 5) + 
  labs(title = "Avg Wind Gust VS Arrival Delay",
       x = "Avg Daily Wind Gust",
       y = "Avg arr_delay")

# Average wind_speed against average arr_delay by date
wind_speed_delay_date = 
  ggplot(average_delay_by_date, aes(x = avg_wind_speed, y = avg_arr_delay)) +
  geom_point(size = 3, aes(x = avg_wind_speed, y = avg_arr_delay, color = month)) +
  geom_smooth(method = "lm", se = FALSE, color = "coral", size = 0.6) +
  stat_cor(method = "pearson", size = 5) + 
  guides(color = 'none') +  
  labs(title = "Avg Wind Speed VS Arrival Delay",
       x = "Avg Daily Wind Speed",
       y = "Avg arr_delay")


(wind_dir_delay_date | wind_speed_delay_date) / wind_gust_delay_date

```

The three scatterplots with fitted lines here illustrated the relationship between wind features and arrival delay. We also calculated the correlation coefficients to confirm the trends we observed. 
According to the first graph, arrival delay was negatively related to wind direction. The correlation coefficient of -0.19 suggested that the negative correlation between wind direction and arrival delay was weak. In the second graph, generally arrival delay increased as wind speed increased, and this positive correlation, though weak, was also shown by its correlation coefficient of 0.1. Through the third graph and the calculated correlation coefficient of 0.11, a weak positive correlation between arrival delay and wind gust could be observed and derived. It was worth noticing that the significantly small p-values calculated for the three correlation coefficients might still be caused by the large sample size (n=72,734) rather than a solid evidence of the correlations. 




### Arrival Delay & Carriers, Temporal Factors

```{r, fig.width=10}

carrier_plot = df_2013 |> 
  group_by(carrier) |> 
  mutate(delay_count = if_else(arr_delay >0 , 1, 0)) |> 
  summarise(delay_count = sum(delay_count)) |> 
  ggplot(aes(x = carrier,
            y = delay_count)) +
  geom_bar(stat = "identity",  fill = viridis(16)) +
  labs(title = "Number of Delays for Each Carrier", 
        x = "Carrier",
        y = "Total Number of Delays")


month_plot = df_2013 |> 
  group_by(month) |> 
  mutate(delay_count = if_else(arr_delay >0 , 1, 0)) |> 
  summarise(delay_count = sum(delay_count)) |> 
  ggplot(aes(x = month,
            y = delay_count, group=1)) +
  geom_point(color="#0aa2fa") +
  geom_line(color="#0aa2fa") +
  labs(title = "Number of Delays for Each Month", 
        x = "Month",
        y = "Total Number of Delays")

carrier_plot | month_plot

```

The carrier `EV` was shown to have the highest number of delays, followed by carrier `B6`. The number of delays for carrier `HA` was found to be the lowest across all carriers. Months such as March and April appeared to have the relatively high number of delays, and the lowest total number of delays occurred in September. 


```{r, message=FALSE, fig.width=9}
df_2013_weekday = read_csv("data/merge_data_2013.csv") |> 
  mutate(delay_count = if_else(arr_delay >0 , 1, 0),
         arrival_date = paste(year,"-",month,"-",day, sep = "")
         ) |>
  mutate(week_date = weekdays(date(arrival_date)),
          week_date = factor(week_date, 
                             levels = c("Monday", "Tuesday", "Wednesday", 
                                        "Thursday", "Friday", "Saturday",
                                        "Sunday"))
  )

weekday_plot = df_2013_weekday |> 
  group_by(week_date) |> 
  summarise(delay_count = sum(delay_count)) |> 
  ggplot(aes(x = week_date,
            y = delay_count, group=1)) +
    geom_point(color="#6f98d1") +
  geom_line(color="#6f98d1") +
  labs(title = "Number of Delays Vs. Weekday", 
        x = "Weekday",
        y = "Total Number of Delays")

weekday_plot_delay = df_2013_weekday |> 
  group_by(week_date) |> 
  summarise(avg_delay = mean(arr_delay)) |> 
  ggplot(aes(x = week_date,
            y = avg_delay, group=1)) +
    geom_point(color="#6f98d1") +
  geom_line(color="#6f98d1") +
  labs(title = "Average Time of Arrival Delays (minutes) Vs. Weekday", 
        x = "Weekday",
        y = "Average Time of Arrival Delays (minutes)")


weekday_plot /weekday_plot_delay
  
  
```

An increasing trend in the number of delays was perceived from Monday to Thursday, followed by a decline on Friday and Saturday. Based on the figures, Saturday had the lowest number of delays and average time of arrival delays. It was also noted that Thursday was found to have the highest number of delays and average time of arrival delays in minutes. 

```{r, fig.width=9}
part_day_plot = df_2013 |> 
  mutate(hour = cut(hour, 
                    breaks = c(-Inf, 5, 12, 17, 21, Inf),
                    labels = c("Late Night", "Morning", "Afternoon", "Evening", "Late Night"),
                    include.lowest = TRUE
                    )) |>
  group_by(hour) |>
  mutate(delay_count = if_else(arr_delay >0 , 1, 0)) |> 
  summarise(delay_count = sum(delay_count),
            avg_delay = mean(arr_delay)) |> 
  ggplot(aes(x = hour,
            y = delay_count, group=1)) +
  geom_point(color="#eb8f2d") +
  geom_line(color="#eb8f2d") +
  labs(title = "Number of Delays Vs. Daytime",
       x = "Daytime",
       y = "Total Number of Delays")


part_day_plot_delay = df_2013 |> 
  mutate(hour = cut(hour, 
                    breaks = c(-Inf, 5, 12, 17, 21, Inf),
                    labels = c("Late Night", "Morning", "Afternoon", "Evening", "Late Night"),
                    include.lowest = TRUE
                    )) |>
  group_by(hour) |>
  mutate(delay_count = if_else(arr_delay >0 , 1, 0)) |> 
   summarise(avg_delay = mean(arr_delay)) |> 
  ggplot(aes(x = hour,
            y = avg_delay, group=1)) +
  geom_point(color="#eb8f2d") +
  geom_line(color="#eb8f2d") +
  labs(title = "Average Time of Arrival Delay (minutes) Vs. Daytime",
       x = "Daytime",
       y = "Average Arrival Delays (minutes)")

part_day_plot / part_day_plot_delay
  
```

Furthermore, we further visualized the changes of number of delays and daytime divided into `Late Night`, `Morning`, `Afternoon`, `Evening`. The number of arrival delays was increased from time during laye night to afternoon, dropping during the evening. We found a largest number of delays during afternoon, with the average time delay of 9 minutes. 

### Insights from Visualization
Based on above visualizations, we anticipated weather factors such as pressure,visibility and wind conditions can be beneficial in predicting flight delay task. Additionally, other features such as information about flight carrier and temporal characteristics are believed to helpful in such prediction task as well. Thus, in the subsequent sections, we would further quantity their correlations by leveraging multiple linear regression.

## Additional Analysis

### Multiple Linear Regression

During the exploratory data analysis (EDA) phase, we identified highly correlated variable pairs. The presence of such pairs indicates a potential for one variable to essentially substitute another in the model, leading to multicollinearity issues. In response to this correlation finding, we assessed each variable within these pairs based on their unique values. Our criterion was to retain the variable with a higher number of unique values, signifying a broader range of information.

As a result, we opted to exclude `wind_speed` and `distance` from the analysis. This decision aims to improve the model's stability and interpretability by addressing multicollinearity issues.

|highly correlated variables  | wind_gust | wind_speed | air_time | distance |
|---|-----------|------------|----------|----------|
|unique cnt|33 |30|462|177|

#### Stepwise Regression

In the linear model context, we've taken crucial preprocessing steps, including standardizing variables, applying one-hot encoding for categorical variables, and removing highly correlated pairs. These actions were taken to improve the model's stability and interpretability.

With increased confidence in the processed dataset, we proceeded with stepwise regression. This method systematically refines the model by iteratively selecting or removing variables based on criteria like the Akaike Information Criterion (AIC). The aim was to precisely capture relationships between flight delay time and various flight information parameters.

By employing stepwise regression, our goal was to ensure that the final model incorporates the most relevant and significant predictors. This approach provides a streamlined and statistically robust representation of the intricate relationships within the data.

The plot below unveils insights into the top 20 estimated values of the predictors. Strikingly, these top 20 predictors exclusively consist of destination dummy variables. Notably, only Daniel K. Inouye International Airport (HNL) stands out with negative estimates, implying that flights to HNL are less likely to experience delays—a delightful discovery. Perhaps both visitors and flight crews couldn't wait to join the beach!

For a more comprehensive overview of the Multiple Linear Regression (MLR) model, please refer to the [Analysis-MLR page](https://yuxinyin0906.github.io/p8105_final_project/MLR.html).In brief, the Adjusted 
$R^2$ stands at a mere 0.110, suggesting that the linear model falls short in adequately addressing delay-related issues.

![Estimates of MLR](image/mlr_estimate.png){width=80%}

### Machine Learning Methods

After the initial acquisition of insights, we wanted to further analyze the data and build predictive models to predict whether a flight is delayed or not based on highly correlated variables. In order to reach the goal, we chose two models, lasso regression and random forest as they bring different strengths to the table. 

#### Lasso Regression

Using the `logistic_reg` function with hyperparameter tuning for the penalty and mixture parameters, we defined the logistic regression model. The glmnet engine was set up to be used by the model. Feature preprocessing included managing data variables, eliminating unnecessary columns, creating dummy variables for categorical predictors, removing zero-variance predictors, and normalizing all predictors.  <br>
<br>
We then conducted a grid search was conducted over a range of penalty values and saved the results for further analysis. The best-performing hyperparameters were selected based on the highest mean Area under the ROC Curve (AUC). By highlighting the penalty value corresponding to the peak AUC, we visualized the AUC-penalty relationship.

![Finding the Best Hyperparameter](image/best_hyperparameter.png){width=80%}

<br>
After that, we trained the logistic regression model with the optimal hyperparameters on the entire dataset, and examined coefficients to identify the top 20 predictors. 

#### Random Forest Regression

The first step we took was detecting and storing the number of available CPU cores in the variable `cores`. This information is essential for optimizing the performance of parallelized computations. <br>
<br>
We then specified hyperparameters for the random forest model and adjusted it during the model training phase. We employed a preprocessing recipe to handle date variables, remove redundant columns, and create dummy variables (`df_other`) for categorical predictors on the training data. <br>
<br>
We used the area under the ROC curve metric to assess the model once hyperparameter tuning was completed on a grid of values for `mtry` and `min_n` <br>
<br>
One of the most important steps was comparing the performance of the random forest model with the previously developed logistic regression model. We visualized sensitivity and specificity metrics on a ROC curve for both models using the training set. 

![Comparison between Two Models](image/rf.png){width=80%}

<br>
By looking at the graph, we could clearly find that the random forest model performs better than the logistic regression model.

Finally, we trained the final random forest model with the optimal hyperparameters and evaluated its performance. As we did for logistic regression model, we also visualized the top 20 predictors to highlight the feature importance of the model.

## Discussion

The objective of the study was to assess the correlation between flight delay duration and various flight-related factors using multiple linear regression. Additionally, the research aimed to forecast flight delays by employing a machine learning framework.

### Variable Selection

Lasso Regression and Random Forest models were applied to analyze the 2013 New York City flight dataset. The accompanying plots shed light on the 'variable importance.' On the right side, the top 20 absolute value estimates of Lasso Regression reveal a predominant presence of destination dummy variables. In contrast, the Random Forest plot on the left tells a distinct tale. Continuous variables such as airtime, pressure, distance, etc., play a more significant role in predicting delays compared to other categorical variables.

![Caption for Image 1](image/lr_varplot.png){width=45%}\hfill
![Caption for Image 2](image/rf_varplot.png){width=45%}

The difference in 'variable importance' between Lasso Regression and Random Forest models can be attributed to the inherent characteristics and mechanisms of each algorithm:

Lasso Regression (LR):

* Preference for Categorical Variables: LR performs variable selection by shrinking the coefficients of less influential variables to exactly zero. It tends to favor sparse solutions and is particularly effective when dealing with high-dimensional data.

Random Forest (RF):

* Ensemble of Decision Trees: RF builds an ensemble of decision trees, each trained on a random subset of the data and features. This ensemble approach allows the model to capture complex interactions and non-linear relationships.

* Handling Mixed Data Types: Unlike LR, RF is well-suited for handling both categorical and continuous variables simultaneously. It can effectively evaluate the importance of features with various data types.

* Robust to Collinearity: RF is less sensitive to multicollinearity, which means it can handle correlated features without necessarily discarding one of them. This can be advantageous when dealing with a mix of categorical and continuous variables.

### Model Performance

2017 New York City flight dataset was used for model evaluation. Notably, the Random Forest model exhibited superior performance across a majority of evaluation metrics. 


|                    | Sensitivity | Specificity | Pos Pred Value | Neg Pred Value | Precision | Recall | F1    | Prevalence | Detection Rate | Detection Prevalence | Balanced Accuracy |
|--------------------|-------------|-------------|-----------------|-----------------|-----------|--------|-------|------------|-----------------|-----------------------|---------------------|
| Logistic           | 0.880       | 0.470       | 0.759           | 0.673           | 0.759     | 0.880  | 0.815 | 0.655      | 0.577           | 0.760                 | 0.675               |
| RF     | 0.971       | 0.710       | 0.864           | 0.928           | 0.864     | 0.971  | 0.915 | 0.655      | 0.637           | 0.736                 | 0.841               |
 

Random Forest might outperform Lasso Regression in predicting flight delays for several reasons:


* Handling Non-linearity: Flight delay prediction is likely influenced by a combination of linear and non-linear factors. RF, being a non-linear model, can handle complex relationships and interactions more effectively than LR, which assumes a linear relationship between the predictors and the response variable.

* Feature Importance: RF inherently assesses the importance of different features in making predictions. This feature importance can be beneficial in scenarios where certain variables play a more crucial role in determining outcomes, and RF can adapt to this by giving more weight to those features. LR, while performing feature selection by shrinking some coefficients to zero, might not capture non-linear relationships or account for intricate feature interactions as comprehensively.

* Robustness to Outliers: RF is generally more robust to outliers and noisy data compared to LR.