Factors linked to plane crashes since 1919

A data analysis, data visualization, and causal inference project

library(tidyverse)
library(lubridate)
library(rdrobust)

The main data source under consideration in this project was acquired from Kaggle, where Kaggle user Abe Caesar Perez had collected a tabular dataset of all historical plane crashes across the world from 1918 to 1922. The primary source of this dataset was the Bureau of Aircraft Accident Archives (B3A), which is an international nonprofit research organzation based in Geneva that aims to conduct research into and publicize institutional knowledge related to aviation accidents in order to reduce their number, frequency, and impact.

The dataset downloaded from Kaggle included information about the aircrafts themselves, the operating companies or organizations (noting that military and government flights are included as well as commercial flights), the flight phase and location at the time of the crash, the passenger and crew manifest, records of fatalities, and the causes of the crashes, if known.

This research project aims to determine if any of these factors, or others, are associated with a greater number of crash incidents or of fatalities. It is important, however, to remember that aviation accidents are themselves quite rare events, and drawing predictions about the general safety of aviation from a dataset that solely investigates accidents would be erroneous - this project does not aim to make predictions about the likelihood of future plane crashes, but to examine the existing record of plane crashes and explore possible causative factors.

plane_crash <- read.csv('plane_crashes.csv')
head(plane_crash, 10)

print(str(plane_crash))

'data.frame':   28536 obs. of  24 variables:
 $ Date            : chr  "1918-05-02" "1918-06-08" "1918-06-11" "1918-06-19" ...
 $ Time            : chr  NA NA NA NA ...
 $ Aircraft        : chr  "De Havilland DH.4" "Handley Page V/1500" "Avro 504" "De Havilland DH.4" ...
 $ Operator        : chr  "United States Signal Corps - USSC" "Handley Page Aircraft Company Ltd" "Royal Air Force - RAF" "United States Signal Corps - USSC" ...
 $ Registration    : chr  "AS-32084" "E4104" "A8544" "AS-32098" ...
 $ Flight.phase    : chr  "Takeoff (climb)" "Takeoff (climb)" "Flight" "Flight" ...
 $ Flight.type     : chr  "Test" "Test" "Training" "Military" ...
 $ Survivors       : chr  "No" "Yes" "Yes" "No" ...
 $ Crash.site      : chr  "Airport (less than 10 km from airport)" "Airport (less than 10 km from airport)" "Plain, Valley" "Airport (less than 10 km from airport)" ...
 $ Schedule        : chr  "Dayton - Dayton" "Cricklewood - Cricklewood" "Abukir - Abukir" "Wright Patterson AFB-Wright Patterson AFB" ...
 $ MSN             : chr  NA NA NA NA ...
 $ YOM             : int  NA 1918 NA NA NA NA NA NA NA 1918 ...
 $ Flight.no.      : logi  NA NA NA NA NA NA ...
 $ Crash.location  : chr  "Dayton-McCook Field Ohio" "Cricklewood London Metropolis" "Abukir (Abu Qir) Alexandria" "Wright-Patterson AFB (Dayton) Ohio" ...
 $ Country         : chr  "United States of America" "United Kingdom" "Egypt" "United States of America" ...
 $ Region          : chr  "North America" "Europe" "Africa" "North America" ...
 $ Crew.on.board   : int  2 6 2 1 NA 1 2 2 2 2 ...
 $ Crew.fatalities : int  2 5 1 1 0 1 2 2 2 0 ...
 $ Pax.on.board    : int  0 0 0 0 NA 0 0 5 0 0 ...
 $ PAX.fatalities  : int  0 0 0 0 0 0 0 5 0 0 ...
 $ Other.fatalities: int  0 0 0 0 0 0 0 0 0 0 ...
 $ Total.fatalities: int  2 5 1 1 0 1 2 7 2 0 ...
 $ Circumstances   : chr  "The single engine airplane departed Dayton-McCook Field for a local test flight. Shortly after takeoff, the air"| __truncated__ "Assembled at Cricklewood Airfield in May 1918, the aircraft departed Cricklewood for its 13th test flight, carr"| __truncated__ "The single engine aircraft was completing a local training flight when it stalled and crashed near Abukir. A cr"| __truncated__ "Lt. Frank Stuart Patterson, son and nephew of the co-founders of National Cash Register, is killed in the crash"| __truncated__ ...
 $ Crash.cause     : chr  "Technical failure" "Technical failure" "Unknown" "Technical failure" ...
