# Bayesian Network for Managing Runway Overruns in Aviation Safety

## Abstract

Runway excursions at landing constitute a major threat to aviation safety. Among them, runway overruns, defined as those occurrences when an aircraft departs the end of a runway, are the most frequent events. Although their occurrence rate is low, the entailed consequences may be very severe in terms of lives and aircraft damage. The main contributing factors to this event and their relationships are studied with the aid of a Bayesian network, while also modeling the nodes’ conditional distributions. Then, inferences and predictions are made for the quantities of interest. The issues uncovered suggest several operational recommendations to reduce the probability of facing a runway overrun when landing.

## I. Introduction

We present here a risk assessment model for ROs based on a Bayesian network (BN) [12

Our ultimate intention is to provide a template that can be used as a starting point for similar applications aimed at assessing the risk of ROs in different scenarios. With the proposed model, the involved stakeholders convert their safety-related data into valuable information, providing suggestions to improve safety, as we shall show. Section II describes the features of RO events. In Sec. II, we provide our BN model for RO risk assessment. We populate the model in Sec. IV with the aid of a case study. We end up with some discussion.

## II. Runway Overruns at Landing

According to Ref. [3], there were 1429 accidents worldwide entailing major (57%) or substantial (43%) aircraft damage in commercial air transport over the period of 1995–2008. Of those, almost 30% (417) were ROs: 34 of which involved fatalities. The likelihood of fatalities in such types of accidents is smaller than in other runway-related events, like incursions. However, due to their higher occurrence rate, ROs cause most aviation casualties during the landing phase. As an example, in the aforementioned period, 712 people were killed as a result of an RO. Moreover, the annual economic impact of ROs to, e.g., U.S. airlines was estimated at around 900 million U.S. dollars in Ref. [17].

In 2006, the Flight Safety Foundation (FSF) conducted a major descriptive study on runway safety [18]. After a comprehensive exploratory data analysis, they confirmed that ROs pose greater threats than other types of runway-related accidents. To mitigate such accidents, the FSF developed the “Safe Landing Guidelines” [18] summarized in Table 1.

The data suggested that the likelihood of an approach-and-landing accident increases greatly if any of the seven guidelines is not met and, even more, if two of them are not satisfied. Except for guideline d, the guidelines refer to human reliability issues (related to runway conditions at landing). The FSF recognizes guideline g as a precursor of having an RO: the risk of an accident is significantly greater if the aircraft is not slowed down to less than 80 kt by the time it reaches the point on the runway when only 2000 ft of landing distance available (LDA) remain. Typically, at 80 kt, the thrust-reverse levers are turned to reverse idle before stowing them at taxi speed.

However, to our knowledge, explicit probabilistic relationships among these variables have not yet been established. In Ref. [19], RO accident data were collected and normalized, pointing out multiple contributing factors to ROs. Reference [20] indicated that the most serious RO events could have been prevented if the Safe Landing Guidelines in Table 1 had been observed. Our approach will incorporate flight data and structured expert judgment to build a BN relating the factors potentially affecting RO events. As a consequence, we propose the adoption of several operational recommendations aimed at reducing the likelihood of having an RO.

## III. Model for Runway Overrun Risk Assessment at Landing

We provide a model to facilitate the evaluation of the probability of an overrun at a runway. We intend it to be a template that can serve as a starting point for similar applications aimed to assess the risk of ROs in different scenarios. Our final goal is to compare several runways or airports in terms of their probabilities of RO at landing as a pointer to identifying more dangerous infrastructures. This may suggest airport operators, airliners, or aviation safety agencies about runways in which special care is needed under identified operational conditions. Our model will relate overrun probabilities with potentially triggering factors, suggesting the introduction of mitigation measures. We focus on assessing the probability of RO, without paying attention to the entailed consequences. Thus, we implicitly assume that we are comparing runways with a similar operational environment, as in the case in Sec. IV. Our final discussion includes pointers on consequence assessment.

We build a BN to predict ROs based on several sources: 1) the contributing factors presented in Table 1, which summarize information available from manufacturers, operators, and aviation safety agencies [18]; 2) some of the results presented in Ref. [21] related to normal operational landing performance; 3) expert judgment concerning relations between the involved variables; and finally, 4) landing data. Our initial training dataset covers 897 operations for the same airline and aircraft type (short- and medium-haul fleets); landings took place at five runways in northern Spain airports with similar operational conditions, and sometimes with tough (rain, fog, etc.) but not extreme (in terms of ice, snow, hail, etc.) weather. Their locations are not disclosed for confidentiality reasons.

As a target variable $Y$, we have chosen the LDR (in feet) when the aircraft speed is 80 kt, as considered in guideline g. Ignoring the presence of other variables in this discussion for the moment, we are interested in estimating $Pr(Y<2000)$, which is used to compute

The numerator probabilities may be estimated from data in the FSF studies mentioned previously based on standard beta-binomial models [22]. The other nine variables included in our BN are listed in Table 2, and they are related to those appearing in the FSF guidelines [18].

We used GeNIe’s [23] built-in Bayesian search learning (BSL) algorithm [24] to build our BN. A first version of it was provided by our experts, which served as the initial solution to the algorithm. Because GeNIe does not support continuous variables, we discretized them using its built-in hierarchical algorithm [25]. We then fed BSL to obtain our solution using GeNIe’s default clustering algorithm for inference [26]. Additional information about learning BN structures can be found in Refs. [27,28].

Several experts checked the soundness of the suggested variable relations in the network in Fig. 1.§ The experts found all (appearing or missing) links pertinent, except for the possible connection between cross and tailwind, which did not seem relevant according to the used data, as confirmed by the slightly higher Bayesian information criterion (BIC) [29] network score computed when removing that link. We decided to leave it as such and check for possible correlations in specific cases, as in Sec. IV. Then, the proposed BN defines a probabilistic model relating factors connected with ROs through

We now discuss several potential model uses derived from Eq. (1), which are later illustrated in Sec. IV. As mentioned, for a given organization, we shall compute the target probability $Pr(Y<2000)$ to support the assessment of its probability of RO. We may use it to rank airliners, airports, or runways under similar environmental (weather at landing or state of the runway) and operational conditions. This probability may be obtained through

The target probability at a given runway $\underline{a}$ is obtained using Bayes’s formula

