1 Chi-Square Analysis: Bird Strike Rates

1.1 Data

Source: Korea Times (Dec 30, 2024), citing Ministry of Land, Infrastructure and Transport. Period: 2019 - August 2024.

airport_data <- data.frame(
  Airport      = c("Muan", "Gimpo", "Jeju", "Gimhae", "Daegu", "Cheongju"),
  BirdStrikes  = c(10, 140, 119, 147, 38, 33),
  TotalFlights = c(11004, 757479, 926699, NA, NA, NA)
)

airports_with_flights <- airport_data[!is.na(airport_data$TotalFlights), ]
airports_with_flights$StrikeRate_pct <- airports_with_flights$BirdStrikes /
  airports_with_flights$TotalFlights * 100

knitr::kable(airports_with_flights, digits = 4,
             col.names = c("Airport", "Bird Strikes", "Total Flights", "Rate (%)"),
             caption = "Bird strike data for Korean airports (2019-Aug 2024)")
Bird strike data for Korean airports (2019-Aug 2024)
Airport Bird Strikes Total Flights Rate (%)
Muan 10 11004 0.0909
Gimpo 140 757479 0.0185
Jeju 119 926699 0.0128

1.2 Contingency Table (Muan vs. Gimpo + Jeju)

muan_strikes    <- 10
muan_flights    <- 11004
muan_no_strikes <- muan_flights - muan_strikes

other_strikes    <- 140 + 119
other_flights    <- 757479 + 926699
other_no_strikes <- other_flights - other_strikes

contingency <- matrix(
  c(muan_strikes, muan_no_strikes,
    other_strikes, other_no_strikes),
  nrow = 2, byrow = TRUE,
  dimnames = list(
    Airport = c("Muan", "Gimpo + Jeju"),
    Outcome = c("Bird Strike", "No Bird Strike")
  )
)
knitr::kable(contingency, format.args = list(big.mark = ","),
             caption = "2 x 2 contingency table")
2 x 2 contingency table
Bird Strike No Bird Strike
Muan 10 10,994
Gimpo + Jeju 259 1,683,919

1.3 Chi-Square Test

R’s chisq.test() applies Yates’ continuity correction by default for 2x2 tables.

chi_result <- chisq.test(contingency)
chi_result
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  contingency
## X-squared = 34.661, df = 1, p-value = 3.924e-09

Manuscript values: \(\chi^2 = 34.66\), \(df = 1\), \(p < 0.001\).

print(paste0("Computed chi-square: ", round(chi_result$statistic, 2), "  (manuscript: 34.66)"))
## [1] "Computed chi-square: 34.66  (manuscript: 34.66)"
print(paste0("P-value: ", format(chi_result$p.value, scientific = TRUE, digits = 2), "  (manuscript: < 0.001)"))
## [1] "P-value: 3.9e-09  (manuscript: < 0.001)"

1.4 Odds Ratio and 95% Confidence Interval

odds_muan  <- muan_strikes / muan_no_strikes
odds_other <- other_strikes / other_no_strikes
OR         <- odds_muan / odds_other

se_log_OR <- sqrt(1/muan_strikes + 1/muan_no_strikes +
                  1/other_strikes + 1/other_no_strikes)
CI_lower  <- exp(log(OR) - 1.96 * se_log_OR)
CI_upper  <- exp(log(OR) + 1.96 * se_log_OR)

print(paste0("Odds Ratio: ", round(OR, 2)))
## [1] "Odds Ratio: 5.91"
print(paste0("95% CI: [", round(CI_lower, 2), ", ", round(CI_upper, 2), "]"))
## [1] "95% CI: [3.14, 11.13]"

Manuscript values: OR = 5.91, 95% CI: 3.14-11.13.

1.5 Rate Comparison

muan_rate  <- muan_strikes / muan_flights * 100
other_rate <- other_strikes / other_flights * 100
rate_ratio <- muan_rate / other_rate