NULL

keeps <- c("Date", "Time", "Aircraft", "Operator", "Registration", "Flight.phase", "Flight.type", "Survivors", "Crash.site", "Schedule", "YOM", "Country", "Region", "Crew.on.board", "Crew.fatalities", "Pax.on.board", "PAX.fatalities", "Total.fatalities", "Crash.cause")
plane_crash_update_1 <- subset(plane_crash, select = keeps)
# remove any rows from the data frame with missing data and see how many observations are left: 
crash_data <- plane_crash_update_1[complete.cases(plane_crash_update_1), ]

crash_data$Date <- ymd(crash_data$Date) # convert Date column to date format using the package `lubridate`

crash_data$Crash.cause <- as.factor(crash_data$Crash.cause)
print(str(crash_data))

'data.frame':   11408 obs. of  19 variables:
 $ Date            : Date, format: "1919-05-01" "1919-06-15" ...
 $ Time            : chr  "5H 30M 0S" "8H 40M 0S" "8H 45M 0S" "13H 0M 0S" ...
 $ Aircraft        : chr  "De Havilland DH.9" "Vickers FB.27 Vimy Commercial" "De Havilland DH.4" "Vickers Viking (Serie I/II/III & IV)" ...
 $ Operator        : chr  "Aircraft Transport %26 Travel - AT%26T" "Vickers-Armstrongs Ltd" "Aircraft Transport %26 Travel - AT%26T" "Vickers-Armstrongs Ltd" ...
 $ Registration    : chr  "G-EAAA" "G-EAAR" "G-EAEW" "G-EAOV" ...
 $ Flight.phase    : chr  "Flight" "Landing (descent or approach)" "Landing (descent or approach)" "Flight" ...
 $ Flight.type     : chr  "Postal (mail)" "Private" "Scheduled Revenue Flight" "Positioning" ...
 $ Survivors       : chr  "No" "Yes" "Yes" "No" ...
 $ Crash.site      : chr  "Plain, Valley" "Plain, Valley" "Airport (less than 10 km from airport)" "Plain, Valley" ...
 $ Schedule        : chr  "Hendon - Bournemouth" "Saint John's - Clifden" "Croydon – Newcastle" "Brooklands - Paris" ...
 $ YOM             : int  1919 1919 1919 1919 1919 1920 1919 1919 1920 1920 ...
 $ Country         : chr  "United Kingdom" "Ireland" "United Kingdom" "France" ...
 $ Region          : chr  "Europe" "Europe" "Europe" "Europe" ...
 $ Crew.on.board   : int  1 2 1 1 4 1 4 1 1 0 ...
 $ Crew.fatalities : int  1 0 1 1 0 0 0 0 0 0 ...
 $ Pax.on.board    : int  0 0 1 0 1 3 3 2 0 0 ...
 $ PAX.fatalities  : int  0 0 0 0 0 0 0 0 0 0 ...
 $ Total.fatalities: int  1 0 1 1 0 0 0 0 0 0 ...
 $ Crash.cause     : Factor w/ 6 levels "Human factor",..: 5 5 3 5 3 3 1 3 3 5 ...
NULL

