Cálculo de potencia

Introducción

Se ha utilizado el ejemplo de los efectos de dos dietas diferentes sobre el peso de una población de ratones. Dado que en este ejemplo ilustrativo se tiene acceso a la población, se sabe que, de hecho, existe una diferencia sustancial (alrededor del 10%) entre los pesos medios de las dos poblaciones:

library(dplyr)
dat <- read.csv("mice_pheno.csv") #Previously downloaded 

controlPopulation <- filter(dat,Sex == "F" & Diet == "chow") %>%  
  select(Bodyweight) %>% unlist

hfPopulation <- filter(dat,Sex == "F" & Diet == "hf") %>%  
  select(Bodyweight) %>% unlist

mu_hf <- mean(hfPopulation)
mu_control <- mean(controlPopulation)
print(mu_hf - mu_control)
## [1] 2.375517
print((mu_hf - mu_control)/mu_control * 100) #percent increase
## [1] 9.942157

También se ha visto que, en algunos casos, cuando se toma una muestra y se realiza una prueba t, no siempre se obtiene un valor p menor que 0.05. Por ejemplo, aquí hay un caso en el que se toma una muestra de 5 ratones y no se logra una significancia estadística al nivel de 0.05:

set.seed(1)
N <- 5
hf <- sample(hfPopulation,N)
control <- sample(controlPopulation,N)
t.test(hf,control)$p.value
## [1] 0.5806661

¿Se cometió un error? Al no rechazar la hipótesis nula, ¿se está diciendo que la dieta no tiene ningún efecto? La respuesta a esta pregunta es no. Todo lo que se puede decir es que no se rechaza la hipótesis nula. Pero esto no implica necesariamente que la nula sea verdadera. El problema es que, en este caso particular, no se tiene suficiente potencia, un término que ahora se va a definir. Si está realizando una investigación científica, es muy probable que tenga que hacer un cálculo de potencia en algún momento. En muchos casos, es una obligación ética, ya que puede ayudar a evitar sacrificar ratones innecesariamente o limitar el número de sujetos humanos expuestos a riesgos potenciales en un estudio. Aquí se explica qué significa poder estadístico.

Tipos de error

Siempre que se realiza una prueba estadística, se es consciente de que se puede cometer un error. Esta es la razón por la que los valores p no son 0. Bajo el valor nulo, siempre existe una posibilidad positiva, quizás muy pequeña, pero positiva, de que se rechace el valor nulo cuando es cierto. Si el valor p es 0.05, sucederá 1 de cada 20 veces. Este error es llamado error de tipo I por los estadísticos.

Un error de tipo I se define como rechazar el nulo cuando no deberíamos, también se conoce como falso positivo. Entonces, ¿por qué se usam 0.05? ¿No se debería usar 0.000001 para estar realmente seguros? La razón por la que no se usan cortes infinitesimales para evitar errores de tipo I a toda costa es que hay otro error que se puede cometer: no rechazar el nulo cuando se debería. Esto se denomina error de tipo II o falso negativo. El análisis del código R anterior muestra un ejemplo de un falso negativo: no se rechazó la hipótesis nula (en el nivel 0.05) y, debido a que se conoce y se puede observar las verdaderas medias de la población, se sabe que de hecho hay una diferencia. Si se hubiéra utilizado un valor de p de corte de 0,25, no se habría cometido este error. Sin embargo, en general, ¿el investigador se siente cómodo con una tasa de error de tipo I de 1 en 4? Normalmente no lo está.

Los límites de 0,05 y 0,01 son arbitrarios

La mayoría de las revistas y agencias reguladoras insisten frecuentemente en que los resultados sean significativos en los niveles de 0.01 o 0.05. Por supuesto, no hay nada especial en estos números aparte del hecho de que algunos de los primeros artículos sobre valores p utilizaron estos valores como ejemplos. Parte del objetivo de este libro es brindar a los lectores una buena comprensión de los valores p y los intervalos de confianza para que estas elecciones puedan juzgarse de manera informada. Desafortunadamente, en la ciencia, estos límites se aplican de manera algo inconsciente, pero ese tema es parte de un debate complicado.

Cálculo de potencia

La potencia es la probabilidad de rechazar el nulo cuando el nulo es falso. Por supuesto, “cuando el nulo es falso” es una declaración complicada porque puede ser falsa de muchas maneras.