print(paste0("Muan rate: ", round(muan_rate, 4), "%"))
## [1] "Muan rate: 0.0909%"
print(paste0("Combined rate: ", round(other_rate, 4), "%"))
## [1] "Combined rate: 0.0154%"
print(paste0("Rate ratio: ", round(rate_ratio, 1), "x"))
## [1] "Rate ratio: 5.9x"

How the combined rate is calculated: \[\text{Combined rate} = \frac{140 + 119}{757{,}479 + 926{,}699} = \frac{259}{1{,}684{,}178} = 0.0154\% \approx 0.015\%\]

1.6 Expected vs. Observed

expected_muan <- muan_flights * (other_strikes / other_flights)

print(paste0("Observed strikes at Muan: ", muan_strikes))
## [1] "Observed strikes at Muan: 10"
print(paste0("Expected at comparison rate: ", round(expected_muan, 1)))
## [1] "Expected at comparison rate: 1.7"
print(paste0("Excess: ", round((muan_strikes - expected_muan) / expected_muan * 100, 0), "%"))
## [1] "Excess: 491%"

2 Bird Strike Incidents at South Korean Airports (Figure 3)

Data from 2022-2024, sourced from Seou1. Values extracted from the published figure.

bird_data <- data.frame(
  Year    = rep(2022:2024, each = 5),
  Airport = rep(c("Gimhae", "Gimpo", "Incheon", "Jeju", "Muan"), 3),
  Rate_per_10k = c(
    4.2, 7.8, 3.0, 2.8, 14.3,
    3.9, 6.3, 3.8, 2.5, 10.3,
    3.5, 5.7, 3.0, 2.3, 22.0
  ),
  Incidents = c(
    58, 47, 35, 48, 1,
    50, 50, 55, 42, 2,
    53, 103, 42, 38, 7
  ),
  stringsAsFactors = FALSE
)

knitr::kable(bird_data, digits = 1,
             caption = "Bird strike data by airport and year (2022-2024)")
Bird strike data by airport and year (2022-2024)
Year Airport Rate_per_10k Incidents
2022 Gimhae 4.2 58
2022 Gimpo 7.8 47
2022 Incheon 3.0 35
2022 Jeju 2.8 48
2022 Muan 14.3 1
2023 Gimhae 3.9 50
2023 Gimpo 6.3 50
2023 Incheon 3.8 55
2023 Jeju 2.5 42
2023 Muan 10.3 2
2024 Gimhae 3.5 53
2024 Gimpo 5.7 103
2024 Incheon 3.0 42
2024 Jeju 2.3 38
2024 Muan 22.0 7 
airport_colors <- c(
  "Gimhae"  = "#2B5F5F",
  "Gimpo"   = "#5C7A3A",
  "Incheon" = "#8CA63A",
  "Jeju"    = "#C4CC70",
  "Muan"    = "#1B2A4A"
)

airports <- c("Gimhae", "Gimpo", "Incheon", "Jeju", "Muan")
years <- 2022:2024
pch_val <- 8

par(mfrow = c(1, 2), mar = c(4, 5, 3, 1), oma = c(0, 0, 0, 8), xpd = FALSE)

plot(NULL, xlim = c(2021.8, 2024.2), ylim = c(0, 24),
     xlab = "Year", ylab = "Incidents per 10k Operations",
     main = "A    Bird-Strike Incidents per 10,000\n     Operations by Airport",
     xaxt = "n", las = 1, cex.main = 1.0)
axis(1, at = years, labels = years)
abline(h = seq(0, 20, 5), col = "grey90", lty = 1)

for (ap in airports) {
  sub <- bird_data[bird_data$Airport == ap, ]
  lwd_val <- if (ap == "Muan") 2.5 else 1.5
  lines(sub$Year, sub$Rate_per_10k, col = airport_colors[ap], lwd = lwd_val)
  points(sub$Year, sub$Rate_per_10k, col = airport_colors[ap], pch = pch_val, cex = 1.2)
}

