3
$\begingroup$

I have an experiments that counts how many events happened during a certain period of time. The experiments is repeated many times in different "runs" (I have 65 runs).

Here for example the list of observed events for every run:

obs = [48, 88, 9, 38, 84, 84, 8, 14, 47, 6, 202, 156, 69, 460, 75, 145, 243, 390, 3, 541, 613, 57, 353, 231, 5, 360, 484, 705, 745, 61, 516, 437, 71, 96, 183, 1121, 405, 66, 1000, 995, 271, 93, 368, 10, 182, 509, 1044, 0, 208, 732, 23, 817, 415, 238, 111, 1765, 188, 1129, 283, 89, 1550, 1974, 471, 281, 825]

and here the duration of every run (the number of events is proportional to the time of observation):

times = [6.25, 11.91, 0.66, 4.22, 9.00, 10.87, 0.68, 1.60, 6.65, 0.74, 21.82, 19.96, 9.96, 55.58, 7.32, 17.10, 32.03, 52.81, 0.35, 68.68, 75.22, 6.77, 48.46, 30.07, 1.02, 45.11, 66.18, 89.72, 95.19, 8.43, 69.69, 59.28, 7.45, 12.42, 24.83, 151.14, 45.64, 7.16, 132.74, 126.59, 35.19, 12.45, 46.30, 0.85, 24.19, 64.07, 136.96, 10.55, 26.52, 93.41, 2.97, 108.62, 57.15, 29.44, 13.52, 235.91, 24.08, 146.92, 39.54, 10.21, 202.20, 249.35, 57.50, 35.37, 104.50]

There is no constraint on the total number of events. There is no correlation between runs.

I want to test the null hypothesis "the rate is the same for every run (two side test: the incompatibility could be for excess or deficit)".

From the data I now:

mean_obs_rate = sum(obs) / sum(times) = 0.0352

So I can compute the expected number of events under the null hypothesis, using the observed mean rate exp = mean_obs_rate * times:

exp = [48.31, 92.00, 5.12, 32.61, 69.55, 83.99, 5.22, 12.33, 51.35, 5.74, 168.58, 154.16, 76.92, 429.37, 56.57, 132.13, 247.46, 407.94, 2.69, 530.54, 581.06, 52.33, 374.35, 232.29, 7.90, 348.46, 511.26, 693.06, 735.33, 65.08, 538.34, 457.92, 57.55, 95.91, 191.79, 1167.55, 352.56, 55.28, 1025.39, 977.87, 271.81, 96.17, 357.68, 6.60, 186.90, 494.95, 1058.01, 81.54, 204.90, 721.59, 22.92, 839.10, 441.52, 227.41, 104.41, 1822.38, 186.02, 1134.93, 305.44, 78.85, 1562.01, 1926.26, 444.23, 273.23, 807.27]

Then I have tried with the likelihood ratio using as a model a product of Poissonian distribution:

$$ -2 \log \frac{\sup_\lambda \text{Pois}(n_1|\lambda t_1)\cdots\text{Pois}(n_{65}|\lambda t_{65})}{\sup_{\{\lambda_i\}} \text{Pois}(n_1|\lambda_1 t_1)\cdots\text{Pois}(n_{65}|\lambda_{65}t_{65})} $$

The two hypotheses are nested (you can write $\lambda_i = \lambda + \delta\lambda_i$ and the numerator is the special case with $\delta\lambda_i=0$)

The denominator is maximized with $\lambda_i t_i = n_i$, while the numerator is maximized with $\lambda t_i = t_i\times$ mean_obs_rate.

After some simplifications I have found that this is equal to:

$$ 2 \sum_i n_i \log\frac{n_i}{e_i} $$ with $e_i = \lambda t_i$. This is exactly the G-test. That should be asymptotically distributed as a $\chi^2$ distribution with 65 - 1 dof (even if generating pseudo experiments under the null-hypothesis I found that there is a small bias). Small number of the G-test should be good compatibility, while large should be bad compatibility with the null-hypothesis. Is this correct, taking into account the fact that I want a two-sided test?

I was trying to see what is the contribution to the G-test of every run ($n_i \log\frac{n_i}{e_i}$), and I looked to one where the number of observed event is 0 and expected is 81.54. It should give a large contribution to the G-test quantity, since its contribution should be proportional to the probability to observe 0 events when the expected is 81.45, but actually its contribution is 0 (just taking the limit for $n_i\to 0$ in the $G$ definition). That is completely the opposite compared to a simpler $\chi^2$-test $(n_i - e_i)^2 / e_i = (0 - 81.54)^2 / 81.54 = 81.54$.

So the question: is this ok? How should I implement runs where I observe 0 events? Why runs with 0 observed events count nothing in the final value of the G-test? Should they contribute in the number of degree of freedom of the asymptotic $\chi^2$ distribution?

$\endgroup$
1
$\begingroup$

The likelihood ratio

$$ -2\log\Lambda = -2 \log \frac{\sup_\lambda \text{Pois}(n_1|\lambda t_1)\cdots\text{Pois}(n_{65}|\lambda t_{65})}{\sup_{\{\lambda_i\}} \text{Pois}(n_1|\lambda_1 t_1)\cdots\text{Pois}(n_{65}|\lambda_{65}t_{65})} $$

can be simplified to:

$$ -2\log\Lambda = 2\sum (e_i - n_i + n_i\log(n_i/e_i)) $$

then using $\sum e_i=\sum n_i$, that hold in this case, the expression can be further simplified to:

$$ -2\log\Lambda = 2\sum n_i\log(n_i/e_i) $$

The point is that the last simplification holds when doing the sum, not when considering each single term of the summation. So, if I want to see the contribution of each term to the G-test I have to consider the quantity $2 (e_i - n_i + n_i\log(n_i/e_i))$, which has a minimum for $e_i=n_i$.

In the figure the dotted line is the expression above, while the dotted lines are the chi-square $(e_i - n_i)^2 / e_i$ for different values of $n$ ($O$ in the figure) as a function of $e_i$ ($E$ in the figure).

gtest_vs_chisquare

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.