crash_hours <- vector()
crash_minutes <- vector()
for (row in 1:nrow(crash_data)) {
  hour <- str_split_i(crash_data[row,2], " ", 1)
  new_hour <- gsub('[H]', '', hour)
  minute <- str_split_i(crash_data[row,2], " ", 2)
  new_minute <- gsub('[M]', '', minute)
  crash_hours <- c(crash_hours, as.numeric(new_hour))
  crash_minutes <- c(crash_minutes, as.numeric(new_minute))
}

crash_data$Hour <- crash_hours
crash_data$Minute <- crash_minutes
print(str(crash_data))

'data.frame':   11408 obs. of  21 variables:
 $ Date            : Date, format: "1919-05-01" "1919-06-15" ...
 $ Time            : chr  "5H 30M 0S" "8H 40M 0S" "8H 45M 0S" "13H 0M 0S" ...
 $ Aircraft        : chr  "De Havilland DH.9" "Vickers FB.27 Vimy Commercial" "De Havilland DH.4" "Vickers Viking (Serie I/II/III & IV)" ...
 $ Operator        : chr  "Aircraft Transport %26 Travel - AT%26T" "Vickers-Armstrongs Ltd" "Aircraft Transport %26 Travel - AT%26T" "Vickers-Armstrongs Ltd" ...
 $ Registration    : chr  "G-EAAA" "G-EAAR" "G-EAEW" "G-EAOV" ...
 $ Flight.phase    : chr  "Flight" "Landing (descent or approach)" "Landing (descent or approach)" "Flight" ...
 $ Flight.type     : chr  "Postal (mail)" "Private" "Scheduled Revenue Flight" "Positioning" ...
 $ Survivors       : chr  "No" "Yes" "Yes" "No" ...
 $ Crash.site      : chr  "Plain, Valley" "Plain, Valley" "Airport (less than 10 km from airport)" "Plain, Valley" ...
 $ Schedule        : chr  "Hendon - Bournemouth" "Saint John's - Clifden" "Croydon – Newcastle" "Brooklands - Paris" ...
 $ YOM             : int  1919 1919 1919 1919 1919 1920 1919 1919 1920 1920 ...
 $ Country         : chr  "United Kingdom" "Ireland" "United Kingdom" "France" ...
 $ Region          : chr  "Europe" "Europe" "Europe" "Europe" ...
 $ Crew.on.board   : int  1 2 1 1 4 1 4 1 1 0 ...
 $ Crew.fatalities : int  1 0 1 1 0 0 0 0 0 0 ...
 $ Pax.on.board    : int  0 0 1 0 1 3 3 2 0 0 ...
 $ PAX.fatalities  : int  0 0 0 0 0 0 0 0 0 0 ...
 $ Total.fatalities: int  1 0 1 1 0 0 0 0 0 0 ...
 $ Crash.cause     : Factor w/ 6 levels "Human factor",..: 5 5 3 5 3 3 1 3 3 5 ...
 $ Hour            : num  5 8 8 13 14 22 15 16 20 12 ...
 $ Minute          : num  30 40 45 0 0 30 0 30 45 0 ...
NULL

years <- vector()
months <- vector()
for (row in 1:nrow(crash_data)) {
  year <- format(crash_data[row,1], "%Y")
  month <- format(crash_data[row,1], "%m")
  years <- c(years, as.numeric(year))
  months <- c(months, as.numeric(month))
}

crash_data$Year <- years
crash_data$Month <- months

print(str(crash_data))

