More info on dataset

df_status(master_data)

##               variable q_zeros p_zeros   q_na  p_na q_inf p_inf      type
## 1              fl_date       0    0.00      0  0.00     0     0 character
## 2    op_unique_carrier       0    0.00      0  0.00     0     0 character
## 3             tail_num       0    0.00    397  0.24     0     0 character
## 4    op_carrier_fl_num       0    0.00      0  0.00     0     0   numeric
## 5               origin       0    0.00      0  0.00     0     0 character
## 6                 dest       0    0.00      0  0.00     0     0 character
## 7             dep_time       0    0.00   3597  2.16     0     0   numeric
## 8            dep_delay    7494    4.51   3769  2.27     0     0   numeric
## 9             taxi_out       0    0.00   3662  2.20     0     0   numeric
## 10             taxi_in       0    0.00   3674  2.21     0     0   numeric
## 11            arr_time       0    0.00   3674  2.21     0     0   numeric
## 12           arr_delay    2711    1.63   3943  2.37     0     0   numeric
## 13           cancelled  162527   97.79      0  0.00     0     0   numeric
## 14 actual_elapsed_time       0    0.00   3849  2.32     0     0   numeric
## 15            air_time       0    0.00   3849  2.32     0     0   numeric
## 16            distance       0    0.00      0  0.00     0     0   numeric
## 17       carrier_delay   15073    9.07 140602 84.60     0     0   numeric
## 18       weather_delay   23853   14.35 140602 84.60     0     0   numeric
## 19           nas_delay   10806    6.50 140602 84.60     0     0   numeric
## 20      security_delay   25539   15.37 140602 84.60     0     0   numeric
## 21 late_aircraft_delay   12707    7.65 140602 84.60     0     0   numeric
##    unique
## 1      89
## 2      17
## 3    4712
## 4    3818
## 5     195
## 6       3
## 7    1297
## 8     698
## 9     152
## 10    109
## 11   1406
## 12    734
## 13      2
## 14    434
## 15    413
## 16    287
## 17    558
## 18    290
## 19    321
## 20     32
## 21    393

Missing values

plot_missing(master_data)

Converting the date in year-month-date format

master_data <- master_data%>%mutate(fl_date=ymd(fl_date))

Generating day of week column

master_data <- master_data%>%mutate(dow=wday(fl_date,label=T))
ggplot(master_data,aes(x=master_data$dow))+geom_histogram(stat='count')

## Warning: Ignoring unknown parameters: binwidth, bins, pad

Filtering on ORD and ATL destination airports

md_pred <- master_data%>%filter(dest %in% c("ORD", "ATL"))

Removing negative arrival and departure delay

md_pred <- md_pred%>%filter(arr_delay> 0 & dep_delay>0)
dim(md_pred)

## [1] 31820    22

Frequency of origin airports

origin <- md_pred%>%group_by(origin)%>%summarise(cnt=n())%>%arrange(desc(cnt))

kable(head(origin,194),'html') %>%
 kable_styling() %>%
 scroll_box(width = "400px", height = "700px")

Binning origin airports using case statement

origin <- origin%>%filter(cnt>=10)%>%mutate(origin_freq=case_when(cnt >=175 & cnt<300 ~ 'Low' ,
                                      cnt>=300 & cnt<500 ~ 'Moderate',
                                      cnt>=500 ~ 'High',
                                      TRUE ~ 'Very Low'))

kable(head(origin,25),'html') %>%
 kable_styling() %>%
 scroll_box(width = "400px", height = "700px")

Joining the two datasets to add a new origin frequency column

md_pred <- left_join(md_pred,origin,by='origin')
md_pred <- md_pred%>%filter(!is.na(origin_freq))

Adding month column

md_pred <- md_pred%>%mutate(Month=lubridate::month(ymd(fl_date),label=T))

Adding levels to categorical columns

md_pred <- md_pred %>% 
  mutate_at(vars(op_unique_carrier,dest,origin_freq, Month,dow), funs(as.factor))

## Warning: funs() is soft deprecated as of dplyr 0.8.0
## please use list() instead
## 
## # Before:
## funs(name = f(.)
## 
## # After: 
## list(name = ~f(.))
## This warning is displayed once per session.

md_pred$Month <- factor(md_pred$Month, levels = c("Jan","Feb","Mar"), labels =c(1,2,3),
                        ordered = F )
md_pred$dow <- factor(md_pred$dow, levels = c("Mon","Tue","Wed","Thu","Fri","Sat","Sun"), 
                        labels =c(1,2,3,4,5,6,7),ordered = F )