plot(NULL, xlim = c(2021.8, 2024.2), ylim = c(0, 110),
     xlab = "Year", ylab = "Incidents",
     main = "B    Bird-Strike Incidents by Airport",
     xaxt = "n", las = 1, cex.main = 1.0)
axis(1, at = years, labels = years)
abline(h = seq(0, 100, 25), col = "grey90", lty = 1)

for (ap in airports) {
  sub <- bird_data[bird_data$Airport == ap, ]
  lwd_val <- if (ap == "Muan") 2.5 else 1.5
  lines(sub$Year, sub$Incidents, col = airport_colors[ap], lwd = lwd_val)
  points(sub$Year, sub$Incidents, col = airport_colors[ap], pch = pch_val, cex = 1.2)
}

par(xpd = NA)
legend("right", inset = c(-0.45, 0),
       legend = airports,
       col    = airport_colors[airports],
       lwd    = c(1.5, 1.5, 1.5, 1.5, 2.5),
       pch    = pch_val,
       cex    = 0.85, bty = "n",
       title  = "Airport", title.font = 2)

3 Fault Tree Analysis

3.1 Basic Event Probabilities

From published aviation safety literature (see Table 1 in manuscript).

event_names <- c(
  "Hydraulic failure",
  "Pilot reaction delay",
  "High touchdown speed",
  "Low runway friction",
  "Initiating go around",
  "Unavailable thrust"
)

probs <- c(7e-7, 0.145, 0.0227, 0, 0.03, 0.09)
names(probs) <- event_names

knitr::kable(
  data.frame(Event = event_names, Probability = probs),
  caption = "Basic event probabilities (Table 1 in manuscript)"
)
Basic event probabilities (Table 1 in manuscript)
Event Probability
Hydraulic failure Hydraulic failure 7.00e-07
Pilot reaction delay Pilot reaction delay 1.45e-01
High touchdown speed High touchdown speed 2.27e-02
Low runway friction Low runway friction 0.00e+00
Initiating go around Initiating go around 3.00e-02
Unavailable thrust Unavailable thrust 9.00e-02

3.2 Fault Tree Structure (Figure 6)

Top Event: Collision with ILS Antenna Structure (OR gate, Q = 0.0253)
+-- Runway Overrun on Landing (OR gate, Q = 0.0227)
|   +-- Ineffective Braking (AND gate, Q = 1.015e-07)
|   |   +-- [1] Hydraulic Failure (p = 7e-7)
|   |   +-- [2] Pilot Reaction Delay (p = 0.145)
|   +-- Needing Extra Distance (OR gate, Q = 0.0227)
|       +-- [3] High Touchdown Speed (p = 0.0227)
|       +-- [4] Low Runway Friction (p = 0)
+-- Inadequate Pilot Performance (AND gate, Q = 0.0027)
    +-- [5] Initiating Go-Around (p = 0.03)
    +-- Inadequate Climb Performance (OR gate, Q = 0.09)
        +-- [6] Unavailable Engine/Thrust (p = 0.09)
        +-- [1] Hydraulic Failure (p = 7e-7)    <- shared event

3.3 Top Event Probability

Computed using the closed-form gate equations.

calculate_top_event_prob <- function(q_hf, q_prd, q_hts, q_lrf, q_iga, q_uet) {
  q_ib  <- q_hf * q_prd
  q_ned <- 1 - (1 - q_hts) * (1 - q_lrf)
  q_rol <- 1 - (1 - q_ib) * (1 - q_ned)
  q_icp <- 1 - (1 - q_uet) * (1 - q_hf)
  q_ipp <- q_iga * q_icp
  q_top <- 1 - (1 - q_rol) * (1 - q_ipp)
  return(q_top)
}

