Simulation of maximum likelihood ratio test to test two poisson random variables

Question

I have two random poisson variables $x_1$ and $x_2$ with value 10 and 25 respectively. I am interested to use likelihood ratio test to test the null hypothesis: $\lambda_1=\lambda_2$, versus alernate hypthesis $\lambda_1$ not equal to $\lambda_2$.

I want to use simulation to calculate power and alpha values. I would want to do it in R so any reference to R codes will be appreciated. Thanks in advance

Answer 1

This is a particularly ill-formed question.

If by "alpha" you mean Type I error, you need to go back to Square One and get definitions straight. Type I error is not something inherent in the data, or even in the hypothesis; it's a subjectively and externally applied measure of risk. And without the Type I error, you have no reference point from which to calculate Type II error, the complement of power.

Worse yet, it's not clear--after Adam's question--whether you have JUST TWO OBSERVATIONS (10 and 25), or two distributions with means of 10 and 25, and you're looking for a suitable sample size for a balanced test comparing the means. In the first case, all you can do is a likelihood ratio test that gives an approximate p-value; there's no more information to be had in two observations. In the second case, simulation can give some useful results, but you still need a value for the Type I error to get started.

Answer 2

For the simulation let's first choose sample sizes N1 and N2 for the two Poisson samples:

require(lmtest)
N1 = 20; N2 = 15

Generate a random sample and run a likelihood ratio test:

# CODE BLOCK "A"
x = c(rpois(N1, 10), rpois(N2, 25))
group = factor(c(rep('a', N1), rep('b', N2)))
m1 = glm(x ~ 1, family=poisson)
m2 = glm(x ~ group, family=poisson)
(t = lrtest(m1, m2))

Result:

Likelihood ratio test

Model 1: x ~ 1
Model 2: x ~ group
  #Df  LogLik Df  Chisq Pr(>Chisq)    
1   1 -158.26                         
2   2  -93.39  1 129.75  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Now let's run many simulations to see the power for this particular N1 and N2:

s = 1000 # 10*1000 simulations
sigs = NULL
for (i in 1:10) {
  sig = 0
  for (j in i:s) {
    CODE BLOCK "A" COMES HERE
    if (t$Pr[2] <= 0.05) sig = sig + 1
  }
  sigs = c(sigs, sig / s)
}
c(quantile(sigs, c(.025, .5, .975)), mean=mean(sigs), sd=sd(sigs))

Result:

       2.5%         50%       97.5%        mean          sd 
0.991225000 0.995500000 0.999775000 0.995500000 0.003027650

Thus the power for N1 = 20; N2 = 15 is 99%.

You can calculate power for various N1 and N2 values.