We are also interested in computing the aforementioned probabilities under certain operational conditions, as we will outline in the following. Such quantities could provide crews with crucial warnings about potentially dangerous situations. For brevity, we omit most computational details because they resemble those provided for $Pr(Y<2000|A=\underline{a})$, and we relegate the technical aspects of the following discussion to Appendix A.

We first assess the target probability at a given runway $\underset{\_}{a}$ under strong wind conditions [say, $Pr(Y<2000|A=\underline{a},B>\underline{b},C>\underline{c})$] for threshold tail and crosswind speeds $\underset{\_}{b}$ and $\underset{\_}{c}$. This helps in setting alarms to the crew, letting them know an eventual important increase in RO probability, given the wind conditions forecast. A similar approach may be used to assess the influence of the reverse thrust (RT) in runway $\underline{a}$ by computing $Pr(Y<2000|A=\underline{a},F\le \underline{f})$ for a certain cutoff maximum reverse thrust (MRT) time $\underset{\_}{f}$.

Another key aspect for operators is the influence of stability, which is assessed through the probability

It has a likelihood ratio form comparing risks at a given runway $\underline{a}$, depending on the approach stability. Should this be bigger than one, we would alert about dangerous conditions associated with an unstabilized approach, suggesting the implementation of procedures to increase situational awareness.

It is also relevant to determine the influence of runway conditions through

When the ratio is much bigger than one, it suggests the strong influence of runway conditions, pointing out toward improving situational awareness among crew if the runway is contaminated.

## IV. Template for Case Studies

We present a case study that may be used as a template for related applications aimed at assessing the risk of ROs in different scenarios with similar operational and environmental conditions. It is undertaken by an airline willing to assess RO risk at three runways designated rwy-01, rwy-02, and rwy-03. Within the framework of RO risk management, the proposed graphical model (Fig. 1) and the models and priors at various nodes identified in the following could be used as starting points for other case studies considered, e.g., for other airports, runways, companies, or countries. The case study also serves to show the operational insights and recommendations that may be derived from our model.¶

We have available a second dataset provided by an operator, which is independent of that used in Sec. III to learn the BN structure. We use it, together with expert judgment, to build the predictive distributions at the network nodes. References [30,31] provided examples in which solely expert judgment was used in building and populating a BN model for, respectively, highway and nuclear plant safety assessments. In turn, Refs. [32,33] illustrated how to learn the BN structure and parameters solely from data. Due to the nature of the RO problem, we found very few critical cases in the period considered. However, we used expert opinion as another source of information to provide forecasts and combined it with the little data available. Reference [34] included additional information about learning the node distributions. Some of the priors we use are diffuse but proper (i.e., with large variance but integrating to one); others reuse knowledge gained from data in other parts of the problem; finally, some other priors model the experts’ knowledge in relation with the physical aspects of the problem. In any case, they are object to a sensitivity analysis in Sec. IV.B.

We therefore used data from $n=266$ landings over a period of 10 months, with all referring to the same airline and aircraft type. They concern the three runways mentioned previously and shown in Fig. 2 (all with LDAs smaller than 7200 ft), as well as similar operational and environmental conditions around the airports. In these runways, should an RO happen, it would just affect the aircraft and/or passengers and crew because there are no buildings nearby. Thus, we pay no attention to consequences because they would be similar for the three involved infrastructures (in terms of fatalities, injuries, and aircraft damage or destruction), and we focus on estimating the probability of the relevant event. See our final Sec. V for a discussion on consequence assessment.

Recalling our analysis in Fig. 1, we start by investigating possible relationships between cross- and tailwind components in this case. Pearson’s correlation coefficients are, respectively, ${r}_{1}=-0.12$, ${r}_{2}=0.04$, and ${r}_{3}=0.21$ if we include zero tailwinds; after removing such null values, the correlations become ${r}_{1}=0.04$, ${r}_{2}=-0.01$, and ${r}_{3}=0.35$. None suggest strong dependence. Therefore, we proceed with our graphical model from Fig. 1 with no direct link between cross- and tailwind nodes.

We illustrate now several model uses supporting relevant operational decisions. We describe in Appendix B the predictive distributions at various nodes deployed to populate model (1). For simplicity, we drop the dependence on data, writing, e.g., $p(A)$ rather than $p(A|\mathrm{data})$. In most nodes, we shall be able to obtain only a sample from the relevant predictive distribution. When pertinent, we display the expected predictive distribution over the data histogram to suggest the appropriate model fit. All computations were performed with R and WinBUGS [35,36].

### A. Using the Model

We have developed a sampling scheme to implement the computations in Sec. III. The reported results were obtained after 1) generating samples of a size of 1,000,000 from the involved predictive distributions of the variables $A,B,\dots ,I,Y$, which are detailed in Appendix B; 2) performing the relevant calculations (probabilities and likelihood ratios) as outlined in Sec. III; and 3) replicating such a scheme 100 times to estimate the corresponding uncertainty.

First, we are interested in the overall probability of having a landing in which the landing distance remaining (LDR) is less than 2000 ft when the aircraft is at 80 kt. Through Algorithm 1, we have obtained $\widehat{Pr}(Y<2000|\mathrm{data})=0.00019$ with an associated standard deviation of $4.43\cdot {10}^{-5}$. Thus, given the load of the incumbent operator, we should expect one landing with high risk of RO at any of the three runways approximately every 400 days. Extrapolating this result to the operations of the airline in its country where it operates (assuming all runways have similar operational environments), we should expect one hazardous operation (in relation to RO) approximately every 26 days. *If this is deemed excessive by the airline, this would require risk mitigation strategies*.

We find out which of the three runways is the most dangerous one. Using Algorithm 1 with ${a}_{k}=\underline{a}$, we get the estimates in Table 3.

Rwy-01 has the largest probability, entailing roughly one hazardous operation every two years. Figures are comparable for rwy-02, whereas the target probability at rwy-03 is approximately 10 times smaller. Rwy-01 is indeed the most complicated runway from an operational viewpoint because 1) it is not equipped with an instrument landing system (ILS) and pilots have to approach it visually; 2) it has a slightly U-shaped valley form; and 3) adverse wind conditions are more frequent (see Figs. 3 and 4). *Eventually, risk management in such a runway would proceed by introducing an ILS and emphasizing wind issues in briefing manuals and arrival charts*.

In regard to the influence of wind conditions, we consider $\underline{b}\in \{10,15,20\}$ and $\underline{c}\in \{0,5,10\}$ as relevant cutoff wind speeds and assess the target probabilities in Table 4. We use the samples generated in Algorithm 1 with ${a}_{k}=\underline{a}$ for such a purpose.