Q_top <- calculate_top_event_prob(
  q_hf = probs[1], q_prd = probs[2], q_hts = probs[3],
  q_lrf = probs[4], q_iga = probs[5], q_uet = probs[6]
)
print(paste0("Top event probability: ", round(Q_top, 8)))
## [1] "Top event probability: 0.02533883"

3.3.1 Intermediate Gate Verification

print(paste0("Ineffective Braking        = ", format(probs[1] * probs[2], scientific = TRUE, digits = 4), "   (Figure 6: 1.015e-07)"))
## [1] "Ineffective Braking        = 1.015e-07   (Figure 6: 1.015e-07)"
print(paste0("Needing Extra Distance     = ", round(1 - (1 - probs[3]) * (1 - probs[4]), 4), "    (Figure 6: 0.0227)"))
## [1] "Needing Extra Distance     = 0.0227    (Figure 6: 0.0227)"
q_icp <- 1 - (1 - probs[6]) * (1 - probs[1])
print(paste0("Inadequate Climb Perf      = ", round(q_icp, 6), "  (Figure 6: 0.09)"))
## [1] "Inadequate Climb Perf      = 0.090001  (Figure 6: 0.09)"
print(paste0("Inadequate Pilot Perf      = ", round(probs[5] * q_icp, 6), "  (Figure 6: 0.0027)"))
## [1] "Inadequate Pilot Perf      = 0.0027  (Figure 6: 0.0027)"
print(paste0("Runway Overrun on Landing  = ", round(1 - (1 - probs[1] * probs[2]) * (1 - 1 + (1 - probs[3]) * (1 - probs[4])), 6), "  (Figure 6: 0.0227)"))
## [1] "Runway Overrun on Landing  = 0.0227  (Figure 6: 0.0227)"
print(paste0("Top Event                  = ", round(Q_top, 6), "  (Figure 6: 0.0253)"))
## [1] "Top Event                  = 0.025339  (Figure 6: 0.0253)"

3.4 Minimal Cut Sets

mcs_list <- list(
  list(events = c(3),    label = "{High touchdown speed}"),
  list(events = c(4),    label = "{Low runway friction}"),
  list(events = c(1, 2), label = "{Hydraulic failure, Pilot reaction delay}"),
  list(events = c(5, 6), label = "{Initiating go-around, Unavailable thrust}"),
  list(events = c(5, 1), label = "{Initiating go-around, Hydraulic failure}")
)

mcs_df <- data.frame(
  MCS = paste0("MCS", seq_along(mcs_list)),
  Events = sapply(mcs_list, `[[`, "label"),
  Probability = sapply(mcs_list, function(m) prod(probs[m$events]))
)
knitr::kable(mcs_df, caption = "Minimal cut sets and their probabilities")
Minimal cut sets and their probabilities
MCS Events Probability
MCS1 {High touchdown speed} 2.27e-02
MCS2 {Low runway friction} 0.00e+00
MCS3 {Hydraulic failure, Pilot reaction delay} 1.00e-07
MCS4 {Initiating go-around, Unavailable thrust} 2.70e-03
MCS5 {Initiating go-around, Hydraulic failure} 0.00e+00

3.5 Importance Measures

The importance measures use the following definitions:

  • Fussell-Vesely (FV): \((Q_{top} - Q_0) \;/\; Q_{top}\)
  • Birnbaum: \(Q_1 - Q_0\)
  • RAW (Risk Achievement Worth): \(Q_1 - Q_{top}\), risk increase if event is guaranteed
  • RRW (Risk Reduction Worth): \(Q_{top} - Q_0\), risk decrease if event is eliminated

where \(Q_1 = P(\text{top} \mid p_i = 1)\) and \(Q_0 = P(\text{top} \mid p_i = 0)\).

calc_Q <- function(p) {
  calculate_top_event_prob(
    q_hf = p[1], q_prd = p[2], q_hts = p[3],
    q_lrf = p[4], q_iga = p[5], q_uet = p[6]
  )
}

