我有两个独立的泊松随机变量，和， 和和. 我想测试与替代方案.
我已经在零假设和替代假设（模型）下得出了最大似然估计，并基于我计算的似然比检验（LRT）统计量（下面给出的 R 代码）。
现在我有兴趣根据以下内容计算测试的功效：
LRTstat如下）。这是我的 R 代码：
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)
从这里开始，我该如何进行功率计算（最好在 R 中）？
您可以使用模拟来做到这一点。
编写一个函数来进行测试并接受 lambdas 和样本大小作为参数（你在上面有一个好的开始）。
现在，对于给定的一组 lambdas 和样本大小，多次运行该函数（R 中的复制函数对此非常有用）。那么功效就是您拒绝原假设的次数的比例，您可以使用 mean 函数来计算比例，并使用 prop.test 给出功效的置信区间。
这是一些示例代码：
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)
我的结果（您会因不同的种子而有所不同，但应该相似）显示 I 型错误率 (alpha) 为 0.0496（95% CI 0.0455-0.0541），接近 0.05，通过增加 10000 可以获得更高的精度在复制命令中。我计算的功效是：9.86% 和 28.6%。直方图不是绝对必要的，但我喜欢看到这些模式。