Clearly, higher probabilities hold as wind conditions worsen: landing under such circumstances tends to increase the probability of an unstabilized approach. This will affect the MRT, the difference $I$, and the height at threshold (see Fig. 1), eventually increasing the risk of having less LDR at 80 kt. If that is the case, the autobrake system must be selected. *Nevertheless, under severe wind conditions (e.g.*, $\underline{b}=20\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{kt}$*and*$\underline{c}=10\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{kt}$), *“going around” and proceeding to the alternate airport would be the right decision to avoid putting safety in jeopardy. This should be briefed to crews*.

We analyze now the influence of RT. Its use is linked with policies to reduce fuel consumption and noise at airports: pilots are requested to operate RT at idle during ground roll, except when safety reasons recommend it. How long pilots operate at MRT will depend on the speed the aircraft needs to be slowed down, which is influenced by factors such as tailwind. RT cutoff values were set at $\underline{f}\in \{0,5,10\}$. Using the samples from Algorithm 1 with ${a}_{k}=\underline{a}$, we obtain Table 5.

Therefore, the longer the MRT has been operated, the higher the target probability. This indicates that, if RT has been used for too long (say, over $\underline{f}=10\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$), the aircraft might have experienced difficulties in decelerating during ground roll, risking an RO. *Thus, this variable could be used as an RO risk precursor: relevant information for flight data monitoring teams*.

We consider now the influence of approach unstability. After running Algorithm 1 with ${a}_{k}=\underline{a}$ and ${e}_{k}=\underline{e}$, we have obtained ratios (9.11, 11.56, 2.59) for the runways, respectively. Thus, more dangerous conditions will arise under an unstabilized approach, especially at rwy-01 and rwy-02. *Recommendations oriented toward increasing situational awareness should be given to crews. Innovative mitigation measures could be a notification from the flight dispatcher to the crew of such high risk level or including a warning in briefing manuals or arrival charts*.

Finally, we discuss the influence of runway conditions. Through Algorithm 1 with ${a}_{k}=\underline{a}$ and ${d}_{k}=\underline{d}$, we obtained ratios (1.08, 0.96, 0.95) for the three runways. They are relatively close to one, suggesting little influence of surface contamination on the risk of runway overrun in this case: possibly because pilots usually take extra precautions when facing adverse runway conditions. *Thus, no specific information needs to be included in briefing manuals or arrival charts*.

### B. Sensitivity Analysis

We have performed a sensitivity analysis over various priors to check for robustness of the results. We have focused on those parameters with larger posterior variances to ascertain the influence of different prior choices. As an illustration, and because of its judgmental nature, we explore here sensitivity to prior variations in relation with height at the threshold. We tried alternative priors for the $w$, $\nu $, and $\mu $, reflecting different precision regarding prior knowledge. We considered the eight combinations of $[\varphi ,\theta ,(\alpha ,\beta )]$ in Table 6.

The results obtained did not show significant changes, suggesting robustness: the estimates for the corresponding predictive distribution varied by less than 5%, whereas the target probabilities differed by less than 2% with respect to those in Table 3. Similar analyses were performed for other model parameters, confirming the robustness of our priors and results.

## V. Conclusions

Although they occur with very low probability, runway overruns constitute a major threat to aviation safety because their implications may be very severe in terms of lives and destroyed aircraft. Several international aviation organizations like the FSF, the International Civil Aviation Organization, or aircraft manufacturers are launching efforts to reduce landing overruns, investigating various safety strategies.

As cited, information available from manufacturers, operators, and safety agencies (The Boeing Company, the FSF, the Federal Aviation Administration, Honeywell, Airbus, etc.); expert judgment; and data were used to propose a BN defining a probabilistic model for assessing RO probabilities. The critical precursor event considered was whether the LDR at 80 kt was less than 2000 ft. Its probability was used to assess the RO risk and rank various aviation systems like runways, airports, and others. The methodology is general, and the probabilistic graph underlying the current network may be used in other cases, except for the possible relation between the crosswind and tailwind, which would need to be checked as explained. The node predictive distributions would also need to be reevaluated in other case studies, but the current proposals may serve as initial models when comparing homogeneous aviation systems. The current working BN can easily integrate new data, accordingly updating the node predictive distributions; see Ref. [34] for a discussion.

As with many other statistical models, one of the main limitations of the current proposal might be the scarcity of data. However, in this case, this is a fortunate circumstance because it means that few accidents have actually occurred. Nevertheless, it is believed that it is not a matter of having enough data but of using them properly in combination with the expert judgment available. Furthermore, as data accumulate, the corresponding predictive distributions would be updated.

The validity of the current methods has been thoroughly assessed. First, the BN built in Fig. 1 was chosen among a set of suitable models as the one with the lowest BIC. Second, each of the predictive distributions was validated for the goodness of fit, displaying (when pertinent) the expected predictive distribution over the data histogram. Some assessments were also made based on the predictive probabilities obtained. To wit, the results in Table 3 were in accordance with the observed data where, for instance, only one out of 80 operations took place under an unstabilized approach at rwy-02.

From the current case study, several general relevant issues were uncovered, including the following:

- The first relevant issue involves the weather conditions. Landing under windy conditions tends to increase the probability of an unstabilized approach.
- The second relevant issue involves the variables with the strongest influence. The variables having stronger influence over the LDR at 80 kt (and hence over RO) were the autobrake system and the difference $I$ at threshold.
- The third relevant issue involves ducking under. The presence of a ducking under effect was alerted against, which is not always a good procedure.
- The fourth relevant issue involves unstabilized approach. As RO accident reports highlight, much more danger arises from unstabilized approaches.
- The fifth relevant issue involves reverse thrust. Excessive use of MRT suggests difficulties for decelerating during ground roll.

Note that the focus has been on evaluating the probability of a leading precursor of RO because homogeneous runways are being dealt with. Thus, attention has not been paid to forecasting the entailed consequences, which would typically include the number of deaths and injured people, the potential destruction of the aircraft, the induced delay, and so on. Forecasting them would require additional modeling efforts. They would then be evaluated through a multiattribute loss function; see Refs. [37,38] for the details of a general framework for aviation safety risk management. With such ingredients, the acceptability of certain risks could be assessed or appropriate countermeasures (possibly referring to signals, markings, lighting, or updated procedures that could potentially reduce overruns) could be identified, minimizing the expected loss associated with their introduction.

