test d'hypothèse à un échantillon pour les proportions

Question

Je cherche une fonction R intégrée qui calcule la puissance d'un test d'hypothèse à un échantillon pour les proportions.

La fonction intégrée power.prop.test ne fait que des tests d'hypothèse à DEUX ÉCHANTILLONS pour les proportions.

La question initiale est la suivante : "Combien de fois devez-vous lancer une pièce de monnaie pour déterminer qu'elle est biaisée ?

p.null <- 0.5            # null hypothesis.

On dit qu'une pièce est " biaisée " si la probabilité de lancer face est soit supérieure à 0,51 ou inférieure à 0,49. Dans le cas contraire, on dit qu'elle est "suffisamment bonne".

delta <- 0.01       

Voici une fonction permettant de lancer N fois une pièce de monnaie biaisée et de renvoyer la proportion de têtes :

biased.coin <- function(delta, N) {
  probs <- runif(N, 0, 1) 
  heads <- probs[probs < 0.5+delta]      
  return(length(heads)/N)
}

Nous fixons alpha et bêta aux valeurs standard. Notre objectif est de calculer N.

alpha = 0.05        # 95% confidence interval
beta = 0.8          # Correctly reject the null hypothesis 80% of time.

La première étape consiste à utiliser une simulation.

Une seule expérience consiste à tirer à pile ou face N fois et à rejeter l'hypothèse nulle si le nombre de têtes s'écarte "trop" de la valeur attendue de N/2.

Nous répétons ensuite l'expérience M fois et comptons combien de fois l'hypothèse nulle est (correctement) rejetée.

M <- 1000

simulate.power <- function(delta, N, p.null, M, alpha) {
   print(paste("Calculating power for N =", N))
   reject <- c()
   se <- sqrt(p.null*(1-p.null))/sqrt(N)   
   for (i in (1:M)) {
       heads <- biased.coin(delta, N)       # perform an experiment
       z <- (heads - p.null)/se             # z-score 
       p.value <- pnorm(-abs(z))            # p-value
       reject[i] <- p.value < alpha/2       # Do we rejct the null? 
   }
   return(sum(reject)/M)        # proportion of time null was rejected.
  }

Ensuite, nous traçons un graphique (lentement, environ 5 minutes) :

ns <- seq(1000, 50000, by=1000)
my.pwr <- c()
for (i in (1:length(ns))) {
  my.pwr[i] <- simulate.power(delta, ns[i], p.null, M, alpha)
}
plot(ns, my.pwr)

D'après le graphique, il semble que le N nécessaire pour une puissance de bêta = 0,8 soit d'environ 20 000.

La simulation est très lente et il serait bon d'avoir une fonction intégrée.

Un peu de bricolage m'a donné ça :

magic <- function(p.null, delta, alpha, N) {
  magic <-power.prop.test(p1=p.null, 
                          p2=p.null+delta,
                          sig.level=alpha,

                          ###################################
                          n=2*N,             # mysterious 2
                          ###################################

                          alternative="two.sided",
                          strict=FALSE)
  return(magic[["power"]])                 
}

Comparons-le à nos données simulées.

pwr.magic <- c()
for (i in (1:length(ns))) {
   pwr.magic[i] <- magic(p.null, delta, alpha, ns[i]) 
}
points(ns, pwr.magic, pch=20)

L'ajustement est bon, mais je ne vois pas pourquoi je devrais multiplier N par deux, afin d'obtenir une puissance d'un échantillon à partir d'un test de proportion de deux échantillons.

Ce serait bien s'il y avait une fonction intégrée qui vous permette de faire un échantillon directement.

Gracias.

Answer 1

Vous pouvez essayer

library(pwr)
h <- ES.h(0.51, 0.5) # Compute effect size h for two proportions
pwr.p.test(h = h, n = NULL, sig.level = 0.05, power = 0.8, alternative = "two.sided")
# proportion power calculation for binomial distribution (arcsine transformation) 

#           h = 0.02000133
#           n = 19619.53
#   sig.level = 0.05
#       power = 0.8
# alternative = two.sided

En passant, une façon d'accélérer votre simulation de manière significative serait d'utiliser rbinom au lieu de runif :

biased.coin2 <- function(delta, N) {
  rbinom(1, N, 0.5 + delta) / N
}