Designing hypothesis test around failure rates

Question

I'm trying to devise a hypothesis test for failure rate data of machines. The gist is that there are some machines in a factory that run all the time. They fail from time to time and are promptly repaired when they do. Now, a "fix" is deployed to some of the machines (treatment group) and we want to see if it improves (decreases) their failure rate. This is similar to another question I asked a while ago: Hypothesis test for machine failure rates. However, I didn't get an answer back then. What has changed now is that I have taken a stab at it myself and ask the community to please take a look and see if its utterly stupid; what blind spots there might be and what might have been done differently. Next, I will describe the hypothesis test I devised:

---

I assume that the time to repair for the machines is negligible compared to the time to failure in general. I also assume that the exponential distribution is a reasonable assumption for the inter-arrival time of failure events (this is a reasonable assumption since we're doing the test on failure rates and constant failure rates imply the exponential distribution).

Let's say we see $n_1$ downtime events in the first group with total run time (machine hours): $t$ and $n_2$ downtime events in the second group with total run time: $s$. The failure rate of the first group becomes: $\frac{n_1}{t}$ and that of the second group becomes $\frac{n_2}{s}$. So, the difference in failure rate between the two groups becomes:

$$d = \frac{n_1}{t}-\frac{n_2}{s}\tag{1}$$

Now, we want to get the p-value. The null hypothesis is that the two groups have the same failure rate and any $d$ (failure rate difference) we see is statistical noise. So, the null hypothesis assumes a unified failure rate across the two groups of:

$$λ_m=(n_1+n_2)/(t+s)$$

Now, given $t$ and $s$ and $\lambda_m$, the number of failures we expect in those intervals; ($N_1$ and $N_2$ respectively) are Poisson distributed with parameters $λ_m t$ and $λ_m s$ (under the null hypothesis). The difference in two rates we will see will be: $\delta=(N_1/t−N_2/s)$ under the null hypothesis.

So to reject the null, we just get the probability that $\delta > d$ and this becomes the p-value. This can be calculated with simulation or through the double summation (across $N_1$ and $N_2$). Both approaches are implemented here: https://github.com/ryu577/stochproc/blob/master/stochproc/hypothesis/rate.py

Answer 1

Let

• $N_1, N_2, N=N_1+N_2$ denote the (random) failures in durations $t_1, t_2, t=t_1+t_2$
• we observe that $N_1=n_1, N_2=n_2, N=n=n_1+n_2$
• let $\lambda$ be the common (unknown) failure rate under null hypothesis $H_0$.

Here is one way to design a hypothesis test:

$ \begin{align} \mathbb{P}(N_i=n_i) &= \frac{(\lambda t_i)^{n_i} e^{-\lambda t_i}}{n_i!} \\ \mathbb{P}(N=n) &= \frac{(\lambda t)^{n} e^{-\lambda t}}{n!} \\ p_1(n_1,n) \triangleq \mathbb{P}(N_1=n_1 | N=n) &= \frac{\mathbb{P}(N_1=n_1, N=n)}{\mathbb{P}(N=n)} \\ &= \frac{\mathbb{P}(N_1=n_1) ~ \mathbb{P}(N=n | N_1=n_1)}{\mathbb{P}(N=n)} \\ &= \frac{\mathbb{P}(N_1=n_1) ~ \mathbb{P}(N_2=n_2)}{\mathbb{P}(N=n)} \\ &= \frac{t_1^{n_1} t_2^{n_2}}{ t^n {n \choose n_1}}~~~\text{(after simplifying)} \end{align} $

Now, $p_1$ should be high when $n_1 / n \approx t_1 / t$ and it should decrease as $n_1/n \rightarrow 0,1$. Any deviations from it can be used to define p-value. For example, if $t_1=t_2$ we expect $n_1/n$ to be close to $1/2$ under $H_0$. Let $\delta = |n_1/n - 1/2|$ and then p-value $p=\sum_{0 \leq k \leq 1/2 - \delta n} p_1(k, n) + \sum_{1/2 + \delta n \leq k \leq n} p_1(k, n)$. When $t_1 \ne t_2$, the (summation) intervals should preferably not be symmetric. Also, for large $n$ there would be approximations for the tails to make this more tractable.

---

In your question, you define the statistic $d$ that you'd like to use for testing (while in my approach it seems to jump out a bit more naturally) - there's nothing wrong with that. What is an issue though is if the distribution of $d$ depends on $\lambda$ (nuisance parameter?). We must get rid of this dependence. You try to do that by substituting it with its MLE, which is only an approximation to its full distribution.