Evaluar el poder de una prueba de normalidad (en R)

Quiero evaluar la precisión de las pruebas de normalidad en diferentes tamaños de muestra en R (me doy cuenta de que las pruebas de normalidad pueden ser engañosas ). Por ejemplo, para ver la prueba de Shapiro-Wilk, estoy realizando la siguiente simulación (además de trazar los resultados) y esperaría que a medida que aumenta el tamaño de la muestra, la probabilidad de rechazar el valor nulo disminuye:

n <- 1000
pvalue_mat <- matrix(NA, ncol = 1, nrow = n)

for(i in 10:n){
    x1 <- rnorm(i, mean = 0, sd = 1)
    pvalue_mat[i,] <- shapiro.test(x1)$p.value
}   

plot(pvalue_mat)

Mi pensamiento sería que a medida que crece el tamaño de la muestra debería haber una tasa de rechazo menor, sin embargo, parece bastante uniforme. Creo que estoy malinterpretando esto, todos y cada uno de los pensamientos son bienvenidos.

Está simulando bajo la hipótesis nula (distribución normal), por lo tanto, la tasa de rechazo tenderá al nivel de significancia esperado. Para evaluar el poder, debe simular bajo cualquier distribución no normal. Hay infinitas posibilidades / escenarios (por ejemplo, distribuciones gamma con asimetría creciente, distribución t con df decreciente, etc.) para elegir, dependiendo del alcance de su estudio.

La comprensión del análisis de poder de las pruebas de hipótesis estadísticas puede mejorarse llevando a cabo algunas y observando de cerca los resultados.

Por diseño, una prueba de tamaño está destinada a rechazar la hipótesis nula con una posibilidad de al menos cuando el nulo es verdadero (su tasa de falsos positivos esperada ). $\alpha$α$\alpha$ Cuando tenemos la capacidad (o el lujo) de elegir entre procedimientos alternativos con esta propiedad, preferimos aquellos que (a) realmente se acercan a la tasa nominal de falsos positivos y (b) tienen posibilidades relativamente más altas de rechazar la hipótesis nula cuando es no es verdad.

El segundo criterio requiere que estipulemos de qué manera (s) y por cuánto el nulo no es cierto. En los casos de libros de texto esto es fácil, porque las alternativas tienen un alcance limitado y están claramente especificadas. Con pruebas de distribución como Shapiro-Wilk, las alternativas son mucho más vagas: son "no normales". Al elegir entre las pruebas de distribución, entonces, es probable que el analista tenga que realizar su propio estudio de poder único para evaluar qué tan bien funcionan las pruebas contra hipótesis alternativas más específicas que son preocupantes en el problema en cuestión.

Un ejemplo motivado por la respuesta de Michael Mayer postula que la distribución alternativa puede tener cualidades similares a las de la familia de distribuciones t de Student. Esta familia, parametrizada por un número (así como por ubicación y escala) incluye en el límite de grandes las distribuciones normales.$\nu\ge 1$$\nu$

En cualquier situación, ya sea evaluando el tamaño real de la prueba o su potencia , debemos generar muestras independientes a partir de una distribución específica, ejecutar la prueba en cada muestra y encontrar la velocidad a la que rechaza la hipótesis nula. Sin embargo, hay más información disponible en cualquier resultado de la prueba: su valor P. Al retener el conjunto de valores P producidos durante una simulación de este tipo, luego podemos evaluar la tasa a la que la prueba rechazaría el valor nulo para cualquier valor de nos interese. El corazón del análisis de potencia, entonces, es una subrutina que genera esta distribución del valor P (ya sea por simulación, como se acaba de describir, o, ocasionalmente, con una fórmula teórica). Aquí hay un ejemplo codificado . Sus argumentos incluyen$\alpha$R

• rdist, el nombre de una función para producir una muestra aleatoria de alguna distribución

• n, el tamaño de las muestras a solicitar rdist

• n.iter, el número de tales muestras para obtener

• ..., los parámetros opcionales a los que se debe pasar rdist(como los grados de libertad )$\nu$