plot_events <- c("Initiating go around", "Hydraulic failure",
                 "Pilot reaction delay", "Unavailable thrust",
                 "High touchdown speed")

n_ev <- length(probs)
Q1 <- Q0 <- numeric(n_ev)
names(Q1) <- names(Q0) <- event_names

for (i in 1:n_ev) {
  p_tmp    <- probs; p_tmp[i] <- 1; Q1[i] <- calc_Q(p_tmp)
  p_tmp    <- probs; p_tmp[i] <- 0; Q0[i] <- calc_Q(p_tmp)
}

birnbaum <- Q1 - Q0
fv       <- (Q_top - Q0) / Q_top
raw_diff <- Q1 - Q_top
rrw_diff <- Q_top - Q0

importance_df <- data.frame(
  Event    = event_names,
  FV       = fv,
  RRW      = rrw_diff,
  RAW      = raw_diff,
  Birnbaum = birnbaum
)

idx <- match(plot_events, event_names)
knitr::kable(importance_df[idx, ], digits = 6, row.names = FALSE,
             caption = "Importance measures for basic events (difference-based RAW/RRW)")
Importance measures for basic events (difference-based RAW/RRW)
Event FV RRW RAW Birnbaum
Initiating go around 0.104138 0.002639 0.085319 0.087958
Hydraulic failure 0.000005 0.000000 0.164137 0.164138
Pilot reaction delay 0.000004 0.000000 0.000001 0.000001
Unavailable thrust 0.104137 0.002639 0.026680 0.029319
High touchdown speed 0.893439 0.022639 0.974661 0.997300

3.5.1 Q1 and Q0 Detail

q_detail <- data.frame(
  Event = event_names[idx],
  p_i   = probs[idx],
  Q1    = Q1[idx],
  Q0    = Q0[idx],
  `Q1 - Q_top` = (Q1 - Q_top)[idx],
  `Q_top - Q0` = (Q_top - Q0)[idx],
  check.names = FALSE
)
knitr::kable(q_detail, digits = 8, row.names = FALSE,
             caption = "Detailed Q values for each event")
Detailed Q values for each event
Event p_i Q1 Q0 Q1 - Q_top Q_top - Q0
Initiating go around 3.00e-02 0.11065771 0.02270010 0.08531889 0.00263873
Hydraulic failure 7.00e-07 0.18947625 0.02533871 0.16413742 0.00000012
Pilot reaction delay 1.45e-01 0.02533941 0.02533873 0.00000058 0.00000010
Unavailable thrust 9.00e-02 0.05201910 0.02270012 0.02668027 0.00263871
High touchdown speed 2.27e-02 1.00000000 0.00270012 0.97466117 0.02263871

3.5.2 Importance Rankings (1 = most important)

rank_df <- data.frame(
  Event    = event_names[idx],
  FV       = rank(-fv[idx]),
  Birnbaum = rank(-birnbaum[idx]),
  RAW      = rank(-raw_diff[idx]),
  RRW      = rank(-rrw_diff[idx])
)
knitr::kable(rank_df, row.names = FALSE,
             caption = "Importance rankings among the 5 plotted events (1 = most important)")
Importance rankings among the 5 plotted events (1 = most important)
Event FV Birnbaum RAW RRW
Initiating go around 2 3 3 2
Hydraulic failure 4 2 2 4
Pilot reaction delay 5 5 5 5
Unavailable thrust 3 4 4 3
High touchdown speed 1 1 1 1

3.6 Spider Plot (Figure 7 Reproduction)

data_matrix <- matrix(c(
  fv[idx], birnbaum[idx], raw_diff[idx], rrw_diff[idx]
), nrow = length(idx), ncol = 4)
rownames(data_matrix) <- plot_events
colnames(data_matrix) <- c("FV", "Birnbaum", "RAW", "RRW")

