The problem #

The Backblaze company is a well-known cloud-based storage company, who publicizes fleet-related data every quarter. In the following blog entry, we will be using different methods to analyze, illustrate, model, verify and explain the data in a realistic setup.

The data #

Before doing anything, we know very well that understanding the data is a half way to success. We will thus spend a good part to properly understand the data - what happens and why that happens. The actual data is available on Backblaze’s site. The data is downloaded as CSV’s for every quarter. Every quarterly data yields a ~10GB CSV file. This already gives the approximate understanding of the weight of the data we’ll be working width. The data we will be working with is very heavy and will often become our bottleneck in our analyses. It is the exact reason we’ll fail to work with the full dataset (many years) at once. The data has the form of data_Q<N>_<YYYY>.csv, whith N is the quarter number 1,2,3,4 and YYYY - the year. In order to work with the data in the most efficient way, we will be converting the CSV’s into parquet files. The columns of the data consists of the so-called main info, indicating the date, the location, the datacenter, … and the so-called SMART values, standing for Self-Monitoring, Analysis and Reporting Technology. These are values reported by the hardware to monitor and indicate different quantities/values. We note that the values are often not absolute, that is, there is no absolute scale, specifically throughout different manufacturers. This thus means we must interpret them with a grain of salt. The complete list of the columns is as follows:

[ "date", "serial_number", "model", "capacity_bytes", 
"failure", "datacenter", "cluster_id", 
"vault_id", "pod_id", "pod_slot_num", "is_legacy_format", ..."smart_<smart_id>"]

EDA #

The dataset #

We will be using polars to work with the tabular dataframes. In particular, it’s lazy execution functionality, that works great with such large datasets.

Let’s inspect the format of the data we’re working with

The original/raw data is on the daily level - each entry is a daily per unit event that logs the status of the drive. A drive can be healthy and the entry is logged - this is the default/normal behavior. The second option is the failure - denoted by a flag, that indicates that the drive has failed on this day. Another option is that the entry for a certain unit is simply absent. This can be due to the maintenance, preventive actions or a another type of logging issue (all according to backblaze). This issue may or not induce a full censoring.

Intermediary DF’s #

The goal of the first phase of the data exploration is to understand and “standardize” different cases of unusual behavior of logs. To do that, let us roughly define the vocab first. We will call a unit to be one hard drive defined by a unique serial number. We will define the spell or episode to be a continious logged period of time. A gap is the time between two consecutive spells. For example, if we observe 2 distinct periods of continous period of existing entries separated by a gap of two days, we say that the two spells/episodes (each spell within a unit is numerated) is separated by a two-day gap.

We will create intermediary dataframes containing information about spells and units, as well as a final “clean” dataframe, which will be decisive when working with statistical models and estimators.

Let’s introduce the main dataframes. The first is the spell dataframe with the following spells.collect().head():

where each row instead of being a unit-day, it is a spell per unit. We added some new informative columns - the missing_days, indicating, how much days without log days are missing. We note that we have defined some tolerance GAP_TOLERANCE, under which the gap is not considered to be a gap (e.g. in a streak of Mon Tue Wed Fri, we consider Thu to have a log entry).

The second dataframe units has introduced multiple columns

namely, flaky - indicating whether a unit has had many short gaps (the many has a hard-coded threshold), interrupted - whether a unit has had too many observation spells, resurrected - whether a unit has failed and then reappeared in the time window. These features will be the ones determining how “clean”/suited they are for using them in later estimators. We can show

The way we define the thresholds and the allowed gaps are fully arbitrary - we create them based on our interpretation of the process.

The SMART features #

The last “feature” to inspect is the availability of the SMART features, namely, how “available” is a feature per model (as the SMART features are generally manufacturer/model dependent). We can illustrate it as a heatmap (see the expand option below).

As a final exploration step, we can look at the distribution of the fleet by models as a function of time to get a more complete picture. We will keep the first/most common $k$ models to simplify the visualization and interpretation.

First evaluations - rates & Kaplan-Meier #

The very first quantity to compute is to aggregate and compute the (annualized) failure rate, which we define as

Without doing any statistical tests, we see that there is heterogeneity across different models. Similarly, we can perform a similar thing across other units, e.g. data centers, to confirm statistically significant differences to eventually include it in the later feature set.

Kaplan-Meier & survival analysis #