md_pred$origin_freq <- factor(md_pred$origin_freq, levels = c("Very Low","Low","Moderate","High"), 
                                labels =c(1,2,3,4) )

md_pred$dest <- factor(md_pred$dest, levels = c("ORD","ATL"), 
                                labels =c(1,2) )

md_pred$op_unique_carrier <- factor(md_pred$op_unique_carrier, 
                                    levels = c("F9","EV","MQ","NK","OH","OO","VX","WN","YV","YX","AA", "AS", "DL", "9E", "B6", "UA"), 
                                    labels =c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16) )

df_status(md_pred)

##               variable q_zeros p_zeros  q_na  p_na q_inf p_inf      type
## 1              fl_date       0    0.00     0  0.00     0     0      Date
## 2    op_unique_carrier       0    0.00     0  0.00     0     0    factor
## 3             tail_num       0    0.00     0  0.00     0     0 character
## 4    op_carrier_fl_num       0    0.00     0  0.00     0     0   numeric
## 5               origin       0    0.00     0  0.00     0     0 character
## 6                 dest       0    0.00     0  0.00     0     0    factor
## 7             dep_time       0    0.00     0  0.00     0     0   numeric
## 8            dep_delay       0    0.00     0  0.00     0     0   numeric
## 9             taxi_out       0    0.00     0  0.00     0     0   numeric
## 10             taxi_in       0    0.00     0  0.00     0     0   numeric
## 11            arr_time       0    0.00     0  0.00     0     0   numeric
## 12           arr_delay       0    0.00     0  0.00     0     0   numeric
## 13           cancelled   31766  100.00     0  0.00     0     0   numeric
## 14 actual_elapsed_time       0    0.00     0  0.00     0     0   numeric
## 15            air_time       0    0.00     0  0.00     0     0   numeric
## 16            distance       0    0.00     0  0.00     0     0   numeric
## 17       carrier_delay   10871   34.22 10411 32.77     0     0   numeric
## 18       weather_delay   19618   61.76 10411 32.77     0     0   numeric
## 19           nas_delay   10773   33.91 10411 32.77     0     0   numeric
## 20      security_delay   21303   67.06 10411 32.77     0     0   numeric
## 21 late_aircraft_delay    8496   26.75 10411 32.77     0     0   numeric
## 22                 dow       0    0.00     0  0.00     0     0    factor
## 23                 cnt       0    0.00     0  0.00     0     0   integer
## 24         origin_freq       0    0.00     0  0.00     0     0    factor
## 25               Month       0    0.00     0  0.00     0     0    factor
##    unique
## 1      89
## 2      16
## 3    4230
## 4    3103
## 5     182
## 6       2
## 7    1263
## 8     660
## 9     147
## 10    100
## 11   1331
## 12    668
## 13      1
## 14    331
## 15    308
## 16    275
## 17    557
## 18    289
## 19    321
## 20     32
## 21    393
## 22      7
## 23    128
## 24      4
## 25      3

Selecting final columns

md_pred <-  md_pred%>%select(2,6,8,12,14,16:22,24,25)

We see there are 5 factor columns and 9 numerical columns mentioned with extensive information about both the variable types. Many variables are right skewed along with arrival delay. The skewness needs to be addressed.

Data preparation

recipe_obj <- recipe(arr_delay ~ ., data = md_pred) %>%
    step_zv(all_predictors()) %>%
    step_YeoJohnson(nas_delay,late_aircraft_delay,security_delay,weather_delay,carrier_delay,arr_delay)%>%
    prep()

md_pred_final dataset

md_pred_final <- bake(recipe_obj, new_data = md_pred)
skimr::skim(md_pred_final)

