# training = 2022-2023, test = 2024
# path names
flight.folder <- "~/Downloads/Flight Punctuality Data UK, 2022-2024"
weather.folder <- "~/Downloads/Weather data"
set.seed(7406)
# 1. read CAA flight files ----------------------------------------------------------------------------
flight.files <- list.files(flight.folder, pattern = "\\.csv$", full.names = TRUE, recursive = TRUE)
flight.list <- list()
for (i in 1:length(flight.files)){
flight.list[[i]] <- read.csv(file = flight.files[i], sep = ",", header=T);
}
flights <- do.call(rbind, flight.list)
print(dim(flights)) #115754 rows
## [1] 115754 25
# 2. clean flight data ---------------------------------------------------------------------------------
flights <- flights[
#keep three target airport
flights$reporting_airport %in% c("HEATHROW", "GATWICK", "MANCHESTER") &
#keep scheduled flights only, removed chartered flights
flights$scheduled_charter == "S" &
#remove rows with zero match flights
!is.na(flights$number_flights_matched) & flights$number_flights_matched > 0,
]
print(table(flights$reporting_airport))
##
## GATWICK HEATHROW MANCHESTER
## 9568 11225 10462
# 3. read weather files ---------------------------------------------------------------------------------
weather.files <- list.files(weather.folder, pattern = "\\.csv$", full.names = TRUE, recursive = TRUE)
print(length(weather.files)) #4 files
## [1] 4
weather.list <- list()
for (i in 1:length(weather.files)){
weather.list[[i]] <- read.csv(file = weather.files[i], sep = ",", header=T);
}
weather <- do.call(rbind, weather.list)
print(dim(weather)) #3288 rows
## [1] 3288 33
# 4. standardize airport names to match CAA naming------------------------------------------------------
weather$airport <- trimws(weather$name)
weather$airport[weather$airport == "London Heathrow Airport, United Kingdom"] <- "HEATHROW"
weather$airport[weather$airport == "London Gatwick Airport, United Kingdom"] <- "GATWICK"
weather$airport[weather$airport == "Manchester Airport, United Kingdom"] <- "MANCHESTER"
weather <- weather[weather$airport %in% c("HEATHROW", "GATWICK", "MANCHESTER"), ]
print(table(weather$airport))
##
## GATWICK HEATHROW MANCHESTER
## 1096 1096 1096
# 5. create a reporting period matching the CAA data format---------------------------------------------
weather$datetime <- as.Date(weather$datetime,format = "%Y-%m-%d")
weather$reporting_period <- as.numeric(format(weather$datetime, "%Y%m"))
# 6. feature engineering weather variables--------------------------------------------------------------
# calculates monthly averages (airport-month-weather)
monthly.mean <- aggregate(
cbind(temp, humidity, precip, snow, windgust, windspeed,
sealevelpressure, cloudcover, visibility) ~ airport + reporting_period,
data = weather, FUN = mean)
# calculates monthly standard deviations
# addresses TA's suggestion..new features will help keep info within month weather variation
# selected these variables because their monthly average could hide anomalies; day to day variation within a month could affect operations
# EX. precip - month can have few days of heavy rain; windspeed - month could have calm weather for most of the month but a coupole of high wind days which affect operations
monthly.variability <- aggregate(
cbind(precip, snow, windspeed, visibility) ~ airport + reporting_period,
data = weather, FUN = sd)
#rename variables to std dev
names(monthly.variability)[3:6] <- c("precip_sd", "snow_sd", "windspeed_sd", "visibility_sd")
# calculates max windspeed during each airport month
# addresses TA's concern that monthly averages can hide severe wind days; strong winds could potentially affect operations
monthly.maxwind <- aggregate(
windspeed ~ airport + reporting_period,
data = weather, FUN = max)
names(monthly.maxwind)[3] <- "max_windspeed"
#combine all monthly weather features
weather.monthly <- merge(monthly.mean, monthly.variability,
by = c("airport", "reporting_period"))
weather.monthly <- merge(weather.monthly, monthly.maxwind,
by = c("airport", "reporting_period"))
print(dim(weather.monthly)) #108 rows, 16 cols
## [1] 108 16
print(head(weather.monthly))
## airport reporting_period temp humidity precip snow windgust windspeed
## 1 GATWICK 202201 40.28065 87.38710 0.03338710 0 22.51613 8.935484
## 2 GATWICK 202202 45.24286 80.10357 0.09128571 0 28.11786 15.496429
## 3 GATWICK 202203 46.03871 79.05484 0.04970968 0 21.18710 11.922581
## 4 GATWICK 202204 48.82667 72.83000 0.02213333 0 22.14667 12.553333
## 5 GATWICK 202205 55.88065 77.26452 0.06735484 0 20.04839 10.858065
## 6 GATWICK 202206 60.04000 74.29000 0.02836667 0 21.31667 11.270000
## sealevelpressure cloudcover visibility precip_sd snow_sd windspeed_sd
## 1 1025.848 60.05806 10.47419 0.09471208 0 4.134211
## 2 1016.196 65.84286 15.46786 0.15053063 0 5.817183
## 3 1021.823 58.57097 11.51290 0.11369175 0 3.620931
## 4 1015.723 57.15000 17.64667 0.05636013 0 4.636120
## 5 1018.445 65.20323 17.01290 0.14251048 0 3.168993
## 6 1015.913 54.00000 16.48000 0.06356885 0 2.798910
## visibility_sd max_windspeed
## 1 5.449157 17.5
## 2 5.227354 34.9
## 3 5.357098 21.3
## 4 3.839516 26.5
## 5 5.356351 18.0
## 6 4.215930 17.0
# 7. merge flights with monthly weather by airport + month----------------------------------------------------
flights.merged <- merge(flights, weather.monthly,
by.x = c("reporting_airport", "reporting_period"),
by.y = c("airport", "reporting_period"))
print(dim(flights.merged)) #31255x39
## [1] 31255 39
# create airport indicator variables and identifier
heathrow = I(flights.merged$reporting_airport == "HEATHROW")
gatwick = I(flights.merged$reporting_airport == "GATWICK")
#manchester is the reference airport
# creates airport month identifier
airport_month = interaction(flights.merged$reporting_airport,flights.merged$reporting_period,drop = TRUE)
# combine variables
flightdata = data.frame(heathrow,gatwick,airport_month,flights.merged)
#converts airport indiciators to 0 and 1
flightdata$heathrow = as.numeric(flightdata$heathrow)
flightdata$gatwick = as.numeric(flightdata$gatwick)
# 8. splits data before creating the response variable------------------------------------------------------
# split before calculating the median
# train = 2022-2023, test = 2024 (time based split)
flight.train = flightdata[flightdata$reporting_period <= 202312,]
flight.test = flightdata[flightdata$reporting_period >= 202401,]
# calculates median average delay threshold from training data (2022-2023)
# addresses TA'as concerns about data leakage because 2024 values aren't included in the threshold
delay.threshold = median(flight.train$average_delay_mins,na.rm = TRUE)
# binary response: above training median = 1, at or below training median = 0
flight.train$delayed = as.numeric(flight.train$average_delay_mins >delay.threshold)
#same threshold applied to 2024
flight.test$delayed = as.numeric(flight.test$average_delay_mins >delay.threshold)
# define predictors for models
continuous.predictors = c(
"number_flights_matched",
"temp",
"humidity",
"precip",
"snow",
"windspeed",
"sealevelpressure",
"cloudcover",
"visibility",
"precip_sd",
"snow_sd",
"windspeed_sd",
"visibility_sd",
"max_windspeed"
)
dummy.predictors = c("heathrow","gatwick")
#combines predictor lists
predictor.names = c(continuous.predictors,dummy.predictors)
#final list
required.model = c("delayed","reporting_period","airport_month",predictor.names)
print(nrow(flight.train)) #training observations = 20396
## [1] 20396
print(nrow(flight.test)) #test observations - 10859
## [1] 10859
print(table(flight.train$delayed)) #0 = 10199, 1=10197
##
## 0 1
## 10199 10197
print(table(flight.test$delayed)) #0 = 6458, 1 = 4401
##
## 0 1
## 6458 4401
# EDA ---------------------------------------------------------------------------------------------------------------------------
#boxplots
par(mfrow = c(3, 3), mar = c(4, 4, 3, 1))
boxplot(number_flights_matched ~ delayed, data = flight.train, main = "Flights Matched", xlab = "Delay Class", ylab = "Number of Flights")
boxplot(temp ~ delayed, data = flight.train, main = "Temperature", xlab = "Delay Class", ylab = "Temperature")
boxplot(humidity ~ delayed, data = flight.train, main = "Humidity", xlab = "Delay Class", ylab = "Humidity")
boxplot(precip ~ delayed, data = flight.train, main = "Precipitation", xlab = "Delay Class", ylab = "Precipitation")
boxplot(snow ~ delayed, data = flight.train, main = "Snow", xlab = "Delay Class", ylab = "Snow")
boxplot(windspeed ~ delayed, data = flight.train, main = "Wind Speed", xlab = "Delay Class", ylab = "Wind Speed")
boxplot(sealevelpressure ~ delayed, data = flight.train, main = "Pressure", xlab = "Delay Class", ylab = "Sea-Level Pressure")
boxplot(cloudcover ~ delayed, data = flight.train, main = "Cloud Cover", xlab = "Delay Class", ylab = "Cloud Cover")
boxplot(visibility ~ delayed, data = flight.train, main = "Visibility", xlab = "Delay Class", ylab = "Visibility")
par(mfrow = c(1, 1))
# 9. Model -------------------------------------------------------------------------------------------------------------------
library(MASS)
library(e1071)
library(class)
library(pROC)
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
library(sandwich)
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
TrainErr <- NULL
TestErr <- NULL
# LDA----------------------------------------------
mod1 <- lda(flight.train[, predictor.names], flight.train$delayed)
pred1 <- predict(mod1, flight.train[, predictor.names])$class
TrainErr <- c(TrainErr, mean(pred1 != flight.train$delayed))
pred1test.output <- predict(mod1, flight.test[, predictor.names])
pred1test <- pred1test.output$class
TestErr <- c(TestErr, mean(pred1test != flight.test$delayed))
tab1 <- table(
Predicted = factor(pred1test, levels = c(0, 1)),
Actual = factor(flight.test$delayed, levels = c(0, 1))
)
print(tab1)
## Actual
## Predicted 0 1
## 0 4569 2268
## 1 1889 2133
# QDA----------------------------------------------
mod2 <- qda(flight.train[, predictor.names], flight.train$delayed)
pred2 <- predict(mod2, flight.train[, predictor.names])$class
TrainErr <- c(TrainErr, mean(pred2 != flight.train$delayed))
pred2test.output <- predict(mod2, flight.test[, predictor.names])
pred2test <- pred2test.output$class
TestErr <- c(TestErr, mean(pred2test != flight.test$delayed))
tab2 <- table(
Predicted = factor(pred2test, levels = c(0, 1)),
Actual = factor(flight.test$delayed, levels = c(0, 1))
)
print(tab2)
## Actual
## Predicted 0 1
## 0 3668 2216
## 1 2790 2185
# NAIVE BAYES----------------------------------------------
nb.train <- flight.train
nb.train$delayed <- factor(nb.train$delayed, levels = c(0, 1))
mod3 <- naiveBayes(nb.train[, predictor.names], nb.train$delayed)
pred3 <- predict(mod3, flight.train[, predictor.names], type = "class")
pred3.numeric <- as.numeric(as.character(pred3))
TrainErr <- c(TrainErr, mean(pred3.numeric != flight.train$delayed))
pred3test <- predict(mod3, flight.test[, predictor.names], type = "class")
pred3test.numeric <- as.numeric(as.character(pred3test))
TestErr <- c(TestErr, mean(pred3test.numeric != flight.test$delayed))
pred3test.prob <- predict(mod3, flight.test[, predictor.names], type = "raw")
tab3 <- table(
Predicted = factor(pred3test.numeric, levels = c(0, 1)),
Actual = factor(flight.test$delayed, levels = c(0, 1))
)
print(tab3)
## Actual
## Predicted 0 1
## 0 2962 1212
## 1 3496 3189
# LOGISTIC REGRESSION----------------------------------------------
model.formula <- as.formula(paste("delayed ~", paste(predictor.names, collapse = " + ")))
mod4 <- glm(model.formula, family = binomial, data = flight.train)
print(summary(mod4))
##
## Call:
## glm(formula = model.formula, family = binomial, data = flight.train)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.118e+01 5.265e+00 13.521 < 2e-16 ***
## number_flights_matched -5.904e-04 1.767e-04 -3.341 0.000836 ***
## temp 7.785e-02 4.055e-03 19.199 < 2e-16 ***
## humidity -5.345e-03 4.638e-03 -1.152 0.249148
## precip -2.984e+00 8.726e-01 -3.419 0.000629 ***
## snow -4.447e+01 5.588e+00 -7.957 1.76e-15 ***
## windspeed -4.189e-02 1.275e-02 -3.285 0.001020 **
## sealevelpressure -7.419e-02 4.937e-03 -15.026 < 2e-16 ***
## cloudcover 5.889e-03 1.304e-03 4.515 6.33e-06 ***
## visibility 3.809e-02 6.542e-03 5.822 5.81e-09 ***
## precip_sd -2.318e+00 5.016e-01 -4.622 3.81e-06 ***
## snow_sd 1.111e+01 1.438e+00 7.728 1.09e-14 ***
## windspeed_sd 2.513e-01 3.370e-02 7.456 8.89e-14 ***
## visibility_sd 6.184e-02 1.396e-02 4.429 9.46e-06 ***
## max_windspeed -4.350e-02 7.705e-03 -5.645 1.65e-08 ***
## heathrow 5.648e-01 8.157e-02 6.924 4.39e-12 ***
## gatwick 4.775e-01 4.980e-02 9.589 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 28275 on 20395 degrees of freedom
## Residual deviance: 26095 on 20379 degrees of freedom
## AIC: 26129
##
## Number of Fisher Scoring iterations: 4
pred4.prob <- predict(mod4, flight.train, type = "response")
pred4 <- as.numeric(pred4.prob >= 0.50)
TrainErr <- c(TrainErr, mean(pred4 != flight.train$delayed))
pred4test.prob <- predict(mod4, flight.test, type = "response")
pred4test <- as.numeric(pred4test.prob >= 0.50)
TestErr <- c(TestErr, mean(pred4test != flight.test$delayed))
tab4 <- table(
Predicted = factor(pred4test, levels = c(0, 1)),
Actual = factor(flight.test$delayed, levels = c(0, 1))
)
print(tab4)
## Actual
## Predicted 0 1
## 0 4535 2223
## 1 1923 2178
# KNN TUNING----------------------------------------------
# addresses TA's concerns of data leakage; chooses the best k without using the 2024 test data
# uses earlier observations 2022-2023 to train KNN
# uses later observations 2022-2023 to compare k values
# selects k w lowest error
# trains finall KNN on all 2022-2023 obs
# evalutes it on 2024 data
# fold 1: train on 2022, validate on Jan-Jun 2023
# fold 2: train on Jan 2022-Jun 2023, validate on Jul-Dec 2023
folds <- data.frame(
train.end = c(202212, 202306),
validation.start = c(202301, 202307),
validation.end = c(202306, 202312)
)
# test values of k from 5 to 75
k.values <- seq(5, 75, by = 5)
knn.validation <- data.frame()
for (f in 1:nrow(folds)) {
# create training and validation datasets
fold.train <- flight.train[flight.train$reporting_period <= folds$train.end[f], ]
fold.validation <- flight.train[
flight.train$reporting_period >= folds$validation.start[f] &
flight.train$reporting_period <= folds$validation.end[f], ]
# calculate scaling values using fold training data only
train.means <- sapply(fold.train[, continuous.predictors], mean, na.rm = TRUE)
train.sds <- sapply(fold.train[, continuous.predictors], sd, na.rm = TRUE)
train.sds[is.na(train.sds) | train.sds == 0] <- 1
# scale predictors and add unscaled airport indicators
x.train <- cbind(
scale(fold.train[, continuous.predictors], center = train.means, scale = train.sds), fold.train[, dummy.predictors])
x.validation <- cbind(
scale(fold.validation[, continuous.predictors], center = train.means, scale = train.sds), fold.validation[, dummy.predictors])
y.train <- factor(fold.train$delayed, levels = c(0, 1))
# test each value of k
for (k in k.values) {
prediction <- knn(train = x.train, test = x.validation, cl = y.train, k = k)
validation.error <- mean(
as.numeric(as.character(prediction)) != fold.validation$delayed)
knn.validation <- rbind(knn.validation, data.frame(Fold = f, K = k, ValidationError = validation.error))}}
# average error across two folds
knn.summary <- aggregate(ValidationError ~ K, data = knn.validation, FUN = mean)
# select k with the lowest error
best.k <- knn.summary$K[which.min(knn.summary$ValidationError)]
# plot average error
plot(
knn.summary$K, knn.summary$ValidationError,
type = "b",
xlab = "Number of Neighbors (k)",
ylab = "Mean Validation Error",
main = "KNN Time-Based Validation"
)
print(knn.summary)
## K ValidationError
## 1 5 0.4040346
## 2 10 0.4006973
## 3 15 0.4006388
## 4 20 0.3986138
## 5 25 0.3992907
## 6 30 0.3984178
## 7 35 0.3953143
## 8 40 0.3956366
## 9 45 0.3954518
## 10 50 0.3954481
## 11 55 0.3967302
## 12 60 0.3986590
## 13 65 0.3997526
## 14 70 0.4002108
## 15 75 0.4006690
print(best.k)
## [1] 35
# FINAL KNN MODEL----------------------------------------------
# estimate scaling using all 2022-2023 training observations
train.means <- sapply(flight.train[, continuous.predictors], mean, na.rm = TRUE)
train.sds <- sapply(flight.train[, continuous.predictors], sd, na.rm = TRUE)
train.sds[is.na(train.sds) | train.sds == 0] <- 1
flight.train.scaled <- flight.train
flight.test.scaled <- flight.test
flight.train.scaled[, continuous.predictors] <- scale(
flight.train[, continuous.predictors],
center = train.means,
scale = train.sds
)
flight.test.scaled[, continuous.predictors] <- scale(
flight.test[, continuous.predictors],
center = train.means,
scale = train.sds
)
x.train <- as.matrix(flight.train.scaled[, predictor.names])
x.test <- as.matrix(flight.test.scaled[, predictor.names])
y.train <- factor(flight.train$delayed, levels = c(0, 1))
# KNN training prediction
pred5 <- knn(x.train, x.train, y.train, k = best.k)
pred5.numeric <- as.numeric(as.character(pred5))
TrainErr <- c(TrainErr, mean(pred5.numeric != flight.train$delayed))
# final KNN test prediction
pred5test <- knn(x.train, x.test, y.train, k = best.k, prob = TRUE)
pred5test.numeric <- as.numeric(as.character(pred5test))
TestErr <- c(TestErr, mean(pred5test.numeric != flight.test$delayed))
tab5 <- table(
Predicted = factor(pred5test.numeric, levels = c(0, 1)),
Actual = factor(flight.test$delayed, levels = c(0, 1))
)
print(tab5)
## Actual
## Predicted 0 1
## 0 3889 2072
## 1 2569 2329
# 10. calculate training and testing error--------------------------------------------------------------------------------------
error.results <- data.frame(
Model = c("LDA", "QDA", "Naive Bayes", "Logistic Regression",
paste0("KNN (k = ", best.k, ")")),
TrainingError = TrainErr,
TestingError = TestErr
)
error.results$TrainingError <- round(error.results$TrainingError, 4)
error.results$TestingError <- round(error.results$TestingError, 4)
print(error.results)
## Model TrainingError TestingError
## 1 LDA 0.3700 0.3828
## 2 QDA 0.3384 0.4610
## 3 Naive Bayes 0.3759 0.4336
## 4 Logistic Regression 0.3702 0.3818
## 5 KNN (k = 35) 0.3000 0.4274
# 11. majority class benchmark--------------------------------------------------------------------------------------------
#benchmark error = 40.53%
#assigned all 2024 observations to lower delay bc delay class = 0 was more common in 2022-2023
majority.class <- as.numeric(names(which.max(table(flight.train$delayed))))
majority.prediction <- rep(majority.class, nrow(flight.test))
majority.error <- mean(majority.prediction != flight.test$delayed)
cat("Majority-class benchmark error:", round(majority.error, 4), "\n")
## Majority-class benchmark error: 0.4053
tab0 <- table(
Predicted = factor(majority.prediction, levels = c(0, 1)),
Actual = factor(flight.test$delayed, levels = c(0, 1))
)
print(tab0)
## Actual
## Predicted 0 1
## 0 6458 4401
## 1 0 0
# 12. performance metrics-------------------------------------------------------------------------------------------------------------------
model.metrics <- function(tab) {
TN <- tab["0", "0"]
FN <- tab["0", "1"]
FP <- tab["1", "0"]
TP <- tab["1", "1"]
sensitivity <- TP / (TP + FN)
specificity <- TN / (TN + FP)
precision <- ifelse(TP + FP == 0, NA, TP / (TP + FP))
f1 <- ifelse(
is.na(precision) | precision + sensitivity == 0,
NA,
2 * precision * sensitivity / (precision + sensitivity)
)
accuracy <- (TP + TN) / sum(tab)
error <- 1 - accuracy
balanced.accuracy <- mean(c(sensitivity, specificity))
return(c(
Error = error,
Accuracy = accuracy,
Sensitivity = sensitivity,
Specificity = specificity,
Precision = precision,
F1 = f1,
BalancedAccuracy = balanced.accuracy
))
}
results <- rbind(
model.metrics(tab0),
model.metrics(tab1),
model.metrics(tab2),
model.metrics(tab3),
model.metrics(tab4),
model.metrics(tab5)
)
rownames(results) <- c(
"Majority Benchmark",
"LDA",
"QDA",
"Naive Bayes",
"Logistic Regression",
paste0("KNN k=", best.k)
)
print(round(results, 4))
## Error Accuracy Sensitivity Specificity Precision F1
## Majority Benchmark 0.4053 0.5947 0.0000 1.0000 NA NA
## LDA 0.3828 0.6172 0.4847 0.7075 0.5303 0.5065
## QDA 0.4610 0.5390 0.4965 0.5680 0.4392 0.4661
## Naive Bayes 0.4336 0.5664 0.7246 0.4587 0.4770 0.5753
## Logistic Regression 0.3818 0.6182 0.4949 0.7022 0.5311 0.5124
## KNN k=35 0.4274 0.5726 0.5292 0.6022 0.4755 0.5009
## BalancedAccuracy
## Majority Benchmark 0.5000
## LDA 0.5961
## QDA 0.5322
## Naive Bayes 0.5916
## Logistic Regression 0.5986
## KNN k=35 0.5657
# 12. confusion matrix------------------------------------------------------------------------------------------------------------------------------------------
confusion.matrices <- list(
Majority_Benchmark = tab0,
LDA = tab1,
QDA = tab2,
Naive_Bayes = tab3,
Logistic_Regression = tab4,
KNN = tab5
)
print(confusion.matrices)
## $Majority_Benchmark
## Actual
## Predicted 0 1
## 0 6458 4401
## 1 0 0
##
## $LDA
## Actual
## Predicted 0 1
## 0 4569 2268
## 1 1889 2133
##
## $QDA
## Actual
## Predicted 0 1
## 0 3668 2216
## 1 2790 2185
##
## $Naive_Bayes
## Actual
## Predicted 0 1
## 0 2962 1212
## 1 3496 3189
##
## $Logistic_Regression
## Actual
## Predicted 0 1
## 0 4535 2223
## 1 1923 2178
##
## $KNN
## Actual
## Predicted 0 1
## 0 3889 2072
## 1 2569 2329