'data.frame':   11408 obs. of  23 variables:
 $ Date            : Date, format: "1919-05-01" "1919-06-15" ...
 $ Time            : chr  "5H 30M 0S" "8H 40M 0S" "8H 45M 0S" "13H 0M 0S" ...
 $ Aircraft        : chr  "De Havilland DH.9" "Vickers FB.27 Vimy Commercial" "De Havilland DH.4" "Vickers Viking (Serie I/II/III & IV)" ...
 $ Operator        : chr  "Aircraft Transport %26 Travel - AT%26T" "Vickers-Armstrongs Ltd" "Aircraft Transport %26 Travel - AT%26T" "Vickers-Armstrongs Ltd" ...
 $ Registration    : chr  "G-EAAA" "G-EAAR" "G-EAEW" "G-EAOV" ...
 $ Flight.phase    : chr  "Flight" "Landing (descent or approach)" "Landing (descent or approach)" "Flight" ...
 $ Flight.type     : chr  "Postal (mail)" "Private" "Scheduled Revenue Flight" "Positioning" ...
 $ Survivors       : chr  "No" "Yes" "Yes" "No" ...
 $ Crash.site      : chr  "Plain, Valley" "Plain, Valley" "Airport (less than 10 km from airport)" "Plain, Valley" ...
 $ Schedule        : chr  "Hendon - Bournemouth" "Saint John's - Clifden" "Croydon – Newcastle" "Brooklands - Paris" ...
 $ YOM             : int  1919 1919 1919 1919 1919 1920 1919 1919 1920 1920 ...
 $ Country         : chr  "United Kingdom" "Ireland" "United Kingdom" "France" ...
 $ Region          : chr  "Europe" "Europe" "Europe" "Europe" ...
 $ Crew.on.board   : int  1 2 1 1 4 1 4 1 1 0 ...
 $ Crew.fatalities : int  1 0 1 1 0 0 0 0 0 0 ...
 $ Pax.on.board    : int  0 0 1 0 1 3 3 2 0 0 ...
 $ PAX.fatalities  : int  0 0 0 0 0 0 0 0 0 0 ...
 $ Total.fatalities: int  1 0 1 1 0 0 0 0 0 0 ...
 $ Crash.cause     : Factor w/ 6 levels "Human factor",..: 5 5 3 5 3 3 1 3 3 5 ...
 $ Hour            : num  5 8 8 13 14 22 15 16 20 12 ...
 $ Minute          : num  30 40 45 0 0 30 0 30 45 0 ...
 $ Year            : num  1919 1919 1919 1919 1920 ...
 $ Month           : num  5 6 10 12 2 6 6 7 8 11 ...
NULL

In this next section, the relationship between number of total fatalities and the different crash causes and flight phases are considered. The first step is to pull vectors of fatality numbers for each crash cause from the full crash_data dataframe. These are pulled and named below. The number of crashes in the dataset (with complete details) for each crash cause can also be determined in this step.

fatalities_technical <- crash_data$Total.fatalities[crash_data$Crash.cause == "Technical failure"]
print(length(fatalities_technical))

[1] 2617

fatalities_weather <- crash_data$Total.fatalities[crash_data$Crash.cause == "Weather"]
print(length(fatalities_weather))

[1] 685

fatalities_unknown <- crash_data$Total.fatalities[crash_data$Crash.cause == "Unknown"]
print(length(fatalities_unknown))

[1] 1273

fatalities_sabotage <- crash_data$Total.fatalities[crash_data$Crash.cause == "Terrorism act, Hijacking, Sabotage"]
print(length(fatalities_sabotage))

[1] 198

fatalities_human <- crash_data$Total.fatalities[crash_data$Crash.cause == "Human factor"]
print(length(fatalities_human))

[1] 6342

fatalities_other <- crash_data$Total.fatalities[crash_data$Crash.cause == "Other causes"]
print(length(fatalities_other))

[1] 293

In this dataset, considering only the crashes that contain full details (no missing data), we have records of 2617 crashes caused by “Technical failure”, 685 crashes caused by “Weather”, 1273 crashes with “Unknown” causes, 198 crashes caused by “Terrorism, Hijacking, Sabotage”, 6342 crashes caused by “Human factors”, and 293 crashes with some other cause not stated previously (“Other causes”).

library(plotly)
p2 <- ggplot(crash_data, aes(x=Crash.cause, y=Total.fatalities, fill = Crash.cause)) + geom_boxplot() + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) + ylim(0,100) + theme(legend.position="none") + labs(title = "Figure 1: Distributions of Fatality Counts by Crash Cause", x = "Stated Cause of Crash", y = "Total fatality count")
ggplotly(p2)

