Designing hypothesis test around failure rates

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

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.