The central quantity to estimate is the survival function $S(t)$ The central quantity in many survival-related tasks is the survival function or any derived quantities. The survival function $S(t)$ is the probability that the unit survives past time $t$. In a slightly more formal setting, $S(t)\coloneqq 1 - F(t)$, where $F(t)$ - the cumulative distribution of $f(t)$, which is the pdf of $f(t)$ - having a failure at time $t$. There are other fundamental quantities, such as the hazard rate $h(t)$ (often denoted $\lambda$), which expresses the rate failure at time $t$. Note that this is not a probability function and does not have the probability properties. Another common quantity is the cumulative hazard $H(t)$, which is the cumulative hazards at up to $t$. The quantities are related via known relations -

$$ S(t) = \exp( - \int_0^t h(s)ds )\;, \qquad S(t) = 1 - F(t)\; \qquad H(t) = \int_0^t h(s)ds$$

The most straightforward way of estimating $S(t)$ is the Kaplan-Meier estimator, which is defined as

where $d_j$ - number of events at $t_j$ and $Y_j = |R(t_j)|$ - the size of the risk set just before $t_j$. Instead of manually implementing the algorithm, we can use the survival analysis tollbox lifelines, which implements many of the known functionality out of the box.

We can attempt to understand the formula of the KM estimator. We will not be a rigorous path - instead, a very handwavy one. We will start by the Nelson-Aalen estimator, which is the estimator of the cumulative hazard rate, given by $\hat{H}(t) = \sum_{s_i\leq t} d_i/n_i$ - a much more intuitive quantity to interpret. Indeed, from the definition of the hazard rate, $H(t) \coloneqq \int_0^t h(s)ds$, which, in the discrete case takes the form $H(t) \simeq \sum_{s\leq t} (\Delta s) h(s)$ in the Riemann sum sense. So the quantity $\frac{\text{number of events}}{\text{number of observed units}}$ is exactly the quasi-probability quantity $\text{p} \simeq h(t)dt$ in the infinitesimal sense! We thus arrived to the (informal) equivalence of $d_i/n_i$ and $\Delta s h(s)$. From that, we write

For the KM, we will be using the dataframe by spells. Actually the difference between using the by-spells df and by-units is almost none and arises only from our assumptions of allowed gaps and bridging. There are, nevertheless, two ways time scales that we construct our KM estimator with. We can use the observed dates - the start and end date of the study. Meaning the only information is the entry and exit times of the study (e.g. 01.01.2022 to 31.12.2025). The second scale is to use the power-on hours of each hard drive. This means a hard drive entering the study (on the 01.01.2022) will already have some “story”, i.e. some elapsed hours e.g. 20k and some exit time. Thus, in the second case, we have a bit more context. Indeed, in the first case, an old drive surviving the whole study period and a recent one will look the same for the model. Let’s compare the two.

As expected, the time scales are different, and the overall appearance of the estimates are similar. The most apparent difference is the size of the confidence interval, which is much wider in the hour’s case. Let’s attempt to understand the intuition behind that. When we’re working with the annual calendar scale, all of the units start together and evolve together. If a failure occurs at some hour $H$, this makes a contribution/gives info to the whole fleet having the age of $H$ (which we have a lot of). If, however, we use the “birth” hours, the events are much more sparse. So when a drive fails at some $H'$, there are actually very few fails at that period of time. In other words, the information about the (life)time is much more scattered in this case, than in the previous one, explaining the difference in the widths of CI’s.

Aalen-Johanses & competing events #

Up to now, we’ve considered an exlusive failure event - either a hard drive (unit) died or not. A very common case are multiple failure modes. Namely, there are multiple ways a hard drive may “die”. Although this is a very natural scenario in many cases (a different failure cause), it is less so in this situation. Yet, there are cases of decomissions of hard drives, where many hard within a same location (vault) are removed - decomissioned, for example, 1.2k at once. This is nothing but a heuristic and may or not become come up to be a reasonable one!

In order to consider the additional failure mode, the simple KM estimator is not suitable anymore. A natural extension to it is the Aalen-Johnson estimator. It estimates the so-called Cumulative Incidence Function (CIF), which is a modified cumulative function that gives the probability $F(t) = P(T < t, J=j)$ where $J=j$ symbolically represents the failure mode. The analytical expression of this cumulative (incidence) function is given by