Los parámetros restantes controlan la visualización de los resultados; se incluyen principalmente como una conveniencia para generar las cifras en esta respuesta.

sim <- function(rdist, n, n.iter, prefix="",
                breaks=seq(0, 1, length.out=20), alpha=0.05,
                plot=TRUE, ...) {

  # The simulated P-values.
  # NB: The optional arguments "..." are passed to `rdist` to specify
  #     its parameters (if any).
  x <- apply(matrix(rdist(n*n.iter, ...), ncol=n.iter), 2, 
             function(y) shapiro.test(y)$p.value)

  # The histogram of P-values, if requested.
  if (plot) {
    power <- mean(x <= alpha)
    round.n <- 1+ceiling(log(1 + n.iter * power * (1-power), base=10) / 2)
    hist(x[x <= max(breaks)], xlab=paste("P value (n=", n, ")", sep=""), 
         breaks=breaks, 
         main=paste(prefix, "(power=", format(power, digits=round.n), ")", sep=""))
    # Specially color the "significant" part of the histogram
    hist(x[x <= alpha], breaks=breaks, col="#e0404080", add=TRUE)
  }

  # Return the array of P-values for any further processing.
  return(x)
}

Puede ver que el cálculo en realidad solo toma una línea; El resto del código traza el histograma. Para ilustrar, usémoslo para calcular las tasas de falsos positivos esperados. "Tasas" está en plural porque las propiedades de una prueba generalmente varían con el tamaño de la muestra. Dado que es bien sabido que las pruebas de distribución tienen un alto poder frente a alternativas cualitativamente pequeñas cuando los tamaños de muestra son grandes, este estudio se centra en un rango de tamaños de muestra pequeños donde tales pruebas se aplican a menudo en la práctica: generalmente de a Para ahorrar computación vez, solo informo sobre valores de de a$5$$100.$$n$$5$$20.$

n.iter <- 10^5                 # Number of samples to generate
n.spec <- c(5, 10, 20)         # Sample sizes to study
par(mfrow=c(1,length(n.spec))) # Organize subsequent plots into a tableau
system.time(
  invisible(sapply(n.spec, function(n) sim(rnorm, n, n.iter, prefix="DF = Inf ")))
)

Después de especificar los parámetros, este código también es solo una línea. Produce el siguiente resultado:

Esta es la apariencia esperada: los histogramas muestran distribuciones casi uniformes de valores P en todo el rango de a . Con el conjunto tamaño nominal a las simulaciones informan de entre y de los valores de P eran en realidad menor que el umbral: estos son los resultados resaltados en rojo. La cercanía de estas frecuencias al valor nominal atestigua que la prueba de Shapiro-Wilk funciona como se anuncia.$0$$1$$\alpha=0.05,$$.0481$$0.0499$

(Parece haber una tendencia hacia una frecuencia inusualmente alta de valores P cercanos a Esto es poco preocupante, porque en casi todas las aplicaciones los únicos valores P que se observan son o menos).$1$$0.2$

Pasemos ahora a evaluar el poder. El rango completo de valores de para la distribución t de Student se puede estudiar adecuadamente evaluando algunas instancias desde hasta . ¿Cómo sé eso? Realicé algunas ejecuciones preliminares utilizando un número muy pequeño de iteraciones (de a ), lo que no lleva tiempo en absoluto. El código ahora requiere un bucle doble (y en situaciones más complejas a menudo necesitamos bucles triples o cuádruples para acomodar todos los aspectos que necesitamos variar): uno para estudiar cómo varía la potencia con el tamaño de la muestra y otro para estudiar cómo varía con Los grados de libertad. Una vez más, sin embargo, todo se hace en una sola línea de código (la tercera y final):$\nu$$\nu=100$$\nu=1$$100$$1000$

df.spec <- c(64, 16, 4, 2, 1)
par(mfrow=c(length(n.spec), length(df.spec)))
for (n in n.spec) 
  for (df in df.spec)
    tmp <- sim(rt, n, n.iter, prefix=paste("DF =", df, ""), df=df)