Warning: Removed 194 rows containing non-finite outside the scale range
(`stat_boxplot()`).

These are all highly right-skewed distributions, where most of the numbers of total fatalities are closer to zero (less than 10 in all causes!) with many outliers at the higher end of the distribution, meaning that a comparison of the medians is a more robust statistic than comparing averages, as median values are more stable in the presence of outliers than mean values. The distribution of crash fatality counts for crashes caused by “Terrorism, Hijacking or Sabotage” appear to have higher fatality counts than crashes with other causes, with “Weather”-caused crashes having the next highest IQR according to the boxplots shown above.

The next section conducts a test of the medians of these distributions to determine if the fatality counts for “Terrorism, Hijacking, Sabotage-caused crashes are actually significantly different than the median of crashes with all other causes. The coin library includes a non-parametric median hypothesis test used to determine if there are significant differences in the median values between two or more independent samples. The following hypotheses and assumptions are used in this test:

Null Hypothesis (H0): The medians of all groups are equal.

Alternative Hypothesis (H1): At least one of the medians is different from the others.

Assumptions: (1) Independence of samples, (2) Data is ordinal or continuous (3) Random sampling

library(coin)

Loading required package: survival

median_test_result <- median_test(Total.fatalities ~ Crash.cause, data = crash_data)
pvalue <- pvalue(median_test_result)
print(pvalue)

[1] <2.220446e-16

Based on the p-value of the median test, the null hypothesis is rejected, and there are significant differences in the medians of the distributions of fatality counts for crashes with different causes - most likely the higher median for crashes caused by “Terrorism, Hijacking, Sabotage” as observed in the side-by-side boxplots above (Figure 1). This of course aligns with expectations; in aviation crashes with accidental origins - such as technical failures, human error, or weather - the flight crew on board would be taking action to minimize expected harm during a crash landing.

Plot data over time for United States of America to investigate use of RDA

The plot of annual plane crash fatalities in the United States from 1919-2022 shows a marked increasing rate of crash fatalities for roughly the first 60 years of this time period, and decreasing since 1980. The US annual fatalities data will be used for Regression Discontinuity Analysis (RDA) with a cutpoint at 1980, and will consider significance of changes in the distribution of plane crash fatalities before and after this cutpoint. There is also a distinct outlier in 2001, corresponding to the 2001 9-11 terror attacks which involved the hijacking and purposeful crash of four passenger flights in the eastern US. This outlier may be removed in further analysis.

Plot data with cutpoint

In addition, adding the regression lines before and after the cutpoint will further help determine what type of RDA should be applied.

ggplotly(crash_plot_cp_lm)  # display plot

`geom_smooth()` using formula = 'y ~ x'

The scatter plot with regression lines as shown reinforces that these groupings present a good candidate for RDA. The before- and after-groups show a change in slope of the regression line (positive before 1980, negative after 1980) and a change in the intercept with the cutpoint (higher intercept for the regression model with the cutpoint before 1980, lower intercept with the cutpoint after 1980).

The RDA package rdrobust enables an RDA that tests whether the slopes of the regression lines within a bandwidth on either side of the cutpoint are significantly different. To install rdrobust, use install.packages("rdrobust").

Regression Discontinuity Analysis with rdrobust

Select bandwidth

bandwidths$bws #show the bandwidths

      h (left) h (right) b (left) b (right)
mserd 13.78192  13.78192 21.58261  21.58261

h, the main bandwidth, is +/- 13.78192 years b, the bias bandwidth, is +/- 21.58261 years. This bandwidth is used for both the bias-corrected and robust estimators when RDA is performed with rdrobust.

Let’s create a new plot showing the bandwidths as designated.

Plot data with bandwidths showing

bw

[1] 13.78192

ggplotly(p6)

`geom_smooth()` using formula = 'y ~ x'