,where the $S(u^-)$ represents the survival function (corresponding to multiple failure modes) just before $u$. In the analytical case, it is less relevant (left-hand-side limit) than in the discrete case, where it corresponds to the previous time $t_{i-1}$. We can easily show that this (in the case of a single failure) is the KM estimator. Indeed, remember that $S(t) = e^{\int^t h(u)du} = \text{KM}(t)$ meaning that $S'(t)=-h(t)S(t)$

We can compare the curves - of the AJ estimator and KM and get a feeling of how they relate to each other:

We see that the survival curve of the migration event drops roughly at around 2.5k hours, then decreases slowly, compared to the real failure mode. We could approach this rigorously, but here, we will accept (or assume), that the migration failure mode is not significant enough to model. This is fully arbitrary. Also, this could very weell be, because our assumptions were not fully valid, or are not fully valid within our limited dataset span!

Imposing parametric forms #

The next step is to involve some sort of parameters. That is, assumptions about the shape of the distributions, controlled by some parameters. The most straightforward model is the Accelerated Failure Time (AFT). This framework treats the time-to-failure to be “variable” under different covariates. Instead of the well-known Cox regression, where the hazard rate $h(t)$ changes depending on covariates, AFT stretches the time scale. Let $T_i$ - the survival time of the item $i$ r.v. we are modeling. We write it as

which exponentiating gives

with $\beta$ - the regression coefficients, $x$ - the covariates, $\epsilon$ - the deviation/noise of the model and $\sigma$ - the scale of this noise. The random error is what we specify by the model. It enters in $T_{0,i}=\exp(\sigma \epsilon)$, which is the baseline, which we assume to be true before taking into account the covariates. The most common assumptions about $\epsilon$ are Weibull, Lognormal and Log-logistic distributions. In our data, however, it turns out that these “classical” assumptions fail to catch some of the behavior. Some models have low failure rates and the Weibull estimates are not very trustworthy. For example, for a flat-like Weibull estimate, the optimizer decreases $\rho$ and sends $\lambda$ to $\infty$. As a result, we have decided here to employ the Generalized Gamma (GG) distribution

Which is a generalization of the known Weibull, Exponential, Gamma, etc… The former ones are thus sub-models of the GG. We can have a look at how our non-parametric KM compares with the parametric GG:

Altough the curves do not perfectly coincide, the shapes are very well captured using our GG fits. For many cases, the curves are overlapping, which confirms that the model is actually capturing well the decrease in $S(t)$. For the first time, we can actually predict the probability of failures during the lifetime of the hard drive. From there, we can study the quality of the potential predictions. In this specific case study, we’re more interested in describing the data using different methods, not a classification/regression task. This is what leads us into the Bayesian framework!

Applying Bayesian framework #

The dataset has a clear nested structure. Multiple drive models, belonging to a certain manufacturer, sitting in a datacenter, in a pod, etc… Suppose we want to study different models and manufacturers age and fail over time. We will be using a multilevel formulation with a parametric Weibull degradation model.

The Weibull model here is taken solely due to its wide use cases. In order to do a more formal justification, there are plenty of things to check and potentially many other models to compare it against. This, however is slightly out of the scope of this small blog post and we will here simply accept the Weibull model to some degree.

Likelihood #

Before specifying the anatomy of the hierarchy, we will first describe the likelihood. A somewhat similar, yet simpler logic has been applied in one of the previous articles.

The likelihood or in Bayesian notation is given by $\mathbb{P}(\vec{D}| \vec{\Theta})$ - the probability of the observed data given the parameters. Let’s consider two cases. We define the spell to be an episode “fail” event or the censoring event¹.

Case 1: The unit enters in the observed window. The pdf we consider for the r.v. $T$ is $f$. The episode enters in the observation window only if $T>t_0$ some starting time. Thus the pdf $f$ of the episodes under consideration must be conditioned on $T>t_0$. We thus write

Thus, if the unit fails at time $t$, the likelihood is given by

Case 2: The unit enters in the observed window, but no event is happening. In this case, the event instead of $t \in (t,t+\Delta t)$ becomes simply $t < T$ yielding

Thus, if no fail is seen, the likelihood becomes

This if-else structure can be encoded in a indicator $\delta \in \{ 0, 1 \}$, with $0$ representing right-censoring giving the full likelihood to be ²