Many safety-related issues are yet to be explored in the aviation industry. To cite just a few, the potential impact of human behavior on the performance when landing, the influences of other external factors (such as, e.g., the congestion of airports or the pressure of tight schedules), or the attrition or lack of maintenance of some aircraft are still, to a great extent, open questions. New methodological solutions, possibly following similar approaches to the one in this paper, would be necessary to gain better understanding and propose innovative solutions to such problems.

§ GeNIe uses the thickness of the arcs to indicate the strength of influence between the nodes they connect. For ease of viewing, we have further emphasized stronger intensities with darker links.

¶ In what follows, we will emphasize in italics the risk management options suggested.

## Appendix A: Probabilities and Likelihood Ratios

### A.1. Computation of $Pr(Y<2000|A=\underline{a},\text{\hspace{0.17em}}B>\underline{b},\text{\hspace{0.17em}}C>\underline{c})$

Using Bayes’s formula, and canceling common terms, this probability is proportional to

### A.2. Computation of $Pr(Y<2000|A=\underline{a},\text{\hspace{0.17em}}E=\mathrm{Unst})/Pr(Y<2000|A=\underline{a},\text{\hspace{0.17em}}E=\mathrm{Stab})$

After canceling common terms, the numerator $Pr(Y<2000,A=\underline{a},E=\underline{e})$ is approximated through

### A.3. Computation of $Pr(Y<2000|A=\underline{a},\text{\hspace{0.17em}}D=\mathrm{Contaminated})/Pr(Y<2000|A=\underline{a},\text{\hspace{0.17em}}D=\mathrm{Noncontaminated})$

The required probability is proportional to

## Appendix B: Predictive Models at Nodes

Runway usage $p(A)$ is a categorical variable with values $A\in \{1,2,3\}$ referring to the runway used. We employ a Dirichlet-multinomial model with a uniform prior for runway usage [22]. Let $({p}_{1}^{a},{p}_{2}^{a},{p}_{3}^{a})$ be the respective usage proportions of such runways. Then,

We model now the crosswind, given runway $p(B|A)$. All 266 flights entailed nonzero crosswind, with histograms in Fig. 3. For each runway, we use a gamma model $B|A=\ell \sim \mathcal{G}({\nu}_{\ell}^{b},{\nu}_{\ell}^{b}/{\mu}_{\ell}^{b})$ [39].

We respectively consider diffuse, but proper, exponential ${\nu}_{\ell}^{b}\sim \mathcal{E}(0.01)$ and inverse gamma ${\mu}_{\ell}^{b}\sim \mathcal{IG}(1,1)$ priors. There is no closed-form expression for the posteriors, but we may sample from ${\mu}_{\ell}^{b}|\mathrm{data}$ and ${\nu}_{\ell}^{b}|\mathrm{data}$, as summarized in Table B1. Based on such samples, we approximate the predictive distribution $p(b|A=\ell ,\mathrm{data})$, plotted with a dashed line in Fig. 3, suggesting a good fit.

Next, we model the tailwind, given runway $p(C|A)$. Out of the 266 flights, 125 took place with no tailwind (Table B2). Figure 4 presents the histograms for the cases with positive tailwind.

For each runway $\ell $, we model the tailwind distribution $C$ through a mixture

${\mathbb{I}}_{0}^{c}$ describes zero tailwind cases, and ${C}_{\ell}^{+}$ describes the positive ones. Based on a uniform prior for the beta-binomial model [22], we get that, a posteriori,

We deal now with modeling $p(D|A)$, referring to runway conditions at $A$, given the runway. We have available data about weather conditions for the incumbent operations (0 denotes noncontaminated, and 1 denotes contaminated). We consider a runway contaminated if rainfall over the previous 24 h exceeded $8\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}/{\mathrm{m}}^{2}$ [40]. We use a beta-binomial model with a uniform prior for ${p}_{1\ell}^{d}$, with the probability of runway $\ell $ being contaminated when landing. A posteriori, ${p}_{1\ell}^{d}|\mathrm{data}\sim \mathcal{B}e(1+{x}_{\ell}^{d},1+{n}_{\ell}-{x}_{\ell}^{d})$, where ${x}_{\ell}^{d}$ is the number of operations on a contaminated runway. When necessary, we summarize ${p}_{1\ell}^{d}$ through its posterior mean $(1+{x}_{\ell}^{d})/(2+{n}_{\ell})$. Table B4 shows the estimates and the available data in parentheses. The predictive distribution is [41]

Then, we model stability given the runway, crosswind, and tailwind through $p(E|A,B,C)$. We label ${e}_{j}=1$ (0) if the approach is unstabilized (stabilized). We use a logistic regression model for the probabilities ${p}_{\ell j}^{e}=Pr({E}_{\ell j}=1|A=\ell ,B,C)$:

Because there is no closed-form expression for the predictive distribution, we use a composition simulation scheme to obtain the required predictive samples. Symmetric 0.9 predictive intervals [22] covered observations, suggesting good performance of our model.

We focus now on modeling RT, given the runway and stability $p(F|A,E)$. Out of 266 operations, 118 entailed reverse at idle thrust. The remaining 148 landings used MRT during ground roll. Table B6 summarizes these data, factored by landing runway and approach stability. Histograms are shown in Fig. B1.

We distinguish between stabilized and unstabilized approaches. For stabilized ones, we decompose RT $F$ through the mixture

For unstabilized cases, there are no zero values (Fig. B1). We model the MRT through a gamma model, using as priors the posterior distributions obtained for the positive cases in the stabilized approach, i.e.,

We proceed as with the tailwind to obtain a sample from the predictive densities. Overlaid expected predictive densities in Fig. B1 suggest a good fit.

We consider autobrake given the runway $p(G|A)$. We use a Dirichlet-multinomial model with a uniform prior for each runway. Based on the coding (1 denotes no, 2 denotes low, and 3 denotes medium) autobrake, we have

When needed, we use estimates

See Table B9. Note that the autobrake system is mostly used in its “medium” setting at the shortest runway (rwy-01) because it provides the highest deceleration.

Next, we model $p(H|A,E)$, referring to height at the threshold, for a given runway and stability conditions. Figure B2 represents the height $H$ when aircraft overfly the threshold. Ideally, this should be done at 50 ft (vertical dashed line). For each runway, we consider that data come from a mixture of, at most, three gamma distributions: one for values close to zero (too low), and the other two for intermediate to upper values (too high),

