Poisson, Geometric & Survival Models: A Hazard-Based Look at Persistent Failure in Retesting

The Problem #

Quality control and testing actions are always mandatory before delivering a product or a service. In this particular case study, we are considering a some manufacturing data recording during testing. Every row is the result of a testing procedure. The tabular data may look like this

timestamp	id	result	retests
13:32:10	#AA3	pass	1
13:32:20	#AA4	failmode1	1
13:32:30	#AA5	pass	1
13:32:50	#AA4	failmode1	2
13:33:20	#AA4	pass	3

From this tabular data, we may look how many units passed on first try, how much units failed with mode 2 on first, second, … tries etc… When a unit fails some $K$ times, it is scrapped. In our demo, we take $K$ to be 4. Different failure modes depend on a various factors. We may divide this into two main “categories” - the testing procedure or the intrinsic problem of the unit. The examples below use synthetic data generated “by hand” with the same structure as the testing data I want to discuss.

The dataset may be modeled in various, sometimes extremely complicated ways, which we will keep for later discussions, to model both the temporal and spacial dependences. For now, we focus on simple summary statistics and fit different types of models, since sometimes, we all need simple things! Another reason of using a summary table are the significant costs related to work/retrieve the full data.

The goal is to describe the summary statistics using some models and observe the severity of some failure modes. Namely, if some failure mode is intrinsic to the batch/units due to some manufacturing issues, many units will be failing with this mode without “recovering”. We will thus want to somehow quantify this issue.

The summary table will thus look like this

failure_mode	rank	count
SUCCESS	2	29769
FMODE_C	2	351
FMODE_C	3	201
FMODE_C	1	1038

We can illustrate the difference between a “real” problem and a “random” one.

First-order estimation #

As we are working with a short summary table, it is clear we are losing a lot of information and relationships. Indeed, we dropped all of the spatio-temporal dependence. Our first central assumption consists of considering that there exists a unique probability for each failure mode. Namely,

$$ \mathbb{P}(\text{fail} \in (M,r) | \text{fail} \in (M, r-1)) = \mathbb{P}(\text{fail} \in (M,r)) $$

, which is almost always not true due to the persistence of a failure mode, due to the temporal dependence and many other reasons. Yet this assumption is a good starting point, especially for discovering and learning about the processes.

Our assumption assumes a model, where the number of units getting $r$-retested will decay geometrically. Let

$$ C_{M,r} $$

the number of testing counts at failure mode $M$ and rank $r$. Our first-order decay model gives

$$ C_{M,r} \sim C_{M, r-1}q_{M} $$

Where $q_M$ is what we will call here the persistence factor of the mode $M$. This characterizes how good the tested units are “recovered” after failures. We may write for any testing attempt

$$ C_{M,r} \sim C_{M,1} q_{M}^{r-1} $$

A small persistence would mean that after the first retest, there will be few units to be retested more, so the decay will be quick.

The first estimate will therefore be by simply computing the ratio

$$ \hat{q} \coloneqq (\frac{C_{M,K}}{C_{M,1}}) ^{1/(K-1)} $$

which is a great first estimate for what we need - to quantify the persistence effect of a certain mode.

Poisson model #

The Poisson model is a well-known model that models events or count data. This is the reason we will be interested in using it to model our situation. In fact, in our setup, there are important notions of success/failure, number of re-tests, … Intuitively, this signals the usage of Bernoulli, geometric and other distributions. In the first order, however, we will consider the retest ranks simply as binned/histogram events.

Remember our discrete decaying model. In this model, we expect the number of retests $\mu$ to behave as

$$\mu_{M,r} \sim A_{M} q_{M}^{r-1}$$

It is highly tempting to use some $\log$ transformation here, so taking it from both sides and developing gives

$$ \log[ \mu_{M,r}] = \log[ A_{M} q_{M}^{r-1} ] \\ \log[ \mu_{M,r}] = \log[A_M] + (r-1)\log[ q_M ] $$

where we rename the quantities $\log q \mapsto \alpha $ and $\log A \mapsto \beta$, and obtain the expression

$$ \log(\mu_{m,r}) = \alpha_m + \beta_m (r-1) $$

which looks exactly the same as the Poisson regression written in the GLM formalism!