Plotting the regression models for solely the data points within the bandwidth around 1980 shows that the change in slope from positive (before 1980) to negative (after 1980) still holds. Performing the RDA will determine whether this difference is statistically significant and, thus, whether there was a significant change in airplane crash fatality counts in the United States before and after 1980.

Perform RDA before/after 1980

print(rda_v1)

Call: rdrobust

Number of Obs.                   97
BW type                       mserd
Kernel                   Triangular
VCE method                       NN

Number of Obs.                   54           43
Eff. Number of Obs.              13           14
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                  13.782       13.782
BW bias (b)                  21.583       21.583
rho (h/b)                     0.639        0.639
Unique Obs.                      54           43

Printing the rda shows that there are a total of 54 observations before the cutpoint and 43 observations after the cutpoint. The Effective Number of Observations demonstrates that 13 observations were within the bandwidth to the left of the cutpoint and 14 observations were within the bandwidth to the right of the cutpoint.

str(rda_v1)

List of 32
 $ Estimate  : num [1, 1:4] -149 -179 108 128
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:4] "tau.us" "tau.bc" "se.us" "se.rb"
 $ bws       : num [1:2, 1:2] 13.8 21.6 13.8 21.6
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:2] "h" "b"
  .. ..$ : chr [1:2] "left" "right"
 $ coef      : num [1:3, 1] -149 -179 -179
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:3] "Conventional" "Bias-Corrected" "Robust"
  .. ..$ : chr "Coeff"
 $ se        : num [1:3, 1] 108 108 128
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:3] "Conventional" "Bias-Corrected" "Robust"
  .. ..$ : chr "Std. Err."
 $ z         : num [1:3, 1] -1.38 -1.67 -1.4
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:3] "Conventional" "Bias-Corrected" "Robust"
  .. ..$ : chr "z"
 $ pv        : num [1:3, 1] 0.1664 0.0951 0.1606
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:3] "Conventional" "Bias-Corrected" "Robust"
  .. ..$ : chr "P>|z|"
 $ ci        : num [1:3, 1:2] -359.6 -390.3 -430.2 61.9 31.3 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:3] "Conventional" "Bias-Corrected" "Robust"
  .. ..$ : chr [1:2] "CI Lower" "CI Upper"
 $ beta_Y_p_l: num [1:2] 429 13
 $ beta_Y_p_r: num [1:2] 279.811 -0.772
 $ V_cl_l    : num [1:2, 1:2] 8004 987 987 128
 $ V_cl_r    : num [1:2, 1:2] 3560.3 -376.1 -376.1 57.1
 $ V_rb_l    : num [1:2, 1:2] 12137 2394 2394 507
 $ V_rb_r    : num [1:2, 1:2] 4228 -715 -715 207
 $ N         : int [1:2] 54 43
 $ N_h       : int [1:2] 13 14
 $ N_b       : int [1:2] 21 22
 $ M         : int [1:2] 54 43
 $ tau_cl    : num [1:2] 429 280
 $ tau_bc    : num [1:2] 452 272
 $ c         : num 1980
 $ p         : num 1
 $ q         : num 2
 $ bias      : num [1:2] -23.07 7.61
 $ kernel    : chr "Triangular"
 $ all       : NULL
 $ vce       : chr "NN"
 $ bwselect  : chr "mserd"
 $ level     : num 95
 $ masspoints: chr "adjust"
 $ rdmodel   : chr "Sharp RD estimates using local polynomial regression."
 $ beta_covs : NULL
 $ call      : language rdrobust(y = annual_crash_fatalities_usa$total_fatalities, x = annual_crash_fatalities_usa$Year,      c = cutpoint)
 - attr(*, "class")= chr "rdrobust"

Looking at the structure of the rda results object shows where the other important outputs are stored. Checking the coefficient, p-value, and standard error will indicate the size of the treatment effect (grouping into before and after 1980) and whether the effect was statistically significant.