We consider the prior

For unstabilized approaches, pilots tend to fly either slightly below or above 50 ft, according to experts. Therefore, we consider only two mixture terms:

Due to the small sample size, we use as priors the posteriors over the last two components in the stabilized case and their corresponding (renormalized) weights. Posterior estimates are shown in Table B11.

For the predictive $p(h|A=\ell ,E=1,\mathrm{data})$, we proceed through a sampling scheme. The overlaid expected predictive densities in Fig. B2 suggest a good fit.

We deal now with modeling $\mathrm{IAS}-{V}_{\mathrm{app}}$, given the stability of approach $p(I|E)$. Figure B3 represents $I=\mathrm{IAS}-{V}_{\mathrm{app}}$ at threshold, which is factored by the type of approach $E$. Ideally, $I=0$; i.e., the indicated airspeed (IAS) should coincide with ${V}_{\mathrm{app}}$, which is indicated with a vertical dashed line.

For both types of approaches, we assume a normal model $I|E\sim \mathcal{N}(\mu ,{\sigma}^{2})$ with a noninformative prior $p(\mu ,{\sigma}^{2})\propto 1/{\sigma}^{2}$. A posteriori,

Finally, we model the LDR $Y$ at 80 kt, given the relevant predecessors $p(Y|A,D,F,G,H,I)$. We plot in Fig. B4 the LDR when the aircraft is decelerated to 80 kt for each runway, indicating with a vertical dashed line the threshold distance of 2000 ft.

We use a linear regression model

The LDR $Y$ at 80 kt seems strongly influenced by the autobrake $G$ and the difference $I$ at threshold. The positive value of ${\beta}_{3}^{y,m}$ suggests that when, the autobrake is selected at medium, the landing roll will decrease, whereas the LDR will increase. The same reasoning holds for ${\beta}_{3}^{y,l}$, albeit to a lesser extent. With respect to the negative sign of ${\beta}_{5}^{y}$, this suggests that a greater IAS at threshold implies greater kinetic energy and, as a consequence, more runway is required. As for the runway condition $D$, ${\beta}_{1}^{y}$ has large positive values for rwy-02 and rwy-03. This could be due to two factors:

- The runway is severely contaminated at landing (typically with heavy rain), which is a relatively frequent phenomenon in the incumbent locations, especially in winter; this might damp the airplane during ground roll, therefore reducing RO risk.
- When pilots observe adverse weather conditions at destination, they will usually increase their situational awareness, actually leading to safer landings.

The influence of runway conditions seems smaller for rwy-01, possibly suggesting that, in some cases, runway contamination might eventually cause hydroplaning, reducing the efficacy of the deceleration process. Nevertheless, the posterior variance of ${\beta}_{1}^{y}$ implies nonnegligible uncertainty, especially for rwy-02 and rwy-03. The impact of the other variables seems less relevant. For $F$, because ${\beta}_{2}^{y}>0$, this suggests that the LDR at 80 kt is slightly increased by the use of MRT. For $H$, ${\beta}_{4}^{y}>0$ could be an indicator of the tendency of pilots to pitch down before the threshold so as to touch down as soon as possible, a maneuver called ducking under: pilots sometimes descend below the glide slope to have a larger LDA, especially at short runways. However, when ducking under is performed incorrectly, it could have an opposite effect due to the increase in kinetic energy and the ground effect, as well as the reduction in induced drag. Thus, pilots should be advised to avoid an excessive or inappropriate use of ducking under in these runways.

A closed-form expression for the required predictive distribution is not available, but we may draw from the runway predictive distribution. Symmetric 0.9 predictive intervals covered observations, suggesting good performance of our model.

## Acknowledgments

This work is supported by the Spanish Ministry of Economy and Innovation project MTM2017-86875-C3-1-R and the AXA-ICMAT Chair on Adversarial Risk Analysis. Data concerning weather conditions at the case study runways have been provided by the Spanish Agencia Estatal de Meteorologa.

## References

[1] “Statistical Summary of Commercial Jet Airplanes Accidents 1959–2017,” Tech. Rept., The Boeing Company, Worldwide Operations, Seattle, Washington, D.C., Oct. 2018, http://www.boeing.com/resources/boeingdotcom/company/about_bca/pdf/statsum.pdf.

[2] , “A Review of Research on Risk and Safety Modelling in Civil Aviation,”

**Journal of Air Transport Management**, Vol. 14, No. 4, 2008, pp. 213–220. doi:https://doi.org/10.1016/j.jairtraman.2008.04.008[3] “Reducing the Risk of Runway Excursions,” Tech. Rept., Flight Safety Foundation, Alexandria, VA, May 2009, https://flightsafety.org/files/RERR/fsf-runway-excursions-report.pdf.

[4] , “Lateral Runway Excursion upon Landing,”

**Safety First**, Vol. 20, July 2015, pp. 1–14.[5] , “A Review of some Technological Aids to Support the Flight Safety Foundation Runway Safety Initiative (RSI),”

**FSF 61st Annual International Air Safety Seminar (IASS)**, Honolulu, Hawaii, 2008, pp. 1–24.[6] Transportation Research Board, and National Academies of Sciences, Engineering, and Medicine,

**Improved Models for Risk Assessment of Runway Safety Areas**, National Academies Press, Washington, D.C., 2011, https://www.nap.edu/catalog/13635/improved-models-for-risk-assessment-of-runway-safety-areas. doi:https://doi.org/10.17226/13635[7] , “A New Airline Safety Index,”

**Transportation Research Part B: Methodological**, Vol. 38, No. 4, 2004, pp. 369–383. doi:https://doi.org/10.1016/S0191-2615(03)00047-X[8] , “Predicting Risk of Near Midair Collisions in Controlled Airspace,”

**Transportation Research Part B: Methodological**, Vol. 25, No. 4, 1991, pp. 237–252. doi:https://doi.org/10.1016/0191-2615(91)90006-5[9] , “Contrasting Safety Assessments of a Runway Incursion Scenario: Event Sequence Analysis Versus Multi-Agent Dynamic Risk Modelling,”

**Reliability Engineering and System Safety**, Vol. 109, Jan. 2013, pp. 133–149. doi:https://doi.org/10.1016/j.ress.2012.07.002[10] , “Improving Safety of Runway Overrun Through the Correct Numerical Evaluation of Rutting in Cleared and Graded Areas,”