base_angles <- c(pi/2, 0, -pi/2, pi)

color_map <- list(
  "Initiating go around"  = list(line = "#5C1010", fill = "#8B1A1A", lty = 1, lwd = 1.5),
  "Hydraulic failure"     = list(line = "#2B3260", fill = "#555588", lty = 1, lwd = 1.5),
  "Pilot reaction delay"  = list(line = "#5C2060", fill = "#886688", lty = 1, lwd = 1.5),
  "Unavailable thrust"    = list(line = "#6090B8", fill = "#8AB8D8", lty = 1, lwd = 1.5),
  "High touchdown speed"  = list(line = "#1A6B25", fill = "#2D8B3D", lty = 2, lwd = 2.2)
)

val_to_r <- function(v) (v + 0.25) / 1.25

par(mar = c(1, 2, 3, 8), xpd = TRUE)
plot(NULL, xlim = c(-1.15, 1.15), ylim = c(-1.15, 1.15), asp = 1,
     axes = FALSE, xlab = "", ylab = "")

grid_vals <- c(0, 0.25, 0.5, 0.75, 1.0)
for (gv in grid_vals) {
  r <- val_to_r(gv)
  vx <- r * cos(base_angles)
  vy <- r * sin(base_angles)
  polygon(vx, vy, border = "#888888", lwd = 0.6, col = NA)
}

for (a in base_angles) {
  lines(c(0, 1.05 * cos(a)), c(0, 1.05 * sin(a)), col = "#888888", lwd = 0.6)
}

for (gv in grid_vals) {
  r <- val_to_r(gv)
  text(-0.04, r - 0.01, format(gv, nsmall = 2), cex = 0.8, col = "#555555",
       adj = c(1, 0))
}

off <- 1.16
text(0, off, "FV", cex = 1.3, font = 2, col = "#1B1B1B")
text(off + 0.02, 0, "Birnbaum", cex = 1.3, font = 2, col = "#1B1B1B", adj = c(0, 0.5))
text(0, -off, "RAW", cex = 1.3, font = 2, col = "#1B1B1B")
text(-off - 0.02, 0, "RRW", cex = 1.3, font = 2, col = "#1B1B1B", adj = c(1, 0.5))

hts_i <- which(plot_events == "High touchdown speed")
draw_order <- c(hts_i, setdiff(seq_along(plot_events), hts_i))

for (di in draw_order) {
  ev <- plot_events[di]
  cm <- color_map[[ev]]
  vals <- data_matrix[di, ]
  r_vals <- val_to_r(vals)
  xs <- r_vals * cos(base_angles)
  ys <- r_vals * sin(base_angles)

  fill_alpha <- if (ev == "High touchdown speed") 0.10 else 0.07
  polygon(xs, ys, col = adjustcolor(cm$fill, alpha.f = fill_alpha), border = NA)
  lines(c(xs, xs[1]), c(ys, ys[1]), col = cm$line, lwd = cm$lwd, lty = cm$lty)
  points(xs, ys, col = cm$line, pch = 16, cex = 0.8)
}

title("FTA Importance Measures", cex.main = 1.4, font.main = 2, col.main = "#1B1B1B")

legend("right", inset = c(-0.32, 0),
       legend = plot_events,
       col    = sapply(plot_events, function(e) color_map[[e]]$line),
       lwd    = sapply(plot_events, function(e) color_map[[e]]$lwd),
       lty    = sapply(plot_events, function(e) color_map[[e]]$lty),
       cex    = 0.75, bty = "n")

3.7 Sensitivity Analysis (+/-50% Variation)

The manuscript states: “Varying probability estimates by +/-50% did not change the relative ranking of importance measures.”