Frequentist fit #

Before proceeding with a (arguably) more complete version of a Bayesian model, we will fit the Poisson GLM. The result for one of the decay modes can be illustrated as below

After fitting a Poisson model, there are multiple ways of assessing the quality of the fitting method/fit. One of the most common things is the post-hoc check the Poissonian equidispersion assumption, which is the known property of variance being equal to the mean.

We say that the data is overdispersed (the most common case in real-world data) when the variance of the data - how the observed data deviates from the prediction/mean is not of the order of the predicted mean (Var=Mean).

The most common way to do that is to use Pearson’s residuals ¹. The Pearson’s residual is defined as

$$ r_i^P = \frac{y_i - \hat{\mu}_i}{\sqrt{ \hat{\mu}_i }} $$

The variance of the residual can be obtained by taking the variance of $y_i$, which is $\mu_i$ (Poissonian property), and the denominator gives a $1/\mu$ scaling, giving that $\text{Var}(r_i^P) \approx 1$. Same for the expectation of $r_i^P$. We can then sum the squares of $(r_i^P)^2$ and sum over $i$’s, and it turns out (known properties of the Pearson’s residuals) that the sum denoted as Pearson’s statistic $X_P^2$ has a mean of $n-p$ under the equidispersion assumption. Thus, if the ratio

$$ X_P^2 /(n-p) $$

is strongly higher/lower than $1$, the data is over/under dispersed. In our case, the data is overdispersed, which is expected for a read industrial data. The most straightforward solution for that would be to use the quasi-Poissonian family, which will “artificially inflate” the variance or to use the negative Binomial distribution to fit the data.

Bayesian approach #

After a simple statsmodels model definition and fit, we will use the Bayesian approach to fit the same model. (as it is known to be much more intuitive and slightly scientifically superior).

Since the models implies a different Poisson regression for each failure mode $m$, this technically becomes a multilevel model. The model for every dataframe record $i$ is given by

$$ C_i \sim \text{Pois}(\mu_i) \\ \log(\mu_{i}) = \alpha_{m[i]} + \beta_{m[i]}(r-1) $$

where $C_i$ - the count we’re trying to model. We sample the model and we obtain the corresponding coefficients. The decay rate is obtained through $q_m = e^{\beta_m}$.

When inspecting the parameter distributions, it is natural to assess how likely it is that one of the failure mode is more/less persistent than the other. Within the Bayesian framework one is naturally led to compute it directly from distributions as $q_1 - q_2$ and compare it to $0$.

As mentionned before, after fitting a Poisson-like model, it is handy to appraise the overdispersion, which is again nicely computed using full distributions. We compare the deviation of predicted vs observed and predicted vs generated by the model. The two statistics are $T_\text{obs}$ and $T_\text{ppc}$ (posterior predictive check). A large ratio of obs to ppc signifies that the observed variance is greater than the one predicted by the model.

Next steps #

Based on evidence on overdispersion, we may consider an alternative, for example the negative binomial model. In our case, we will proceed with a slightly more involved dataset, that will be the base to a more informative geometric model.

The geometric model #

We now return to the previous dataset, that has the form

timestamp	id	result
13:32:10	#AA3	pass
13:32:20	#AA4	failmode1
13:32:30	#AA5	pass

Let’s recall the geometric model! The question it models “the number of trials needed to achieve the first success in a series of independent Bernoulli trials”, or simply put, “the number of attempts before a success”. Clearly, this question feels somewhat similar to our problem. Note that we clearly oversimplify our problem - by e.g. assuming that the trials are independent.

Let’s define the parameters and set up the problem. For a unit with id $i$, let’s define $F_i$ as the $\text{number of failed attempts before success or scrap}$. Our end goal would be to esimate the probability $p$, that we assume to be dependent on the mode. So we define $p_{m_i}$ - the probability of pass given the unit $i$ that belongs to mode $m$. There are multiple ways of defining the mode $m$, that the unit $i$ belongs to. Here, we will define as the first failure mode.

For example, for a given unit belonging to mode $m$,

$$ \mathbb{P}(F=0) = p_m \newline \mathbb{P}(F=1) = p_m(1-p_m)\newline \mathbb{P}(F=2) = p_m(1-p_m)^2 \newline ... $$