Check the coefficients:

rda_v1$coef

                   Coeff
Conventional   -148.8218
Bias-Corrected -179.4990
Robust         -179.4990

The regression slope (within the bandwidth) before the cutpoint was 148.8218 less than the slope after the cutpoint. When the larger (bias or b) bandwidth was used (21.583 years instead of 13.782 years), the difference in slopes before and after the cutpoint was larger (179.4990).

The standard error of the coefficient gives an indication of how meaningful this difference is:

rda_v1$se

               Std. Err.
Conventional    107.5366
Bias-Corrected  107.5366
Robust          127.9225

-148.8218 +/- a standard error of 107.5366 does not lead to a difference equal to or close to zero. Thus, even when taking the standard error into account, there is still a relatively large difference in the regression slopes before and after the cutpoint.

The p-value will indicate whether this difference is statistically significant:

rda_v1$pv

                    P>|z|
Conventional   0.16638353
Bias-Corrected 0.09507959
Robust         0.16056156

Here, the p-value being larger than the standard significance level of 0.05 indicates that the difference in regression slopes before and after the cutpoint, but within the bandwidth, is not statistically significant.

Repeat RDA removing the 2001 9-11 terror attack fatality counts as an outlier

bw_v2

[1] 15.66757

ggplotly(p9)

`geom_smooth()` using formula = 'y ~ x'

print(rda_v2)

Call: rdrobust

Number of Obs.                   96
BW type                       mserd
Kernel                   Triangular
VCE method                       NN

Number of Obs.                   54           42
Eff. Number of Obs.              15           16
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                  15.668       15.668
BW bias (b)                  25.230       25.230
rho (h/b)                     0.621        0.621
Unique Obs.                      54           42

<!-- rnb-source-end -->


<!-- rnb-output-end -->

<!-- rnb-output-begin eyJkYXRhIjoiRXJyb3I6IGF0dGVtcHQgdG8gdXNlIHplcm8tbGVuZ3RoIHZhcmlhYmxlIG5hbWVcbiJ9 -->

Error: attempt to use zero-length variable name

<!-- rnb-output-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->

Check the coefficients:

<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-output-begin eyJkYXRhIjoiXG48IS0tIHJuYi1zb3VyY2UtYmVnaW4gZXlKa1lYUmhJam9pWUdCZ2NseHVjbVJoWDNZeUpHTnZaV1pjYmx4dVlHQmdJbjA9IC0tPlxuXG5gYGByXG5yZGFfdjIkY29lZlxuXG5gYGBcblxuPCEtLSBybmItc291cmNlLWVuZCAtLT5cbiJ9 -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxucmRhX3YyJGNvZWZcblxuYGBgIn0= -->