\(\Delta \equiv \mu_Y - \mu_X\) podría ser cualquier cosa y la potencia realmente depende de este parámetro. También depende del error estándar de las estimaciones, que a su vez depende del tamaño de la muestra y de las desviaciones estándar de la población. En la práctica, no se conocen, por lo que generalmente se informa la potencia para varios valores plausibles de \(\Delta\), \(\sigma_X\), $_Y $ y varios tamaños de muestra.

La teoría estadística da fórmulas para calcular la potencia. El paquete pwr realiza estos cálculos por usted. Aquí se ilustrán los conceptos detrás de la potencia codificando simulaciones en R.

Supóngase que el tamaño de la muestra es:

N <- 12

y se rechazará la hipótesis nula en:

alpha <- 0.05

¿Cuál es la potencia con estos datos en particular? Se calculará esta probabilidad volviendo a ejecutar el ejercicio muchas veces y calculando la proporción de veces que se rechaza la hipótesis nula. Específicamente, se ejecutará:

B <- 2000

simulaciones. La simulación es la siguiente: se toma una muestra de tamaño \(N\) de los grupos de control y de tratamiento, se realiza una prueba t comparando estos dos e informamos si el valor p es menor que “alfa” o no. Se escribe una función que hace esto:

reject <- function(N, alpha=0.05){
   hf <- sample(hfPopulation,N) 
   control <- sample(controlPopulation,N)
   pval <- t.test(hf,control)$p.value
   pval < alpha
}

A continuación se muestra un ejemplo de una simulación para un tamaño de muestra de 12. La llamada a “rechazar” responde a la pregunta “¿Se rechaza?”

reject(12)
## [1] TRUE

Ahora se puede usar la función replicate para hacer estoB veces.

rejections <- replicate(B,reject(N))

La potencia es solo la proporción de veces que se rechaa correctamente. Entonces, con \(N=12\) la potencia poder es solo:

mean(rejections)
## [1] 0.2155

Esto explica por qué la prueba t no rechazó cuando se sabía que el valor nulo era falso. Con un tamaño de muestra de solo 12, la potencia es aproximadamente del 23%. Para proteger contra falsos positivos en el nivel de 0.05, se había establecido el umbral en un nivel lo suficientemente alto que resultó en muchos errores de tipo II.

Hay que mejorar la potencia con N. Se usará la función sapply, que aplica una función a cada uno de los elementos de un vector. Se quiere repetir lo anterior para el siguiente tamaño de muestra:

Ns <- seq(5, 50, 5)

Entonces se usa apply así:

power <- sapply(Ns,function(N){
  rejections <- replicate(B, reject(N))
  mean(rejections)
  })

Para cada una de las tres simulaciones, el código anterior devuelve la proporción de veces que se rechaza. No es sorprendente que la potencia aumente con N:

plot(Ns, power, type="b")
Power plotted against sample size.

Power plotted against sample size.

De manera similar, si se cambia el nivel “alfa” en el que se rechaza, el poder cambia. Cuanto menor quiera que sea la probabilidad de error de tipo I, menos potencia se tiene. Otra forma de decir esto es que se intercambian los dos tipos de error. Se puede ver esto escribiendo un código similar, pero manteniendo \(N\) fijo y considerando varios valores de alfa:

N <- 30
alphas <- c(0.1,0.05,0.01,0.001,0.0001)
power <- sapply(alphas,function(alpha){
  rejections <- replicate(B,reject(N,alpha=alpha))
  mean(rejections)
})
plot(alphas, power, xlab="alpha", type="b", log="x")
Power plotted against cut-off.

Power plotted against cut-off.

Tenga en cuenta que el eje x en este último gráfico está en la escala logarítmica.

No hay un potencia “correcta” o un nivel alfa “correcto”, pero es importante que se comprenda lo que significa cada uno.

Para ver esto claramente, se puede crear una gráfica con curvas de potencia versus N. Se muestran varias curvas en la misma gráfica con color que represente el nivel alfa.

Los valores p son arbitrarios bajo la hipótesis alternativa

Otra consecuencia de lo que se ha aprendido sobre la potencia es que los valores p son algo arbitrarios cuando la hipótesis nula no es cierta y, por lo tanto, la hipótesis alternativa es verdadera (la diferencia entre las medias poblacionales no es cero).