compute_all_measures <- function(p) {
  Q <- calc_Q(p)
  n <- length(p)
  q1 <- q0 <- numeric(n)
  for (i in 1:n) {
    pt <- p; pt[i] <- 1; q1[i] <- calc_Q(pt)
    pt <- p; pt[i] <- 0; q0[i] <- calc_Q(pt)
  }
  list(
    Q_top    = Q,
    Birnbaum = q1 - q0,
    FV       = (Q - q0) / Q,
    RAW      = q1 - Q,
    RRW      = Q - q0
  )
}

get_top1 <- function(m, idx) {
  sapply(list(m$Birnbaum[idx], m$FV[idx], m$RAW[idx], m$RRW[idx]),
         function(x) plot_events[which.max(x)])
}

baseline     <- compute_all_measures(probs)
baseline_top1 <- get_top1(baseline, idx)

results <- data.frame(
  Perturbation = character(),
  Q_top = numeric(),
  Top1_Stable = character(),
  stringsAsFactors = FALSE
)

for (i in 1:n_ev) {
  if (probs[i] == 0) next
  for (mult in c(0.5, 1.5)) {
    p_v <- probs
    p_v[i] <- min(1, max(0, probs[i] * mult))
    res <- compute_all_measures(p_v)
    tops <- get_top1(res, idx)
    stable <- ifelse(identical(tops, baseline_top1), "Yes", "CHANGED")
    results <- rbind(results, data.frame(
      Perturbation = paste0(event_names[i], " x", mult),
      Q_top = round(res$Q_top, 6),
      Top1_Stable = stable,
      stringsAsFactors = FALSE
    ))
  }
}

knitr::kable(results,
             col.names = c("Perturbation", "Q_top", "#1 Rank Stable?"),
             caption = "Sensitivity of top-ranked event to +/-50% perturbations (among 5 plotted events)")
Sensitivity of top-ranked event to +/-50% perturbations (among 5 plotted events)
Perturbation Q_top #1 Rank Stable?
Hydraulic failure Hydraulic failure x0.5 0.025339 Yes
Hydraulic failure1 Hydraulic failure x1.5 0.025339 Yes
Hydraulic failure2 Pilot reaction delay x0.5 0.025339 Yes
Hydraulic failure3 Pilot reaction delay x1.5 0.025339 Yes
Hydraulic failure4 High touchdown speed x0.5 0.014019 Yes
Hydraulic failure5 High touchdown speed x1.5 0.036658 Yes
Hydraulic failure6 Initiating go around x0.5 0.024019 Yes
Hydraulic failure7 Initiating go around x1.5 0.026658 Yes
Hydraulic failure8 Unavailable thrust x0.5 0.024019 Yes
Hydraulic failure9 Unavailable thrust x1.5 0.026658 Yes

3.7.1 Tornado Plot

tornado <- data.frame(Event = character(), Q_low = numeric(),
                      Q_high = numeric(), stringsAsFactors = FALSE)
for (i in 1:n_ev) {
  if (probs[i] == 0) next
  p_lo <- probs; p_lo[i] <- probs[i] * 0.5
  p_hi <- probs; p_hi[i] <- probs[i] * 1.5
  tornado <- rbind(tornado, data.frame(
    Event = event_names[i],
    Q_low  = calc_Q(p_lo),
    Q_high = calc_Q(p_hi),
    stringsAsFactors = FALSE
  ))
}
tornado$Swing <- tornado$Q_high - tornado$Q_low
tornado <- tornado[order(tornado$Swing, decreasing = TRUE), ]

par(mar = c(5, 12, 4, 2))
bp <- barplot(
  tornado$Swing, horiz = TRUE, names.arg = tornado$Event,
  las = 1, col = "steelblue", border = "white",
  main = expression("Sensitivity of Q"[top]*" to "*pm*"50% Variation"),
  xlab = expression(Delta * " Q"[top]),
  xlim = c(0, max(tornado$Swing) * 1.15)
)
text(tornado$Swing, bp, labels = round(tornado$Swing, 5), pos = 4, cex = 0.8)

