analyse de puissance pour le test binomial

Sep 02 2020

J'essaie de tester si le sex-ratio de certains individus échantillonnés diffère significativement du sex-ratio attendu de 1. J'ai n = 64, dont femme = 34 et homme = 30.

J'ai effectué un test binomial:

succ <- c(34,30) 
binom.test(succ,0.5)

data:  succ
number of successes = 34, number of trials = 64, p-value = 0.708
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.4023098 0.6572035
sample estimates:
probability of success 
               0.53125 

Je voudrais calculer la puissance statistique de ce test, et je sais que puissance = 1-β , où β est l'erreur de type II.

Je suis confus en lisant cette explication . Je ne comprends pas comment adapter cette formule (pour différents choix de n) à mon cas:

enn = 1:2000
critical = qbinom(.025, enn, .5)
beta = pbinom(enn-critical,enn,.55) - pbinom(critical-1,enn,.55)

Ce que j'ai fait était

1-(pbinom(34,64,0.5)- pbinom(30, 64, .5))
[1] 0.7410237

mais je ne suis pas sûr qu'il soit correct d'utiliser 0,5 comme probabilité. De plus, j'ai essayé une méthode différente, et j'obtiens un résultat complètement différent

pwr.p.test(ES.h(.53125,.5),n=64, power=NULL, alternative = "two.sided")

 proportion power calculation for binomial distribution (arcsine transformation) 

              h = 0.06254076
              n = 64
      sig.level = 0.05
          power = 0.07913605
    alternative = two.sided

L'un de ces deux tests est-il correct et pourquoi?

Merci pour ton aide!

Réponses

5 BruceET Sep 02 2020 at 04:48

Afin de trouver le «pouvoir», vous devez avoir une alternative spécifique en tête. Supposons que votre hypothèse nulle soit$H_0: p = 0.5$ vs. $H_a: p > 0.5,$$p = P(\mathrm{Female}).$ Supposons également que vous ayez $n = 64$ et vous voulez la puissance d'un test au niveau $\alpha = 0.05$ contre l'alternative spécifique $p = 0.6.$

Pour un test binomial exact, vous devez trouver la valeur critique $c$ tel que $P(X \ge c\,|\,n=64, p=.5)$ est maximisé, mais toujours en dessous $0.05.$Dans R, où dbinom, pbinomet qbinomdésignent respectivement le binomial PDF, CDF et la fonction quantile (CDF inverse), nous voyons que la valeur critique est$c = 40.$ Notez qu'en raison de la discrétion des distributions binomiales, le soi-disant `` niveau de 5% '' rejette en fait avec probabilité $P(\mathrm{Rej}\, H_0 | H_0\, \mathrm{True}) \approx 3\%.$

qbinom(.95, 64, .5)
[1] 39
sum(dbinom(39:64, 64, .5))
[1] 0.05171094
sum(dbinom(40:64, 64, .5))
[1] 0.02997059
1 - pbinom(39, 64, .5)
[1] 0.02997059

Ensuite, la puissance de ce test par rapport à la valeur alternative $p = 0.6$ est donné par $P(X \ge 40\,|\,n=64, p=0.6) = 0.3927.$

1 - pbinom(39, 64, .6)
[1] 0.392654

Nous pouvons faire une `` courbe de puissance '' pour ce test en regardant une séquence de valeurs alternatives p.aentre$0.5$ et $.75.$ Le premier bloc de code R ci-dessous forme la ligne noire continue dans le graphique ci-dessous.

p.a = seq(.50, .75, by=.01)
p.rej = 1 - pbinom(39, 64, p.a)
plot(p.a, p.rej, type="l", main="Power Curve")
 abline(h=c(.03,1), col="green2")

Si nous regardons un niveau $\alpha = 0.05$ test de $H_0: p = 0.5$ vs $H_a: p > 0.5$ avec $n = 256$ sujets, alors la valeur critique est $c = 141,$ la probabilité de rejet lorsque $H_0$ est vrai est $0.046,$ et le pouvoir contre diverses valeurs alternatives de $p$ est supérieure, comme le montre la ligne bleue en pointillé sur la figure.

c.256 = qbinom(.95, 256, .5); c.256
[1] 141
1 - pbinom(c.256, 256, .5)
[1] 0.04565604
p.rej.256 = 1 - pbinom(c.256, 256, p.a)
lines(p.a, p.rej.256, col="blue", lty="dotted")

Remarques: Parce que$n = 64$est suffisamment grand pour utiliser des approximations normales, vous voudrez peut-être essayer d'utiliser des approximations normales. Un inconvénient est que cela ignore le problème de la discrétion, il peut donc sembler que votre test rejette exactement 5% du temps lorsque$H_0$est vrai. En outre, vous devez utiliser une correction de continuité pour de meilleurs résultats.

Un calcul pertinent pour le niveau de signification dans R est:

1 - pnorm(39.5, 32, 4)
[1] 0.03039636

(Approximative) la puissance est $0.3895:$

mu.a = 64*.6;  sg.a = sqrt(64*.6*.4)
mu.a; sg.a
[1] 38.4
[1] 3.919184

1 - pnorm(39.5, mu.a, sg.a)     # Using NORM(mu.a, sg.a)
[1] 0.3894815
1 - pnorm((39.5 - mu.a)/sg.a)   # Standardizing and using NORM(0,1).
[1] 0.3894815