Cálculo de potência para teste de razão de verossimilhança

8

Eu tenho duas variáveis ​​aleatórias independentes de poisson, e , com e . Desejo testar versus a alternativa .X1X2X1Pois(λ1)X2Pois(λ2)H0:λ1=λ2H1:λ1λ2

Eu já deduzi estimativas de máxima verossimilhança sob hipótese nula e alternativa (modelo) e com base naquelas que calculei a estatística do teste de razão de verossimilhança (LRT) (códigos R fornecidos abaixo).

Agora, estou interessado em calcular a potência do teste com base em:

  1. Alfa corrigida (erro do tipo 1) = 0,05.
  2. Usando diferentes tamanhos de amostra (n), diga n = 5, 10, 20, 50, 100.
  3. Combinação diferente de e , que alterará as estatísticas do LRT (calculadas conforme abaixo).λ1λ2LRTstat

Aqui está o meu código 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)

A partir daqui, como eu poderia proceder com o cálculo da potência (de preferência em R)?

Adão
fonte

Respostas:

9

Você pode fazer isso usando simulação.

Escreva uma função que faça seu teste e aceite os lambdas e o (s) tamanho (s) da amostra como argumento (você teve um bom começo acima).

Agora, para um determinado conjunto de lambdas e tamanhos de amostra, execute a função várias vezes (a função de replicação em R é ótima para isso). Então, o poder é apenas a proporção de vezes que você rejeita a hipótese nula; você pode usar a função média para calcular a proporção e propor o teste para fornecer um intervalo de confiança no poder.

Aqui está um exemplo de código:

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)

Meus resultados (você diferirá com uma semente diferente, mas deve ser semelhante) mostraram uma taxa de erro do tipo I (alfa) de 0,0496 (IC 95% 0,0455-0,0541) que é próxima de 0,05, mais precisão pode ser obtida aumentando-se 10000 no comando replicar. Os poderes que calculei foram: 9,86% e 28,6%. Os histogramas não são estritamente necessários, mas eu gosto de ver os padrões.

Greg Snow
fonte
Eu criei uma função (LRT.POIS) com parâmetros nSim, Lambda1, Lambda2, mas daqui em diante estou meio que perdida.
Adam
2
Eu adicionei algum código de exemplo para mostrar o processo básico.
Greg Snow