## Skim summary statistics
##  n obs: 31766 
##  n variables: 14 
## 
## ── Variable type:factor ───────────────────────────────────────────────────
##           variable missing complete     n n_unique
##               dest       0    31766 31766        2
##                dow       0    31766 31766        7
##              Month       0    31766 31766        3
##  op_unique_carrier       0    31766 31766       16
##        origin_freq       0    31766 31766        4
##                            top_counts ordered
##             2: 16957, 1: 14809, NA: 0   FALSE
##    1: 5878, 3: 5038, 4: 4810, 2: 4567   FALSE
##     1: 12062, 3: 9934, 2: 9770, NA: 0   FALSE
##  13: 9467, 6: 5874, 11: 3297, 3: 2943   FALSE
##    1: 9499, 4: 8874, 3: 7011, 2: 6382   FALSE
## 
## ── Variable type:numeric ──────────────────────────────────────────────────
##             variable missing complete     n    mean      sd    p0    p25
##  actual_elapsed_time       0    31766 31766 128.61   52.69  37     90   
##            arr_delay       0    31766 31766   3.27    1.22   0.69   2.46
##        carrier_delay   10411    21355 31766   0.85    0.92   0      0   
##            dep_delay       0    31766 31766  58.46   95.01   1     13   
##             distance       0    31766 31766 660.87  463.9   67    334   
##  late_aircraft_delay   10411    21355 31766   2.06    1.83   0      0   
##            nas_delay   10411    21355 31766   0.87    0.94   0      0   
##       security_delay   10411    21355 31766   0.059   2.24   0      0   
##        weather_delay   10411    21355 31766   0.022   0.073  0      0   
##     p50    p75    p100     hist
##  117    156     539    ▆▇▃▁▁▁▁▁
##    3.29   4.13    7.1  ▂▅▇▇▆▃▁▁
##    0      1.78    2.58 ▇▁▁▁▂▂▂▁
##   29     66    1439    ▇▁▁▁▁▁▁▁
##  584    763    4502    ▇▃▁▁▁▁▁▁
##    2.64   3.63    6.78 ▇▁▁▅▃▂▁▁
##    0      1.83    2.71 ▇▁▁▁▂▂▂▁
##    0      0     198    ▇▁▁▁▁▁▁▁
##    0      0       0.27 ▇▁▁▁▁▁▁▁

Most frequent carriers on monthly basis

ggplot(data = md_pred_final) +
  aes(x = op_unique_carrier) +
  geom_bar(color = "black", fill = "lightblue") + facet_grid(Month~.)

Average delay time per month for unique carrier

a <- md_pred_final%>%group_by(Month,op_unique_carrier)%>%summarise(avg=mean(arr_delay))
md_pred_final%>%group_by(Month)%>%summarise(avg=mean(arr_delay))

## # A tibble: 3 x 2
##   Month   avg
##   <fct> <dbl>
## 1 1      3.40
## 2 2      3.35
## 3 3      3.03

p <- ggplot(data = a,aes(x = op_unique_carrier,y=avg) )+
  geom_bar(stat = 'identity',fill='#807dba') + facet_grid(Month~.)

p + geom_hline(aes(yintercept = avg), data.frame(Month = c(1, 2, 3), avg = c(3.40,3.35, 3.03)))

Carrier frequency on a daily basis

ggplot(data = md_pred_final) +
  aes(x = op_unique_carrier) +
  geom_bar(color = "black", fill = "lightblue") + facet_grid(dow~.)

Average delay time per day for 16 unique carriers

b <- md_pred_final%>%group_by(dow,op_unique_carrier)%>%summarise(avg=mean(arr_delay))
md_pred_final%>%group_by(dow)%>%summarise(avg=mean(arr_delay))

## # A tibble: 7 x 2
##   dow     avg
##   <fct> <dbl>
## 1 1      3.30
## 2 2      3.24
## 3 3      3.47
## 4 4      3.14
## 5 5      3.25
## 6 6      3.07
## 7 7      3.31

p <- ggplot(data = b,aes(x = op_unique_carrier,y=avg) )+
  geom_bar(stat = 'identity',fill='#807dba') + facet_grid(dow~.)

p + geom_hline(aes(yintercept = avg), data.frame(dow = c(1, 2, 3,4,5,6,7), avg = c(3.30, 3.24, 3.47,3.14, 3.25,3.07,3.31)))

Imputing missing values using Random Forest

class(md_pred)

## [1] "tbl_df"     "tbl"        "data.frame"

#md_df <- data.frame(md_pred_final)
#md_pred_imp<- missForest(md_df,maxiter = 2)

#md_pred_impute <- md_pred_imp$ximp
#kable(head(md_pred_impute,50),'html') %>%
 #kable_styling() %>%
 #scroll_box(width = "800px", height = "700px")

Data summarization of imputed dataset

#skimr::skim(md_pred_impute)

Significant variables can be seen from the multiple linear regression model. All the delay variables are highly significant along with wednesday

plot(reg)

QQ plot- Residuals are left skewed. From other 3 graphs we see there are some outliers.

Finding outliers in md_pred_final dataset

 bp<- ggplot(data = md_pred_final) +
  aes(x = arr_delay, y = arr_delay) +
  geom_boxplot(fill = '#9e9ac8') +
  labs(title = 'Box Plot of Arrival Delay') +
  theme_linedraw()

histo <- ggplot(data = md_pred_final) +
  aes(x = arr_delay) +
  geom_histogram(bins = 12, fill = '#9e9ac8') +
  labs(title = 'Histogram of Arrival Delay') +
  theme_linedraw()

plot_grid(bp,histo)

## Warning: Continuous x aesthetic -- did you forget aes(group=...)?

