Calcul de puissance pour le test du rapport de vraisemblance

8

J'ai deux variables aléatoires de poisson indépendantes, et , avec et . Je veux tester contre l'alternative .X1X2X1Pois(λ1)X2Pois(λ2)H0:λ1=λ2H1:λ1λ2

J'ai déjà dérivé des estimations du maximum de vraisemblance sous hypothèse nulle et alternative (modèle), et sur la base de celles-ci, j'ai calculé la statistique du test de rapport de vraisemblance (TLR) (codes R donnés ci-dessous).

Maintenant, je suis intéressé à calculer la puissance du test en fonction de:

  1. Alpha fixe (erreur de type 1) = 0,05.
  2. En utilisant différentes tailles d'échantillon (n), disons n = 5, 10, 20, 50, 100.
  3. Combinaison différente de et , ce qui changera les statistiques LRT (calculées comme ci-dessous).λ1λ2LRTstat

Voici mon code 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)

À partir de là, comment pourrais-je procéder au calcul de la puissance (de préférence en R)?

Adam
la source

Réponses:

9

Vous pouvez le faire en utilisant la simulation.

Écrivez une fonction qui fait votre test et accepte les lambdas et la taille de l'échantillon comme arguments (vous avez un bon début ci-dessus).

Maintenant, pour un ensemble donné de lambdas et de taille d'échantillon, exécutez la fonction plusieurs fois (la fonction de réplication dans R est idéale pour cela). Ensuite, la puissance n'est que la proportion de fois où vous rejetez l'hypothèse nulle, vous pouvez utiliser la fonction moyenne pour calculer la proportion et prop.test pour donner un intervalle de confiance sur la puissance.

Voici un exemple de 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)

Mes résultats (votre différence avec une graine différente, mais devrait être similaire) ont montré un taux d'erreur de type I (alpha) de 0,0496 (IC à 95% 0,0455-0,0541) qui est proche de 0,05, plus de précision peut être obtenue en augmentant le 10000 dans la commande de réplication. Les puissances que j'ai calculées étaient: 9,86% et 28,6%. Les histogrammes ne sont pas strictement nécessaires, mais j'aime voir les motifs.

Greg Snow
la source
J'ai créé une fonction (LRT.POIS) avec des paramètres, nSim, Lambda1, Lambda2, mais à partir de là, je suis un peu perdu.
Adam
2
J'ai ajouté un exemple de code pour montrer le processus de base.
Greg Snow