Given that we have a maximum retest number $K$, the likelihood can be written as

$$ L_i (p_ {m_i}) = p_{m_i}^{s_i} (1 - p_{m_i} )^{f_i} $$

where we define $f_i = \text{ the number of failed attempts}$ and

$$ s_i= \begin{cases} 1 & \text{ passed} \\\\ 0 & \text{ failed more than } K \text{ times} \end{cases} $$

and the total likelihood by mode

$$ L_m(p_m) = \prod_{i: m_i \in m} p_{m_i}^{s_i} (1 - p_{m_i} )^{f_i} $$

from which we define the total number of success belonging to the mode

$$ S_m = \sum_{i: m_i \in m } s_{m_i} $$

and the total number of failed attempts as

$$ F_m = \sum_{i: m_i \in m } f_{m_i} $$

from which the log-likelihood is given by

$$ \text{loglik}_m = S_m \log(p_m) + F_m \log( 1- p_m) $$

Finding the maximum gives the maximum likelihood estimator of $p_m$ as

$$ \hat{p}_m = \frac{S_m}{F_m+S_m} $$

Short note: In order to fit the probabilities of a given failure mode, we have shown that a natural distribution is the Geometrical one. We have mentioned that it follows from sequential Bernoulli trials. Therefore we may even fit the Bernoulli distribution directly without passing through the Geometrical distribution.

We model our Geometrical distribution for bayesian inference as

$$ \ell_i = s_i \log(p_{m_i}) + f_{m_i} \log( 1- p_{m_i}) $$

where we define the prior for $p$ as $p_m \sim \text{logit} (\mathcal{N}(0,1))$. Here, as opposed to a “usual” model specification, we have an explicit log-likelihood expression, which will be defined through the PyMC’s Potential() functionality.

Before fitting, we will slightly aggregate the data - we take only unique ID’s of the units and the computed values $s_i$, $f_i$ and the assigned mode $m_i$. Within the model, we estimate the per-mode probability

Comparing the two baseline models #

Our problem of test/re-test has been evaluated/fitted using two models. It is clear that the underlying process of the related random variables are different.

The poisson distribution describes the number of events that will occur at some known rate. And in this case, we took the event to be the retest. Within these assumptions the events are independent and have constant rates (which are the assumptions) we are willing to accept.

The geometric distribution is based on the Bernoulli trials and answers a question of “how many trials until a success”. Although we know that “retesting” and event “trial until success” are related, we are still asking a different question.

We can plot the estimated parameters of the two methods side-by-side.

We see that the estimators are different for some of the failure modes. We see that the Poisson model probabilities are (almost) always inferior to the Geometric ones. The differences are sometimes considerable. The FMODE_D, however, is similar using both models, which was the initial goal of the problem - to spot significantly persistent modes! To attempt to justify the differences we must consider each of them separately.

The success mode, clearly, does not make sense to be estimated with the geometric distribution. Indeed, we recall that the question is “what is the expected retries before success given the unit belongs to a failure mode $m$”. Where $m$ was defined as the first mode the unit occurred with. So if a unit belongs to the SUCCESS mode, we do not expect to have re-tries.

The failure mode D is the mode that was intended to be caught as the persistent mode, and both models highly agree with a high value of $p_{m}$.

Failure modes A (B,C). The Poisson model estimates how the count of a mode changes across retest ranks. It does not necessarily know whether the same units are carrying the same mode across time. A unit can start with FMODE_A, then later fail with FMODE_B or FMODE_C. In that case, the histogram count of FMODE_A may drop quickly, even if units that initially had FMODE_A continue to fail under different labels.

The geometric model, on the other hand, assigns each unit to a single mode, such as the first observed failure mode. It then asks whether units in that group keep failing. Therefore, a mode can have low Poisson persistence but high geometric persistence if the failure label disappears from the histogram while the units themselves continue failing under other modes. This is the case, where the unit is actually defect and is flagged by belonging to some mode.

The Interactive notebook #

Since the notebook’s functionality relies heavily on PyMC and large data, it is problematic to run the marimo notebook in WASM embedded environment. For this reason, the notebook is directly available on my github, where a single notebook is available.