Un pequeño estudio de este cuadro proporciona una buena intuición sobre el poder. Me gustaría llamar la atención sobre sus aspectos más destacados y útiles:

• A medida que los grados de libertad se reducen de a la izquierda a a la derecha, cada vez más valores P son pequeños, lo que demuestra que aumenta el poder de discriminar estas distribuciones de una distribución Normal. (El poder se cuantifica en cada título de la trama: es igual a la proporción del área del histograma que es roja).$\nu=64$$\nu=1$

• A medida que el tamaño de la muestra aumenta de en la fila superior a en la parte inferior, la potencia también aumenta.$n=5$$n=20$

• Observe cómo a medida que la distribución alternativa difiere más de la distribución nula y aumenta el tamaño de la muestra, los valores P comienzan a acumularse a la izquierda, pero todavía hay una "cola" de ellos que se extiende hasta . Esto es característico de los estudios de poder. Muestra que la prueba es una apuesta : incluso cuando la hipótesis nula se viola flagrantemente e incluso cuando el tamaño de nuestra muestra es razonablemente grande, nuestra prueba formal puede no producir un resultado significativo.$1$

• Incluso en el caso extremo en la parte inferior derecha, donde se extrae una muestra de de una distribución de t de Student con grado de libertad (una distribución de Cauchy), la potencia no es : hay una probabilidad de que una muestra de variantes iid Cauchy no se considerará significativamente diferente de Normal a un nivel de (es decir, con un confianza).$20$$1$$1$$100 - 86.57 = 13\%$$20$$5\%$$95\%$

• Podríamos evaluar la potencia en cualquier valor de que elijamos coloreando más o menos barras de estos histogramas. Por ejemplo, para evaluar la potencia en , colorea las dos barras de la izquierda en cada histograma y estima su área como una fracción del total.$\alpha$$\alpha=0.10$

  (Esto no funcionará demasiado bien para valores de más pequeños que con esta figura. En la práctica, uno limitaría los histogramas a valores P solo en el rango que se usaría, tal vez de a , y muéstrelos con suficiente detalle para permitir una evaluación visual del poder hasta o incluso . (Para eso está la opción ). El procesamiento posterior de los resultados de la simulación puede proporcionar aún más detalles.$\alpha$$0.05$$0$$20\%$$\alpha=0.01$$\alpha=0.005$breakssim

Es divertido que se pueda deducir tanto de lo que, en efecto, equivale a tres líneas de código: una para simular muestras de iid de una distribución específica, otra para aplicarla a una matriz de distribuciones nulas y la tercera para aplicarla a Un conjunto de distribuciones alternativas. Estos son los tres pasos que intervienen en cualquier análisis de poder: el resto es solo resumir e interpretar los resultados.

(Más que un comentario, quizás no una respuesta completa)

[I] esperaría que a medida que el tamaño de la muestra aumenta, la probabilidad de rechazar el valor nulo disminuye

Dejando de lado las consideraciones de las pruebas sesgadas (que no son infrecuentes en cuanto a la bondad de ajuste, por lo que vale la pena mencionarlas), hay tres situaciones relacionadas con la tasa de rechazo que uno podría considerar:

1) la tasa de rechazo al simular desde el nulo (como parece estar haciendo en su pregunta)

Aquí, la tasa de rechazo debe estar en o cerca del nivel de significancia, por lo tanto, no, si mantiene constante el nivel de significancia, la tasa de rechazo no disminuye a medida que aumenta n , sino que se mantiene en / cerca de .$\alpha$

2) la tasa de rechazo al simular desde alguna alternativa

Aquí la tasa de rechazo debería aumentar a medida que aumenta n .

3) la tasa de rechazo para alguna recopilación de datos reales

Prácticamente, el valor nulo nunca es realmente verdadero, y los datos reales tendrán una mezcla de cantidades de no normalidad (según lo medido por el estadístico de prueba). Si el grado de no normalidad no está relacionado con el tamaño de la muestra, la tasa de rechazo debería aumentar a medida que n aumenta.

De hecho, en ninguna de estas situaciones deberíamos ver que la tasa de rechazo disminuye con el tamaño de la muestra.