Cuando la hipótesis alternativa es cierta, se puede hacer un valor p tan pequeño como se quiera simplemente aumentando el tamaño de la muestra (suponiendo que se tiene una población infinita de la cual tomar muestras). Se puede mostrar esta propiedad de los valores p extrayendo muestras cada vez más grandes de la población y calculando los valores p. Esto funciona porque, en nuestro caso, se sabe que la hipótesis alternativa es cierta, ya que se tiene acceso a las poblaciones y se puede calcular la diferencia en sus medias.

Primero escriba una función que devuelva un valor p para un tamaño de muestra dado $ N $:

calculatePvalue <- function(N) {
   hf <- sample(hfPopulation,N) 
   control <- sample(controlPopulation,N)
   t.test(hf,control)$p.value
}

Aquí se tiene un límite de 200 para la población con una dieta alta en grasas, pero se puede ver el efecto mucho antes de llegar a 200.

Para cada tamaño de muestra, se calacularán algunos valores p. Se puede hacer esto repitiendo cada valor de \(N\) varias veces.

Ns <- seq(10,200,by=10)
Ns_rep <- rep(Ns, each=10)

Nuevamente se usa sapply para ejecutar las simulaciones:

pvalues <- sapply(Ns_rep, calculatePvalue)

Ahora se pueden graficar los 10 valores p que se generaron para cada tamaño de muestra:

plot(Ns_rep, pvalues, log="y", xlab="sample size",
     ylab="p-values")
abline(h=c(.01, .05), col="red", lwd=2)
p-values from random samples at varying sample size. The actual value of the p-values decreases as we increase sample size whenever the alternative hypothesis is true.

p-values from random samples at varying sample size. The actual value of the p-values decreases as we increase sample size whenever the alternative hypothesis is true.

Tenga en cuenta que el eje y es una escala logarítmica y que los valores p muestran una tendencia descendente hasta \(10^{-8}\) a medida que aumenta el tamaño de la muestra. Los límites estándar de 0.01 y 0.05 se indican con líneas rojas horizontales.

Es importante recordar que los valores p no son más interesantes ya que se vuelven muy, muy pequeños. Una vez que se está convencido de rechazar la hipótesis nula en un umbral que se considera razonable, tener un valor p aún más pequeño solo significa que se han muestreado más ratones de los necesarios. Tener un tamaño de muestra más grande ayuda a aumentar la precisión de nuestra estimación de la diferencia \(\Delta\), pero el hecho de que el valor p sea muy pequeño es solo una consecuencia natural de las matemáticas de la prueba. Los valores p se vuelven cada vez más pequeños con el aumento del tamaño de la muestra porque el numerador del estadístico t tiene \(\sqrt{N}\) (para grupos de igual tamaño, y un efecto similar ocurre cuando \(M\neq{N}\)). Por lo tanto, si \(\Delta\) es distinto de cero, el estadístico t aumentará con \(N\).

Por lo tanto, una mejor estadística para informar es el tamaño del efecto con un intervalo de confianza o alguna estadística que le dé al lector una idea del cambio en una escala significativa. Podemos informar el tamaño del efecto como un porcentaje dividiendo la diferencia y el intervalo de confianza por la media de la población de control:

N <- 12
hf <- sample(hfPopulation, N)
control <- sample(controlPopulation, N)
diff <- mean(hf) - mean(control)
diff / mean(control) * 100
## [1] 5.196293
t.test(hf, control)$conf.int / mean(control) * 100
## [1] -7.775307 18.167892
## attr(,"conf.level")
## [1] 0.95

Además, podemos informar una estadística llamada [d de Cohen] (https://en.wikipedia.org/wiki/Effect_size#Cohen.27s_d), que es la diferencia entre los grupos dividida por la desviación estándar combinada de los dos grupos.

sd_pool <- sqrt(((N-1)*var(hf) + (N-1)*var(control))/(2*N - 2))
diff / sd_pool
## [1] 0.3397886

Esto indica cuántas desviaciones estándar de los datos tiene la media del grupo de dieta alta en grasas del grupo de control. Bajo la hipótesis alternativa, a diferencia del estadístico t que se garantiza que aumentará, el tamaño del efecto y la d de Cohen serán más precisos.