**Safety Science**, Vol. 62, Feb. 2014, pp. 326–338. doi:https://doi.org/10.1016/j.ssci.2013.09.008[11] , “Runway Veer-Off Accidents: Quantitative Risk Assessment and Risk Reduction Measures,”

**Safety Science**, Vol. 104, April 2018, pp. 157–163. doi:https://doi.org/10.1016/j.ssci.2018.01.010[12] ,

**Bayesian Networks and Decision Graphs**, 2nd ed., Springer, New York, 2007. doi:https://doi.org/10.1007/978-0-387-68282-2[13] , “Reverend Bayes on Inference Engines: A Distributed Hierarchical Approach,”

**Second National Conference on Artificial Intelligence**, Carnegie–Mellon Univ., Pittsburgh, PA, 1982, pp. 133–136.[14] , “Fusion, Propagation, and Structuring in Belief Networks,”

**Artificial Intelligence**, Vol. 29, No. 3, 1986, pp. 241–288. doi:https://doi.org/10.1016/0004-3702(86)90072-X[15] , “Towards a Causal Model for Air Transport Safety—An Ongoing Research Project,”

**Safety Science**, Vol. 44, No. 8, 2006, pp. 657–673.[16] , “Prediction of Aircraft Safety Incidents Using Bayesian Inference and Hierarchical Structures,”

**Safety Science**, Vol. 104, April 2018, pp. 216–230.[17] “Reducing the Risk of Runway Incursions and Excursions,” Honeywell Aerospace, Alexandria, VA, 2012, http://www.willis.se/documents/Willis_Aerospace_Aviation_Products_Market_Review_2012.pdf.

[18] , “Keys to a Safe Arrival,”

**AeroSafety World**, Oct. 2011, pp. 14–17.[19] , “An Improved Methodology for Assessing Risk in Aircraft Operations at Airports Applied to Runway Overruns,”

**Safety Science**, Vol. 42, No. 10, 2004, pp. 891–905. doi:https://doi.org/10.1016/j.ssci.2004.04.002[20] , “Overrun Breakdown,”

**AeroSafety World**, Nov. 2012.[21] , “Study of Normal Operational Landing Performance on Subsonic, Civil, Narrow-Body Jet Aircraft During Instrument Landing System Approaches,” Federal Aviation Administration TR DOT/FAA/AR-07/7, March 2007.

[22] ,

**Statistical Decision Theory**,Kendall’s Library of Statistics , Vol. 9, Arnold Publishers, London, 2000.[23] GeNIe (Graphical Network Interface), Software Package, BayesFusion, Pittsburgh, PA, 2019.

[24] , “A Bayesian Method for the Induction of Probabilistic Networks from Data,”

**Machine Learning**, Vol. 9, No. 4, 1992, pp. 309–347. doi:https://doi.org/10.1007/BF00994110[25] ,

**Cluster Analysis for Applications**,Probability and Mathematical Statistics , Vol. 19, Academic Press, New York, 1973. doi:https://doi.org/10.1016/C2013-0-06161-0[26] , “Local Computations with Probabilities on Graphical Structures and Their Application to Expert Systems,”

**Journal of the Royal Statistical Society. Series B (Methodological)**, Vol. 50, No. 2, 1988, pp. 157–224. doi:https://doi.org/10.1111/j.2517-6161.1988.tb01721.x[27] , “Learning Causal Bayesian Network Structures from Experimental Data,”

**Journal of the American Statistical Association**, Vol. 103, No. 482, 2008, pp. 778–789. doi:https://doi.org/10.1198/016214508000000193[28] ,

**Approximation Methods for Efficient Learning of Bayesian Networks**,Frontiers in Artificial Intelligence and Applications , Vol. 168, IOS Press, Amsterdam, 2008, https://books.google.es/books?id=37kO0psfWkcC [retrieved 2019].[29] ,

**Pattern Recognition and Machine Learning**, Springer, New York, 2006.[30] , “Highway and Road Probabilistic Safety Assessment Based on Bayesian Network Models,”

**Computer-Aided Civil and Infrastructure Engineering**, Vol. 32, No. 5, 2017, pp. 379–396. doi:https://doi.org/10.1111/mice.12248[31] , “Designing a Bayesian Network for Preventive Maintenance from Expert Opinions in a Rapid and Reliable Way,”

**Reliability Engineering and System Safety**, Vol. 91, No. 7, 2006, pp. 849–856. doi:https://doi.org/10.1016/j.ress.2005.08.007[32] , “Learning Bayesian Networks from Data: An Information-Theory Based Approach,”

**Artificial Intelligence**, Vol. 137, Nos. 1–2, 2002, pp. 43–90. doi:https://doi.org/10.1016/S0004-3702(02)00191-1[33] , “Exploiting Missing Clinical Data in Bayesian Network Modeling for Predicting Medical Problems,”

**Journal of Biomedical Informatics**, Vol. 41, No. 1, 2008, pp. 1–14. doi:https://doi.org/10.1016/j.jbi.2007.06.001[34] , “Parameter Updating in a Bayes Network,”

**Journal of the American Statistical Association**, Vol. 87, No. 420, 1992, pp. 1109–1115. doi:https://doi.org/10.1080/01621459.1992.10476266[35] “R: A Language and Environment for Statistical Computing,” R Core Team, R Foundation for Statistical Computing, Vienna, 2019, http://www.R-project.org/.

[36] , “WinBUGS—A Bayesian Modelling Framework: Concepts, Structure, and Extensibility,”

**Statistics and Computing**, Vol. 10, No. 4, 2000, pp. 325–337. doi:https://doi.org/10.1023/A:1008929526011[37] , “A Framework for Risk Management Decisions in Aviation Safety at State Level,”

**Reliability Engineering and System Safety**, Vol. 179, Nov. 2018, pp. 74–82. doi:https://doi.org/10.1016/j.ress.2016.12.002[38] , “Forecasting and Assessing Consequences of Aviation Safety Occurrences,”

**Safety Science**, Vol. 111, Jan. 2019, pp. 243–252. doi:https://doi.org/10.1016/j.ssci.2018.07.018[39] , “Mixtures of Gamma Distributions with Applications,”

