Power calculation for likelihood ratio test

Question

I have two independent poisson random variables, $X_1$ and $X_2$, with $X_1 \sim \text{Pois}(\lambda_1)$ and $X_2 \sim \text{Pois}(\lambda_2)$. I want to test $H_0:\, \lambda_1 = \lambda_2$ versus the alternative $H_1:\, \lambda_1 \neq \lambda_2$.

I already derived maximum likelihood estimates under null and alternate hypothesis (model), and based on those I calculated likelihood ratio test (LRT) statistic (R codes given below).

Now I am interested to calculate power of the test based on:

1. Fixed alpha (type 1 error) = 0.05.
2. Using different sample sizes (n), say n = 5, 10, 20, 50, 100.
3. Different combination of $\lambda_1$ and $\lambda_2$, which will change LRT statistics (computed as LRTstat below).

Here is my R code:

X1 = rpois(λ1); X2 = rpois(λ2)
Xbar = (X1+X2)/2
LLRNum = dpois(X1, X1) * dpois(X2, X2)
LLRDenom = dpois(X1, Xbar) * dpois(X2, Xbar)
LRTstat = 2*log(LLRNum/LLRDenom)

From here, how could I proceed with power calculation (preferably in R)?

Answer 1

You can do this using simulation.

Write a function that does your test and accepts the lambdas and sample size(s) as arguments (you have a good start above).

Now for a given set of lambdas and sample size(s) run the function a bunch of times (the replicate function in R is great for that). Then the power is just the proportion of times that you reject the null hypothesis, you can use the mean function to compute the proportion and prop.test to give a confidence interval on the power.

Here is some example code:

tmpfunc1 <- function(l1, l2=l1, n1=10, n2=n1) {
    x1 <- rpois(n1, l1)
    x2 <- rpois(n2, l2)
    m1 <- mean(x1)
    m2 <- mean(x2)
    m <- mean( c(x1,x2) )

    ll <- sum( dpois(x1, m1, log=TRUE) ) + sum( dpois(x2, m2, log=TRUE) ) - 
            sum( dpois(x1, m, log=TRUE) ) - sum( dpois(x2, m, log=TRUE) )
    pchisq(2*ll, 1, lower=FALSE)
}

# verify under null n=10

out1 <- replicate(10000, tmpfunc1(3))
mean(out1 <= 0.05)
hist(out1)
prop.test( sum(out1<=0.05), 10000 )$conf.int

# power for l1=3, l2=3.5, n1=n2=10
out2 <- replicate(10000, tmpfunc1(3,3.5))
mean(out2 <= 0.05)
hist(out2)

# power for l1=3, l2=3.5, n1=n2=50
out3 <- replicate(10000, tmpfunc1(3,3.5,n1=50))
mean(out3 <= 0.05)
hist(out3)

My results (your will differ with a different seed, but should be similar) showed a type I error rate (alpha) of 0.0496 (95% CI 0.0455-0.0541) which is close to 0.05, more precision can be obtained by increasing the 10000 in the replicate command. The powers I computed were: 9.86% and 28.6%. The histograms are not strictly necessary, but I like seeing the patterns.