We can customly implement this as either a constant sampling term (pm.Potential) or a custom distribution. In any of these cases, we will work in the logspace, so $L(t)$ becomes log-likelihood $\ell(t)$.

As we’ve mentioned, in this Bayesian study, we’re using a Weibull parametrization, meaning a Weibull

and

where $\lambda,k \in (0,+\infty)$ - scale and shape respectively.

Hierarchy #

Let’s look at how we construct the hierarchical model. We denote the model $m$, for which we define $a_m = \log(\lambda_m)$ and $b_m = \log(k_m)$ - the two Weibull per-model parameters.

For $a_m$, we define

with $j$ - the manufacturer’s index. So $\mu_j$ - the model’s grand mean, tight to the manufacturer, from which models are allowed to deviate. $u_m$, $\tau$ determine how far they’re allowed to deviate.

Similarly, for $b_m$, we have

with $\rho v_m$ playing the same role as $\tau u_m$ - how much the per-model parameter can deviate from the grand mean, in this case, $\nu$. Note that there’s is no manufacturer’s level sharing. The reason is the actual meaning of the $k$ and $\lambda$. $\lambda$ - the scale parameter represents how fast/late the object fails ($S$ drops at early/late $t$’s - see interactive marimo below), which is lowkey what we want to estimate.

The shape, however, is an intrinsic feature of a wear-out process. Thus, there are no particular reasons to make it vary per model and even manufacturer - we let it fluctuate it around the value, representing the pooling in the multilevel model.

Visualizing #

Once the model, the hierarchy and the priors are specified, one can sample. The issue, however, is the complexity of the geometry, and, unfortunately, the size of the dataset needed to be fitted is not the main problem. The main problem encountered is the geometry that is not easy to sample from, confirmed with very high $\hat{R}$ values.

Since I found it pretty interesting and challenging (one of the reasons this took me some time and why I left many further concepts unimplemented), we will be discussing this in separate notes.

Let us assume we now sampled successfully, confirmed by all the main metrics, and let’s look at one of the parameter levels - $a_m$ and $b_m$ for Toshiba models.

We remember the definition of $a_m$ and $b_m$, which are related to $\lambda_m, k_m$ via the exponential. The observed $\lambda_m$’s interquantile range ranges $[8.3, 10.85]$, which corresponds to $10-90$ years of lifetime. It is true, that the values look overestimated. This, however, are just point estimates, and do not represent all of the estimated models.

Yet we can still attempt to understand this (over)estimation. Possible first-order reasons include:

1. Incorrect definitions of the model - spells which have died and the “resurrected” does not correspond to an actual event. To mitigate this, one needs a better understanding of the lifetimes’ processes.
2. Increase the data window
3. Reconsider the priors and the overall hierarchy of the model - adding/removing levels
4. Reconsider the base (Weibull) model

Continuing #

Recalling - the goal is to show a proof-of-concept and not scoring by a full model selection. The next (possibly post-hoc) steps would be to perform a prior-predictive check to reconsider the priors, levels, Model. Additionally, to include new models and compare fit between each other. Based on quantitative and qualitative results, we may continue to dig into this direction or completely change of direction by including other notions - fraility, Cox, … or even stepping outside of the Bayesian framework.

Hopefully, to be continued…

Concluding #

This small entry shows that (hard-drive) survival analysis cannot really be treated as a simple table of independent lifetimes. The raw backblaze data is messy with hidden operations, since they keep disappearing and re-appearing even when flagged as failed. Therefore, constructing a working dataset requires building a usable survival dataset and performing assumptions, that may break afterwards.

Within the defined spells/episodes, the failure behavior is clearly heterogeneous at both model and manufacturer levels. Although AFR’s is a good first start and even presents uncertainty, qualitative conclusions with reliable uncertainties are able to be drawn only at the Kaplan-Meier curves. The survival curves with the generalized gamma model turned out to be highly quantitatively useful.

Within any use approach, the models/manufacturers (levels) with few failures remain hardly identifiable. The quality, however improves by adding data, making stronger assumptions or chooosing other models (where these assumptions are implicitly made).

The main lesson could be simple: survival analysis does not boil down to the choice of the estimator - whether it is KM, NA etc.. As within any inferential and/or statistical modelling problem, the quality and credibility depends on the operational definitions of main quantities and interpretation. If those are available, one can safely build more complex, possibly hierarchical, bayesian or any other type of statistical models without having negative impact of the data on the model.