Removing the outliers

 md_pred_final <- md_pred_final%>%filter(arr_delay<6.9)

Splitting the data

train <- md_pred_final[sample(1:25411),]
colnames(train)

##  [1] "op_unique_carrier"   "dest"                "dep_delay"          
##  [4] "arr_delay"           "actual_elapsed_time" "distance"           
##  [7] "carrier_delay"       "weather_delay"       "nas_delay"          
## [10] "security_delay"      "late_aircraft_delay" "dow"                
## [13] "origin_freq"         "Month"

test <- md_pred_final[sample(25411:nrow(md_pred_final)),]
colnames(test)

##  [1] "op_unique_carrier"   "dest"                "dep_delay"          
##  [4] "arr_delay"           "actual_elapsed_time" "distance"           
##  [7] "carrier_delay"       "weather_delay"       "nas_delay"          
## [10] "security_delay"      "late_aircraft_delay" "dow"                
## [13] "origin_freq"         "Month"

write.csv(train,file="train.csv")
write.csv(test,file="test.csv")

Using H2O AutoML

h2o.init()

##  Connection successful!
## 
## R is connected to the H2O cluster: 
##     H2O cluster uptime:         19 hours 59 minutes 
##     H2O cluster timezone:       America/New_York 
##     H2O data parsing timezone:  UTC 
##     H2O cluster version:        3.24.0.2 
##     H2O cluster version age:    18 days  
##     H2O cluster name:           H2O_started_from_R_saurabhyelne_lbz280 
##     H2O cluster total nodes:    1 
##     H2O cluster total memory:   1.96 GB 
##     H2O cluster total cores:    4 
##     H2O cluster allowed cores:  4 
##     H2O cluster healthy:        TRUE 
##     H2O Connection ip:          localhost 
##     H2O Connection port:        54321 
##     H2O Connection proxy:       NA 
##     H2O Internal Security:      FALSE 
##     H2O API Extensions:         Amazon S3, XGBoost, Algos, AutoML, Core V3, Core V4 
##     R Version:                  R version 3.4.2 (2017-09-28)

train <- h2o.importFile("train.csv")

## 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |========                                                         |  12%
  |                                                                       
  |=================================================================| 100%

test <- h2o.importFile("test.csv")

## 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |=================================================================| 100%

split_h2o <- h2o.splitFrame(as.h2o(train), ratios = c(0.85))

train_h2o <- split_h2o[[1]]
valid_h2o <- split_h2o[[2]]
test_h2o  <- as.h2o(test)

Model generation

y <- 'arr_delay'
x <- setdiff(names(train),y)

#automl_models_h2o <- h2o.automl(
 # x = x,
  #y = y,
#  training_frame = train_h2o,
 # validation_frame = valid_h2o,
 # leaderboard_frame = test_h2o,
  #max_runtime_secs = 300,
  #nfolds = 5
#)

Best performing models

#automl_models_h2o@leaderboard

#print(paste0("The best performing model is ",automl_models_h2o@leader@model_id))

Description about the best model - GBM_1_AutoML_20190504_073500

#automl_models_h2o@leader

Saving the top performing model to working directory

#h2o.getModel("GBM_1_AutoML_20190504_073500") %>%
#  h2o.saveModel(path = "04_Modeling/h2o_models/")

#h2o_gbm <- h2o.loadModel("04_Modeling/h2o_models/GBM_1_AutoML_20190504_073500")

Prediction by Gradient Boosting Machine model

#predictions <- h2o.predict(h2o_gbm, newdata = as.h2o(test))

#predictions_tbl <- predictions %>% as.tibble()
#predictions_tbl

Performance of the top model

#performance_h2o_gbm <- h2o.performance(h2o_gbm, newdata = as.h2o(test))
#typeof(performance_h2o_gbm)
#performance_h2o_gbm %>% slotNames()
#performance_h2o@metrics
#train_rmse_gbm  <- h2o.rmse(gbm_h2o, train = TRUE)
#test_rmse_gbm   <- h2o.rmse(gbm_h2o, valid = TRUE)
#hold_rmse_gbm   <- h2o.rmse(object = performance_h2o_gbm)
#print(paste0("GBM rmse TRAIN = ", train_rmse_gbm, ", rmse valid = ", test_rmse_gbm, ", rmse HOLDOUT = ",
#             hold_rmse_gbm))

#h2o.varimp_plot(gbm_h2o, num_of_features = NULL)

[1] “GBM rmse TRAIN = 8.34806208521452, rmse valid = 17.6673196637401, rmse HOLDOUT = 6.83885108336663”

                                             **Thank You**