```r
rda_v2$coef

                   Coeff
Conventional   -138.8252
Bias-Corrected -170.7526
Robust         -170.7526

The regression slope (within the bandwidth) before the cutpoint was 138.8252 less than the slope after the cutpoint. When the larger (bias or b) bandwidth was used (25.230 years instead of 15.668 years), the difference in slopes before and after the cutpoint was larger (170.7526).

The standard error of the coefficient gives an indication of how meaningful this difference is:

rda_v2$se

               Std. Err.
Conventional    100.3321
Bias-Corrected  100.3321
Robust          118.0563

-138.8252 +/- a standard error of 100.3321 does not lead to a difference equal to or close to zero. Thus, even when taking the standard error into account, there is still a relatively large difference in the regression slopes before and after the cutpoint.

The p-value will indicate whether this difference is statistically significant:

rda_v2$pv

                   P>|z|
Conventional   0.1664637
Bias-Corrected 0.0887791
Robust         0.1480745

About the same as the first analysis. So, excluding 9-11 as an outlier did not lead to a significant result.

Final analysis - include 9-11, but make 2001 the cutpoint!

Our final RDA will reset the cutpoint to 2001, since air travel and air travel safety procedures changed in many ways after the events of the 9-11 terror attacks.

Repeat the process of grouping the data into before/after 2001, and determining the appropriate bandwidth around this year for RDA.

ggplotly(p11)  # display plot

`geom_smooth()` using formula = 'y ~ x'

The before- and after-groups show a change in slope of the regression line (positive but smaller magnitude, flatter slope, before 2001, negative with steeper slope after 2001) and a change in the intercept with the cutpoint (higher intercept for the regression model with the cutpoint before 2001, lower intercept with the cutpoint after 2001).

RDA before/after 9-11

Use library(robust) to determine bandwidths for RDA around cutpoint at 2001.

bandwidths_v3$bws #show the bandwidths

      h (left) h (right) b (left) b (right)
mserd 8.970999  8.970999 13.38957  13.38957

h, the main bandwidth, is +/- 8.970999 years b, the bias bandwidth, is +/- 13.38957 years. This bandwidth is used for both the bias-corrected and robust estimators when RDA is performed with rdrobust.

Let’s create a new plot showing the bandwidths as designated.

Plot data with bandwidths showing

bw_v3

[1] 8.970999

ggplotly(p12)

`geom_smooth()` using formula = 'y ~ x'

print(rda_v3)

Call: rdrobust

Number of Obs.                   97
BW type                       mserd
Kernel                   Triangular
VCE method                       NN

Number of Obs.                   75           22
Eff. Number of Obs.               8            9
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                   8.971        8.971
BW bias (b)                  13.390       13.390
rho (h/b)                     0.670        0.670
Unique Obs.                      75           22

Check the coefficients:

rda_v3$coef

                  Coeff
Conventional   272.8299
Bias-Corrected 338.8700
Robust         338.8700

The regression slope (within the bandwidth) before the cutpoint was 272.829 more than the slope after the cutpoint. When the larger (bias or b) bandwidth was used (13.390 years instead of 8.971 years), the difference in slopes before and after the cutpoint was larger (338.8700).

The standard error of the coefficient gives an indication of how meaningful this difference is:

rda_v3$se

               Std. Err.
Conventional    253.6951
Bias-Corrected  253.6951
Robust          302.8105

The p-value will indicate whether this difference is statistically significant:

rda_v3$pv

                   P>|z|
Conventional   0.2821848
Bias-Corrected 0.1816352
Robust         0.2631047

Here, the p-value being larger than the standard significance level of 0.05 indicates that the difference in regression slopes before and after the cutpoint, but within the bandwidth, is not statistically significant.

This is somewhat surprising, given the great labor and analytical effort to improve airline safety in the US following the 9-11 attacks, however with closer examination of Fig. 12 above it can be observed that in the decade leading up to 9-11 (within the ~9 year bandwidth), the airline crash fatality counts were already trending downward for that time period, and that the 2001 outlier value with extremely high airline crash fatality counts was indeed still an outlier for the data within the bandwidth before the 2001 cutpoint.

Thus, the improvement in airline safety and reduction in crash fatalities that has been observed since the mid-1970s cannot be associated directly with the time period following the 2001 9-11 terror attacks, and likely has causes with effects spread out in time, such as a gradual improvement in radar and sensor technology that led to better flight path management, and engineering developments that made unscheduled landings less destructive to the aircraft and crew/passengers onboard.

Further research into key innovations in air travel safety technology and air traffic control systems could reveal a more suitable cutpoint selection that does lead to a significant before/after change in airline crash fatalities, which we do indeed observe in the rising and falling fatality counts graphs from Fig. 2 to Fig. 12. Additional investigation into the rate of plane crash incidents (relative to the number of “successful” flights) could also show a similar rising/falling pattern worth investigating, or may show a consistent decline in the rate of plane crashes with fatalities and no additional temporal changes. Research consistently shows that air travel is far safer than other forms of travel, so focusing on the rare outcomes of crashes presents on its face a difficult challenge due to the limitated observational record of such incidents.