3. Exploratory data analysis for failure data

Before fitting any life distribution or computing reliability metrics, you must understand the shape of the failure data you actually collected. Exploratory data analysis (EDA) for time-to-failure (TTF) datasets answers three questions: What does the distribution look like? Where is the central tendency and spread? Which observations are anomalous and deserve investigation? Skipping this step is the most common cause of fitting a Weibull model to data that is actually a mixture of two failure modes, or of reporting a mean time to failure (MTTF) that is dominated by a single misrecorded timestamp.

Why Failure Data Needs Special Treatment

TTF data is inherently non-negative, usually right-skewed, and frequently censored. A bearing that has not failed after 4,000 hours contributes a right-censored observation, not a failure at 4,000 hours. Standard descriptive statistics such as the arithmetic mean assume symmetry and no censoring; applied naively to failure data they mislead. EDA must therefore be paired with a clear record of which observations are failures (event = 1) and which are suspensions (event = 0).

Consider a dataset from an accelerated life test of 40 identical DC motor brushes, run at elevated temperature. Each unit records time_hours and event. Twenty-eight units failed; twelve were still running when the test was suspended. The raw failure times (hours) are:

112, 145, 158, 171, 189, 203, 214, 226, 238, 251,
263, 277, 289, 302, 318, 334, 351, 369, 388, 410,
433, 458, 486, 517, 552, 592, 638, 691

The twelve suspensions occurred at 720, 720, 720, 720, 720, 720, 720, 720, 720, 720, 720, 720 hours (test termination). Note that all suspensions share the same value — a common pattern in fixed-duration tests.

Histograms: Seeing the Shape

A histogram bins the failure times and plots frequency per bin. For the brush data, choose bin width using a rule such as Sturges' formula: k = ceil(log2(n)) + 1. With n = 28 failures, k ≈ 6. A bin width of 100 hours gives bins [100–200), [200–300), …, [600–700). The resulting counts are 3, 6, 7, 5, 4, 2, 1. The histogram shows a pronounced right skew: most failures cluster between 150 and 400 hours, with a long tail extending past 600 hours.

Right skew is the signature of wear-out or fatigue-dominated failure. If the histogram were roughly symmetric, you would suspect a random (exponential-like) failure mode. If it were bimodal — two distinct humps — you would suspect a mixture of two populations, for example infant mortality plus wear-out, and you should split the data by failure mode before fitting.

Boxplots: Quantiles, Spread, and Outlier Candidates

A boxplot summarizes the distribution using five order statistics: minimum, first quartile (Q1), median (Q2), third quartile (Q3), and maximum. For the 28 failure times, sorted, the median is the average of the 14th and 15th values: (302 + 318)/2 = 310 hours. Q1 is the 7th value (214 hours) and Q3 is the 22nd value (486 hours). The interquartile range (IQR) is Q3 − Q1 = 272 hours.

The standard outlier rule flags any observation below Q1 − 1.5·IQR or above Q3 + 1.5·IQR. Here the upper fence is 486 + 1.5 × 272 = 894 hours. No failure time exceeds 894, so the boxplot shows no outliers among the failures. The lower fence is 214 − 408 = −194, which is impossible for TTF data, confirming that the lower tail is bounded by zero.

This is an important lesson: the 1.5·IQR rule is a screening heuristic, not a verdict. A point beyond the fence is an outlier candidate that must be traced to its source — a data-entry error, a unit from a different production lot, a test-rig anomaly, or a genuine rare failure mode. Deleting it without investigation biases the reliability estimate.

Identifying Outliers in Time-to-Failure Data

Outlier detection in TTF data differs from generic tabular outlier detection in three ways. First, the distribution is skewed, so distance-based methods that assume normality (z-scores) are inappropriate. Second, censoring means the largest values may be suspensions, which are not outliers at all. Third, a single extreme failure time can legitimately belong to a different failure mechanism.

A practical workflow:

1. Separate failures from suspensions. Compute quartiles using failures only.
2. Apply the 1.5·IQR fence to failures. Mark candidates.
3. For each candidate, inspect the raw test log: timestamp, unit ID, test cell, operator, and failure mode code.
4. Cross-check against a probability plot. On a Weibull probability plot, a true outlier departs from the fitted line at the extreme; a mixture shows curvature in the middle, not just the tail.
5. Document the disposition: retain with justification, correct the record, or exclude with a stated reason.

Suppose a 29th failure is later found in the log at 1,240 hours, unit M-17, with failure mode "brush fracture" rather than the normal "brush wear." The 1.5·IQR fence flags it (1,240 > 894). Investigation shows a handling crack introduced during installation. This is a legitimate exclusion from the wear-out population, but it must be recorded and analyzed separately as a handling-induced failure.

Reading the Combined View

The histogram and boxplot are complementary. The histogram shows modality and skew; the boxplot shows quantiles and flags candidates. Plot them together, and annotate the sample size, number of failures, number of suspensions, and the fence values. The SVG below summarizes this EDA pipeline for the brush dataset, from raw records through distribution shape assessment to outlier disposition.

Practical Notes on Tooling

In Python, the workflow is compact. Load the data with pandas, split on the event column, and compute quartiles on the failure subset:

import pandas as pd
import numpy as np

df = pd.read_csv("brush_ttf.csv")   # columns: unit_id, time_hours, event
fail = df[df.event == 1]["time_hours"]

q1, q2, q3 = fail.quantile([0.25, 0.50, 0.75])
iqr = q3 - q1
lower, upper = q1 - 1.5*iqr, q3 + 1.5*iqr
candidates = fail[(fail < lower) | (fail > upper)]
print(f"Q1={q1}, Q2={q2}, Q3={q3}, IQR={iqr}")
print(f"fences: {lower:.1f} / {upper:.1f}")
print(candidates)

For the brush data this prints Q1=214.0, Q2=310.0, Q3=486.0, IQR=272.0 and fences -194.0 / 894.0, with an empty candidate list until M-17 is added. Plotting with matplotlib or seaborn produces the histogram and boxplot shown above.

What to Carry Forward

EDA is not a formality before modeling; it determines which model is defensible. From the brush example you now have: a right-skewed failure-time distribution consistent with wear-out, quartile estimates (Q1 = 214 h, median = 310 h, Q3 = 486 h), an upper fence of 894 h, and one documented outlier (M-17) excluded on physical grounds. These numbers feed directly into the next step — selecting and fitting a life distribution — and they give you a sanity check: if a fitted Weibull predicts a median far from 310 hours, something in the fit or the data preparation is wrong.