**Journal of Computational and Graphical Statistics**, Vol. 10, No. 3, 2001, pp. 440–454. doi:https://doi.org/10.1198/106186001317115054[40] [B2] “FSF ALAR Briefing Note 8.5—Wet or Contaminated Runways,” Flight Safety Digest, Aug.–Nov. 2000, pp. 179–184, http://flightsafety.org/files/alar_bn8-5-wetrwy.pdf [retrieved 2019].

[41] , “Bayesian Inference for the Beta-Binomial Distribution via Polynomial Expansions,”

**Journal of Computational and Graphical Statistics**, Vol. 11, No. 1, 2002, pp. 202–207. doi:https://doi.org/10.1198/106186002317375686[42] ,

**Bayesian Data Analysis**, 3rd ed., Chapman and Hall/CRC, Boca Raton, FL, 2013.[43] , “Bayesian and Non-Bayesian Analysis of the Regression Model with Multivariate Student-t Error Terms,”

**Journal of the American Statistical Association**, Vol. 71, No. 354, 1976, pp. 400–405. doi:https://doi.org/10.1080/01621459.1976.10480357

## Tables

Guideline | Description |
---|---|

a | Fly a stabilized approach. |

b | Height at threshold crossing should be at most 50 ft.a |

c | Speed at threshold crossing to be within 10 kt of the reference airspeed ${V}_{\mathrm{ref}}$. |

d | Tailwind should be no more than 10 kt (00 kt) for noncontaminated (contaminated) runways. |

e | Touch down on runway centerline at touchdown point. |

f | After touchdown, promptly transition to the desired deceleration setting. |

g | Speed should be less than 80 kt with 2000 ft of landing distance remaining. |

Variable | Definition |
---|---|

$A$ | Runway. Categorical: identifies the relevant runway. |

$B$ | Crosswind (knots). Nonnegative continuous, discretized to nearest integer. |

$C$ | Tailwind (knots). Same as crosswind. |

$D$ | Runway condition. Categorical: contaminated or not.a |

$E$ | Type of approach. Categorical: stabilized or not. |

$F$ | Maximum reverse thrust (seconds). Continuous. Time of MRT operation during ground roll.b |

$G$ | Autobrake. Categorical: no autobrake, low and medium. |

$H$ | Height at threshold (feett). Nonnegative continuous; equals zero when aircraft touches down. |

$I$ | $\mathrm{IAS}-{V}_{\mathrm{app}}$ (knots). Continuous. Difference between IAS and ${V}_{\mathrm{app}}$ at threshold.c |

$Y$ | LDR at 80 kt (feet). |

For $k=1$ To $N$ | |

1: | Sample ${a}_{k}\sim p(a)$ |

2: | Sample ${b}_{k}\sim p(b|{a}_{k})$, ${c}_{k}\sim p(c|{a}_{k})$, ${d}_{k}\sim p(d|{a}_{k})$, ${g}_{k}\sim p(g|{a}_{k})$ |

3: | Sample ${e}_{k}\sim p(e|{a}_{k},{b}_{k},{c}_{k})$ |

4: | Sample ${f}_{k}\sim p(f|{a}_{k},{e}_{k})$, ${h}_{k}\sim p(h|{a}_{k},{e}_{k})$, ${i}_{k}\sim p(i|{e}_{k})$ |

5: | Sample ${y}_{k}\sim p(y|{a}_{k},{d}_{k},{f}_{k},{g}_{k},{h}_{k},{i}_{k})$ |

6: | $k=k+1$ |

Data | |
---|---|

Rwy-01 | 3.27 (0.56) |

Rwy-02 | 2.47 (0.51) |

Rwy-03 | 0.29 (0.84) |

Rwy-01 | Rwy-02 | Rwy-03 | |||||||
---|---|---|---|---|---|---|---|---|---|

$\underline{b}=10$ | $\underline{b}=15$ | $\underline{b}=20$ | $\underline{b}=10$ | $\underline{b}=15$ | $\underset{\_}{b}=20$ | $\underline{b}=10$ | $\underline{b}=15$ | $\underline{b}=20$ | |

$\underset{\_}{c}=0$ | 7.81 | 9.72 | 11.08 | 1.79 | 1.94 | 2.14 | 0.21 | 0.34 | 0.49 |

$\underset{\_}{c}=5$ | 9.43 | 10.57 | 12.11 | 1.98 | 2.17 | 2.31 | 0.26 | 0.39 | 0.48 |

$\underset{\_}{c}=10$ | 11.84 | 12.81 | 13.44 | 2.77 | 3.95 | 4.72 | 0.60 | 0.85 | 0.99 |

$\underset{\_}{f}=0$ | $\underline{f}=5$ | $\underset{\_}{f}=10$ |
---|---|---|

Rwy-01 | ||

0.36 | 1.31 | 1.72 |

Rwy-02 | ||

0.18 | 0.54 | 0.84 |

Rwy-03 | ||

0.0018 | 0.0077 | 0.016 |

$w\sim \mathrm{Dir}(\varphi )$ | $\nu \sim \mathcal{E}(\theta )$ | $\mu \sim \mathcal{IG}(\alpha ,\beta )$ |
---|---|---|

$\varphi $ | $\theta $ | $(\alpha ,\beta )$ |

$(0.4,0.4,0.2)$ | 0.1 | $((10,10,10),(\mathrm{200,400,600}))$ |

$(0.5,0.4,0.1)$ | 0.2 | $((2.2,2.2,2.2),(\mathrm{40,50,100}))$ |

Rwy-01 | Rwy-02 | Rwy-03 | |
---|---|---|---|

${\nu}_{\ell}^{b}|\mathrm{data}$ | 3.76 (0.50) | 3.96 (0.60) | 6.11 (0.91) |

${\mu}_{\ell}^{b}|\mathrm{data}$ | 6.50 (0.34) | 5.21 (0.29) | 4.03 (0.18) |

Tailwind | Rwy-01 | Rwy-02 | Rwy-03 |
---|---|---|---|

No | 49 | 29 | 47 |

Yes | 55 | 51 | 35 |

Rwy-01 | Rwy-02 | Rwy-03 | |
---|---|---|---|

${\nu}_{\ell}^{c}|\mathrm{data}$ | 2.29 (0.40) | 3.16 (0.58) | 2.78 (0.62) |

${\mu}_{\ell}^{c}|\mathrm{data}$ | 4.49 (0.41) | 4.34 (0.35) | 4.24 (0.44) |

Runwwy conditions | Rwy-01 | Rwy-02 | Rwy-03 |
---|---|---|---|