3.7.2 Detailed Sensitivity: HTS Measures Under Each Perturbation

hts_idx_local <- which(plot_events == "High touchdown speed")

detail_rows <- list()
detail_rows[[1]] <- data.frame(
  Scenario = "Baseline",
  Q_top = round(baseline$Q_top, 6),
  `Birn(HTS)` = round(baseline$Birnbaum[idx[hts_idx_local]], 6),
  `FV(HTS)` = round(baseline$FV[idx[hts_idx_local]], 6),
  `RAW(HTS)` = round(baseline$RAW[idx[hts_idx_local]], 6),
  `RRW(HTS)` = round(baseline$RRW[idx[hts_idx_local]], 6),
  check.names = FALSE, stringsAsFactors = FALSE
)

k <- 2
for (i in 1:n_ev) {
  if (probs[i] == 0) next
  for (mult in c(0.5, 1.5)) {
    p_v <- probs
    p_v[i] <- min(1, max(0, probs[i] * mult))
    res <- compute_all_measures(p_v)
    detail_rows[[k]] <- data.frame(
      Scenario = paste0(event_names[i], " x", mult),
      Q_top = round(res$Q_top, 6),
      `Birn(HTS)` = round(res$Birnbaum[idx[hts_idx_local]], 6),
      `FV(HTS)` = round(res$FV[idx[hts_idx_local]], 6),
      `RAW(HTS)` = round(res$RAW[idx[hts_idx_local]], 6),
      `RRW(HTS)` = round(res$RRW[idx[hts_idx_local]], 6),
      check.names = FALSE, stringsAsFactors = FALSE
    )
    k <- k + 1
  }
}

detail_df <- do.call(rbind, detail_rows)
knitr::kable(detail_df, row.names = FALSE,
             caption = "High touchdown speed measures under each perturbation")
High touchdown speed measures under each perturbation
Scenario Q_top Birn(HTS) FV(HTS) RAW(HTS) RRW(HTS)
Baseline 0.025339 0.99730 0.893439 0.974661 0.022639
Hydraulic failure x0.5 0.025339 0.99730 0.893442 0.974661 0.022639
Hydraulic failure x1.5 0.025339 0.99730 0.893437 0.974661 0.022639
Pilot reaction delay x0.5 0.025339 0.99730 0.893441 0.974661 0.022639
Pilot reaction delay x1.5 0.025339 0.99730 0.893438 0.974661 0.022639
High touchdown speed x0.5 0.014019 0.99730 0.807402 0.985981 0.011319
High touchdown speed x1.5 0.036658 0.99730 0.926343 0.963342 0.033958
Initiating go around x0.5 0.024019 0.99865 0.943791 0.975981 0.022669
Initiating go around x1.5 0.026658 0.99595 0.848072 0.973342 0.022608
Unavailable thrust x0.5 0.024019 0.99865 0.943791 0.975981 0.022669
Unavailable thrust x1.5 0.026658 0.99595 0.848072 0.973342 0.022608

4 Session Info

sessionInfo()
## R version 4.5.1 (2025-06-13 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=English_United Kingdom.utf8 
## [2] LC_CTYPE=English_United Kingdom.utf8   
## [3] LC_MONETARY=English_United Kingdom.utf8
## [4] LC_NUMERIC=C                           
## [5] LC_TIME=English_United Kingdom.utf8    
## 
## time zone: America/Los_Angeles
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.37     R6_2.6.1          fastmap_1.2.0     xfun_0.52        
##  [5] cachem_1.1.0      knitr_1.50        htmltools_0.5.8.1 rmarkdown_2.30   
##  [9] lifecycle_1.0.4   cli_3.6.5         sass_0.4.10       jquerylib_0.1.4  
## [13] compiler_4.5.1    tools_4.5.1       evaluate_1.0.5    bslib_0.9.0      
## [17] yaml_2.3.10       rlang_1.1.6       jsonlite_2.0.0