Noncontaminated (0) | 0.89 (93) | 0.93 (75) | 0.89 (74) |

Contaminated (1) | 0.11 (11) | 0.07 (5) | 0.11 (8) |

Rwy-01 | Rwy-02 | Rwy-03 | |
---|---|---|---|

${\alpha}^{e}|\mathrm{data}$ | $-10.96$ (3.81) | $-8.33$ (3.02) | $-8.47$ (2.64) |

${\beta}_{1}^{e}|\mathrm{data}$ | 0.54 (0.26) | 0.49 (0.29) | 0.59 (0.33) |

${\beta}_{2}^{e}|\mathrm{data}$ | 0.39 (0.23) | 0.11 (0.37) | 0.48 (0.23) |

RT | ||||||
---|---|---|---|---|---|---|

Stabilized | Unstabilized | |||||

Rwy-01 | Rwy-02 | Rwy-03 | Rwy-01 | Rwy-02 | Rwy-03 | |

Idle | 39 | 41 | 37 | 0 | 0 | 0 |

Maximum | 63 | 38 | 42 | 2 | 1 | 3 |

Rwy-01 | Rwy-02 | Rwy-03 | |
---|---|---|---|

${\nu}_{\ell}^{f,s}|\mathrm{data}$ | 5.15 (0.88) | 6.16 (1.35) | 3.35 (0.70) |

${\mu}_{\ell}^{f,s}|\mathrm{data}$ | 9.47 (0.52) | 10.24 (0.70) | 9.45 (0.83) |

Rwy-01 | Rwy-02 | Rwy-03 | |
---|---|---|---|

${\nu}_{\ell}^{f,u}|\mathrm{data}$ | 5.27 (0.88) | 3.95 (0.84) | 3.34 (0.66) |

${\mu}_{\ell}^{f,u}|\mathrm{data}$ | 9.46 (0.52) | 9.98 (0.82) | 9.39 (0.79) |

Autobrake | Rwy-01 | Rwy-02 | Rwy-03 |
---|---|---|---|

No (1) | 0.21 (22) | 0.38 (31) | 0.51 (42) |

Low (2) | 0.12 (12) | 0.28 (22) | 0.27 (22) |

Medium (3) | 0.67 (70) | 0.34 (27) | 0.22 (18) |

Component 1 | Component 2 | Component 3 | |||||||
---|---|---|---|---|---|---|---|---|---|

Runway | ${\widehat{{\nu}_{1}}}^{h,s}$ | ${\widehat{{\mu}_{1}}}^{h,s}$ | ${\widehat{{w}_{1}}}^{h,s}$ | ${\widehat{{\nu}_{2}}}^{h,s}$ | ${\widehat{{\mu}_{2}}}^{h,s}$ | ${\widehat{{w}_{2}}}^{h,s}$ | ${\widehat{{\nu}_{3}}}^{h,s}$ | ${\widehat{{\mu}_{3}}}^{h,s}$ | ${\widehat{{w}_{3}}}^{h,s}$ |

Rwy-01 | 2.00 (0.55) | 9.98 (1.65) | 0.33 (0.07) | 14.03 (11.70) | 40.13 (10.49) | 0.27 (0.19) | 13.19 (10.65) | 63.37 (12.83) | 0.40 (0.20) |

Rwy-02 | 1.74 (0.40) | 9.68 (1.42) | 0.51 (0.09) | 23.71 (18.54) | 26.96 (8.22) | 0.19 (0.14) | 27.40 (14.29) | 44.45 (12.37) | 0.30 (0.16) |

Rwy-03 | 1.76 (0.48) | 10.54 (2.13) | 0.40 (0.12) | 9.57 (10.80) | 33.38 (9.04) | 0.32 (0.19) | 16.24 (20.40) | 68.64 (18.12) | 0.28 (0.19) |

Component 1 | Component 2 | |||||
---|---|---|---|---|---|---|

Runway | ${\widehat{{\nu}_{1}}}^{h,u}$ | ${\widehat{{\mu}_{1}}}^{h,u}$ | ${\widehat{{w}_{1}}}^{h,u}$ | ${\widehat{{\nu}_{2}}}^{h,u}$ | ${\widehat{{\mu}_{2}}}^{h,u}$ | ${\widehat{{w}_{2}}}^{h,u}$ |

Rwy-01 | 14.72 (11.23) | 39.77 (10.21) | 0.67 (0.23) | 12.88 (10.33) | 63.89 (13.21) | 0.33 (0.21) |

Rwy-02 | 24.11 (18.21) | 27.12 (8.09) | 0.71 (0.14) | 27.76 (14.93) | 43.94 (11.99) | 0.29 (0.16) |

Rwy-03 | 9.03 (9.89) | 33.82 (9.17) | 0.80 (0.23) | 16.66 (21.07) | 66.94 (19.25) | 0.20 (0.16) |

Approach | ${n}_{e}$ | ${\overline{i}}_{e}$ | ${s}_{e}^{2}$ |
---|---|---|---|

Stabilized | 260 | 2.47 | 14.97 |

Unstabilized | 6 | 5 | 20.80 |

Rwy-01 | Rwy-02 | Rwy-03 | ||
---|---|---|---|---|

${\alpha}^{y}|\mathrm{data}$ | 3107.63 (53.48) | 3655.32 (84.21) | 3553.83 (56.93) | |

$D$ | ${\beta}_{1}^{y}|\mathrm{data}$ | 10.04 (4.32) | 110.38 (71.65) | 91.68 (57.08) |

$F$ | ${\beta}_{2}^{y}|\mathrm{data}$ | 5.25 (2.41) | 3.65 (1.24) | 4.99 (3.17) |

${G}^{l}$ | ${\beta}_{3}^{y,l}|\mathrm{data}$ | 125.46 (21.04) | 186.97 (28.57) | 171.28 (37.80) |

${G}^{m}$ | ${\beta}_{3}^{y,m}|\mathrm{data}$ | 428.88 (70.88) | 518.36 (81.12) | 556.88 (94.58) |

$H$ | ${\beta}_{4}^{y}|\mathrm{data}$ | 3.72 (1.31) | 8.74 (2.33) | 0.84 (0.56) |

$I$ | ${\beta}_{5}^{y}|\mathrm{data}$ | $-25.34\text{\hspace{0.17em}}(12.79)$ | $-24.11\text{\hspace{0.17em}}(11.99)$ | $-23.16\text{\hspace{0.17em}}(11.06)$ |