Bioestadística práctica.

SECUAH 2017

Marcos Marvá y Fernando San Segundo. Depto. de Física y Matemáticas, Universidad de Alcalá.

Replicación.

La crisis de replicación.

  • Podemos empezar por alguno de los (ahora muy numerosos) artículos que han levantado la liebre, como este de Nature (2016-05).

  • Al principio la discusión era sobre todo académica. Pero ahora anda suelta y empezamos a verla en medios generalistas, como la BBC (2017-02) o EL PAÍS (2017-01).

Aclaración sobre terminología:

  • Al principio se hablaba sobre todo de reproducibilidad. Pero ahora vamos afinando el vocabulario y preferimos distinguir entre replicabilidad y reproducibilidad. Aquí replicar se refiere a la idea intuitiva de repetir el experimento y obtener los mismos resultados.

  • La reproducibilidad, en el sentido que vamos a darle aquí, tiene más que ver con la capacidad de repetir el análisis de los datos generado por el experimento. Es una cuestión básicamente computacional. Vemos la reproducibilidad como una condición necesaria para la replicación. Pero a veces, en experimentos complejos, caros, etc., es casi la única alternativa disponible.

Análisis reproducibles.

Hay una tendencia creciente a favor de vetar las publicaciones científicas que no contengan estos ingredientes (cuando son la base de las conclusiones):

  1. Datos accesibles. Ver por ejemplo la política de Nature
  2. Código del análisis de datos. De nuevo, Nature
  3. Descripción adecuada de los métodos estadísticos: Por ejemplo, en Science


En resumen, no se consideran aceptables los trabajos que con un análisis estadístico no reproducible (por ejemplo, que usen cualquier cosa que se haga a base de clicks de ratón en una serie de menús).

Herramientas para el análisis reproducible.

A estas alturas el proyecto R necesita poca presentación. Pero por si acaso ahí van unos datos.

  • R se usa mucho en general. Aquí y aquí hay comparativas con otros lenguajes. R es el primer lenguaje no generalista que aparece (por delante de lenguajes tan establecidos en el mundo de la técnica como Matlab). Y para ampliar aquí hay un análisis relativamente reciente (2015) de la popularidad de muchos programas estadísticos, que incluye esta gráfica sobre el uso en trabajos académicos (reproduccción parcial; en rojo he añadido el año de la 1ª versión):
  • R no para de crecer. El principal repositorio de librerías de R, llamado CRAN, ha experimentado este tipo de crecimiento en los últimos años:

Bioconductor, el principal repositorio para bioinformática, contiene cerca de 1300 librerías.

¿Quién es quién?

  • R es el motor del sistema. Es un lenguaje de programación orientado al análisis estadístico de datos.

  • R consta de un sistema base, al que podemos añadir una gran cantidad de librerías adicionales para las tareas más variopintas.

  • Rmarkdown es una herramienta de Programación Literaria. Por lo tanto combina en un mismo fichero el código del programa en R con la documentación de ese código. Al procesar (tejer) un fichero Rmarkdown se obtiene un informe que contiene los resultados del análisis. Esta presentación está escrita integramente usando Rmarkdown.

  • RStudio es un programa que hace de intermediario entre nosotros y R, facilitando el trabajo. Rstudio ofrece muchos recursos para que la experiencia de trabajar con R sea muy productiva.

Hay otras opciones para el análisis reproducible, siendo quizá el proyecto Python / Jupyter la más destacada. Pero creemos que en este momento R sigue siendo la opción más atractiva. Además es código abierto. Y del precio ni hablamos...

Un documento Rmarkdown básico.

Para empezar a descubrir el funcionamiento de esas herramientas hemos preparado un documento muy sencillo pero que ilustra las ideas fundamentales de este método de trabajo.

Ejemplo-Rmarkdown.Rmd

Vamos a descargar ese fichero y a jugar un poco con él para entender cómo funciona. Usaremos RStudio para obtener versiones de este informe básico en HTML, MsWord y PDF. Esta última opción es muy interesante con vistas a la publicación, pero requiere la instalación de un programa adicional, un compilador del lenguaje LaTeX. Para Windows se recomienda instalar MikTeX.

Ejemplo extendido. Un análisis exploratorio.

En nuestro siguiente ejemplo, más extenso, trataremos de profundizar más en el uso de R, ilustrando cómo se usa en la exploración inicial de un conjunto de datos. Vamos a basarnos en un ejemplo que aparece en el Capítulo 1 del libro "Linear Models with R" de J.Faraway.

El fichero que hemos preparado para este ejemplo es

EjemploAnalisisExploratorio.Rmd

Descárgalo y guárdalo en tu directorio de trabajo. Vamos a ir ejecutando juntos ese fichero para discutir la exploración de los datos que describe, que puede servir de modelo para otras situaciones.

Volvamos sobre los datos.

El fichero exploratorio que hemos usado corresponde a la fase inicial del análisis de un conjunto de datos. A menudo ese tipo de análisis requiere la colaboración de la persona que realiza el estudio (el científico) con un estadístico o analista de datos. Y, siendo políticamente correctos, no siempre es fácil entenderse.

Recomendamos encarecidamente la lectura de la guía de Jeff Leek sobre cómo trabajar con un estadístico / analista de datos.

En particular, en esa guía se menciona el concepto fundamental de datos limpios, o tidy data sets. No podemos entrar en detalle, pero al menos queremos indicar tres principios básicos sobre las tablas de datos:

  1. Cada variable forma una columna (y sólo una) de la tabla.
  2. Cada observación forma una fila de la tabla.
  3. Cada tipo de unidad observacional aparece en su propia tabla.

Para aprender más sobre R y Rmarkdown.

Las siguientes listas muestran algunos recursos básicos. Pero si queréis información más detallada sobre algún campo en concreto, poneos en contacto con nosotros y ampliaremos esta información.

Bibliografía

  • Claro, para refrescar la Estadística básica (la de primero) y para aprender R a la vez, nuestro favorito no puede ser otro que PostData 1.0.

  • Un libro razonablemente bueno para iniciarse en R es The Art of R Programming de Norman Matloff. Tiene capítulos irregulares pero el promedio es bueno. Por si lo veis por ahí, creo que es mejor que el de Dalgaard, que se cita a menudo.

  • Estamos leyendo (y nos va gustando mucho) el libro Data Analysis for the Life Sciences de Irizarry y Love.

  • Aunque no sea específicamente sobre R, quisiera recomendar el libro Bioinformatics Data Skills de V. Buffalo, como un excelente punto de inicio a la bioinformática.

Cursografía / Webografría

Cursos.

Ya hay literalmente cientos de cursos (ver aquí, aquí y aquí).

Tutoriales y Blogs.

Hablemos de p-valores.

Una llamada de atención.

In February 2014, George Cobb, Professor Emeritus of Mathematics and Statistics at Mount Holyoke College, posed these questions to an ASA discussion forum:

Q: Why do so many colleges and grad schools teach p = 0.05?

  A: Because that’s still what the scientific community and journal editors use.

Q: Why do so many people still use p = 0.05?

  A: Because that’s what they were taught in college or grad school.

El abuso del p-valor.

  1. P-values can indicate how incompatible the data are with a specified statistical model.
  2. P-values do not measure the probability that the studied hypothesis is true, or the probability that the data were produced by random chance alone.
  3. Scientific conclusions and business or policy decisions should not be based only on whether a p-value passes a specific threshold (e.g. 0.05).
  4. Proper inference requires full reporting and transparency.
  5. A p-value, or statistical significance, does not measure the size of an effect or the importance of a result.

P-valor y tamaño del efecto

  • La cuarta advertencia es una nueva llamada a la reproducibilidad, y un aviso contra el p-hacking.

  • Las tres primeras se deben simplemente a conceptos estadísticos mal entendidos (o explicados, nostra culpa...). O sea, que se deben a que la gente piensa cosas que son erróneas.

  • Pero vamos a fijarnos en particular en el quinto punto, que nos parece más sutil, porque se debe a cosas en las que muchas veces la gente no piensa. Para ello vamos a valernos de una ventaja fundamental de usar R.

  • R es, además, un lenguaje de programación. Y en particular permite hacer muy fácilmente simulaciones que nos permiten probar in silico nuestras ideas sobre el comportamiento estadístico de los datos. Veámoslo. Puedes ejecutar el código a la vez, para comparar con nuestros resultados.

Para empezar nuestra simulación vamos a elegir una semilla del generador de números aleatorios (yo siempre uso el año para que nadie me acuse de seed-hacking).

set.seed(2017)

Ahora vamos a construir dos poblaciones, ambas con distribuciones normales, que representan las alturas medias de los hombres en dos ciudades A y B.

poblacion1 = rnorm(n = 1000000, mean = 171.1, sd = 1)

poblacion2 = rnorm(n = 1000000, mean = 171.3, sd = 1)

Nuestra hipótesis alternativa es que las alturas son distintas. Para empezar un estudio comparativo tomamos muestras de tamaño \(n = 10\) de ambas poblaciones.

n = 10

muestra1 = sample(poblacion1, n, replace = TRUE)
muestra2 = sample(poblacion2, n, replace = TRUE)

Ahora hacemos, correctamente, una prueba t (bilateral) para obtener el p-valor del contraste:

contraste01 = t.test(x = muestra1, y = muestra2, 
                     alternative="two.sided", var.equal = TRUE)

contraste01$p.value 
## [1] 0.7123433

Hemos supuesto que las varianzas de ambas poblaciones son iguales (de hecho lo son), pero cambiar eso no afectaría al resultado. Como ves el p-valor 0.71 indica que no tenemos evidencias para rechazar la hipótesis nula.

Pero nosotros, infatigables, seguimos creyendo que hay una diferencia y tomamos una muestra más grande, repitiendo el análisis:

n = 100

muestra1 = sample(poblacion1, n, replace = TRUE)
muestra2 = sample(poblacion2, n, replace = TRUE)

contraste02 = t.test(x = muestra1, y = muestra2, 
                     alternative="two.sided", var.equal = TRUE)

contraste02$p.value
## [1] 0.7649759

¿Todavía no? No pasa nada, el premio es de los que perseveran y se esfuerzan. Tomemos una muestra aún más grande y repitamos el análisis.

n = 1000

muestra1 = sample(poblacion1, n, replace = TRUE)
muestra2 = sample(poblacion2, n, replace = TRUE)

contraste03 = t.test(x = muestra1, y = muestra2, 
                     alternative="two.sided", var.equal = TRUE)

contraste03$p.value
## [1] 0.0004295796

Ahora el p-valor es 4.3 × 10-4. ¡Altamente significativo! Si ya lo decía yo... ya podemos ir a la prensa y decir: "un estudio basado en una amplia muestra confirma que los habitantes de A son significativamente más altos que los de B..."

¿Y entonces?

  • El resultado se debe, naturalmente, a que mientras haya una diferencia entre las medias (¿y cuándo no la habría?), si tomamos una muestra suficientemente grande, siempre acabaremos por "detectar" la diferencia, por irrelevante que sea esa diferencia desde el punto de vista científico. Es una advertencia que no debe olvidarse nunca, pero especialmente en la era del Big Data.

  • ¿Remedios? Una buena manera de evitar esto es desplazando el foco del p-valor a la estimación del efecto real, por ejemplo mediante un intervalo de confianza. El intervalo de confianza para la diferencia de medias en el último contraste muestra que hemos detectado una diferencia máxima en altura de 0.24cm (¡irrelevante!):

contraste03$conf.int
## [1] -0.23773123 -0.06783195
## attr(,"conf.level")
## [1] 0.95

Algunas ideas sobre gráficos.

Por supuesto, el problema del anterior contraste se podría detectar con un gráfico bien hecho. Incluso con el gráfico de barras típico, denostado por muchos estadísticos:

plot of chunk unnamed-chunk-11

La diferencia es tan insignificante que las barras de error apenas se distinguen...

Pero entonces, ¿por qué a los estadísticos no les gustan los gráficos de barras como este? Ese es nuestro siguiente tema de discusión.

Lo primero que debemos reconocer es que ese tipo de gráficos son extremadamente populares. Aparecen en casi cualquier número de tu revista favorita. Un ejemplo reciente de Science:

Y sin embargo muchos estadísticos detestan estos gráficos y los denominan peyorativamente "dynamite plots". Para entender una de las razones fundamentales por las que no nos gustan vamos a hacer otro ejercicio de simulación (basado en Harrell2008 ).

Vamos a trabajar con los datos de un fichero llamado dynaPLotData.csv que puedes descargar de este enlace. Los descargamos en el directorio de trabajo (no lo olvides). Y para cargar esos datos en R ejecutamos:

datosEstudio = read.table("dynaPLotData.csv", header = TRUE, sep = ",")

Veamos el comienzo de la tabla de datos:

head(datosEstudio)
##     plasma   group   when gender
## 1 42.45219 placebo before      M
## 2 46.12612 placebo before      F
## 3 34.88628 placebo before      F
## 4 42.48579 placebo before      M
## 5 44.83357 placebo before      F
## 6 33.98735 placebo before      M

Como puedes ver, cada fila contiene una observación de cinco variables:

  1. La medida de cierta cantidad en plasma, que es la variable objeto de interés en este estudio.
  2. Un factor que identifica si esa observación corresponde al grupo tratado o al placebo.
  3. Otro factor when que nos dice si el valor de plasma es previo o posterior al tratamiento/placebo.
  4. El género del paciente, un factor con niveles M y F.

La columna when, toma estos valores (con las frecuencias que se indican):

table(datosEstudio$when)
## 
##  after before 
##     36     36

Para empezar vamos a seleccionar las observaciones hechas antes del tratamiento, usando subset para hacer esto:

antes = subset(datosEstudio, when=="before")

Veamos la composición de los grupos:

table(antes$group)
## 
##   placebo treatment 
##        14        22

Vamos a hacer una prueba \(t\) para ver si hay diferencia significativa de los valores de plasma entre ambos grupos (siempre dentro del conjunto de datos antes).

El resultado es:

t.test(plasma ~ group, data = antes, alternative="two.sided")
## 
##  Welch Two Sample t-test
## 
## data:  plasma by group
## t = 0.11886, df = 28.166, p-value = 0.9062
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.246574  3.646648
## sample estimates:
##   mean in group placebo mean in group treatment 
##                42.37184                42.17180

Y claramente muestra que no hay evidencia a favor de una diferencia significativa (de lo contrario tendríamos un problema, claro). Veamos el diagrama de columnas con barras de error correspondiente.

plot of chunk unnamed-chunk-19

El diagrama parece confirmar el resultado de la prueba \(t\). Este diagrama está hecho con la librería ggplot de R. No vamos a explicar el código, por falta de tiempo, pero si vas a usar R en serio en algún momento tendrás que aprender a usar ggplot.

A continuación vamos a hacer un análisis similar para las observaciones después de aplicar el tratamiento:

despues = subset(datosEstudio, when=="after")

La prueba \(t\) muestra una diferencia muy significativa entre ambos grupos:

t.test(plasma ~ group, data = despues, alternative="two.sided")
## 
##  Welch Two Sample t-test
## 
## data:  plasma by group
## t = -3.7311, df = 21.087, p-value = 0.001227
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -136.00708  -38.67083
## sample estimates:
##   mean in group placebo mean in group treatment 
##                 52.4068                139.7458

Y el diagrama correspondiente a los grupos tras el tratamiento parece enviar el mismo mensaje:

plot of chunk unnamed-chunk-22

En muchas publicaciones la información que contienen estos gráficos es parecida a la de este ejemplo.

Pero ahora vamos a hacer otra cosa. Vamos a añadir los datos reales a ese diagrama.

plot of chunk unnamed-chunk-23

¡Sorpresa! La respuesta al tratamiento depende de la variable género de tal manera que la media de las mujeres se ha disparado y la de los hombres se ha hundido.

Este ejemplo ha tratado de ilustrar una de las limitaciones de la forma habitual de usar los gráficos con barras: no proporcionan apenas información sobre la distribución de los datos. Hay alternativas mejores, desde los dot plots (diagramas de puntos) para muestras pequeñas a los boxplots o violin plots para muestras más grandes. Un ejemplo de estos últimos en el último número de Science:

Por cierto, el artículo de Science que hemos citado al principio de esta discusión incluye varios ejemplos en los que se ha incorporado la muestra al diagrama.

Un violin plot de los datos del grupo tratamiento en el anterior ejemplo habría mostrado claramente que estamos ante una distribución bimodal (a la izquierda aparece el grupo placebo):

library(vioplot)
vioplot(studyA$plasma[studyA$group=="placebo"], wex = 0.5, col="orange")
vioplot(studyA$plasma[studyA$group=="treatment"], wex = 0.5, col="orange")

plot of chunk unnamed-chunk-24

Confiamos en que esta discusión al menos te haya puesto en guardia frente a esos gráficos. Para ampliar esta discusión puedes consultar por ejemplo la información expuesta aquí y aquí.

Pero ya que estamos:

  • ¿qué significan las barras que aparecen en ellos? ¿Son errores estándar (se)? ¿Son intervalos de confianza?

  • Y en cualquier caso, ¿entendemos bien la razón teórica para hacer contrastes comparando intervalos?

Creemos que en general hay bastante confusión sobre este último punto, así que vamos a tratar de abordarlo en el siguiente apartado.

Usando los intervalos de confianza para hacer contrastes.

Un problema que afecta también a los gráficos con barras de error es que sugieren el uso de intervalos de confianza, uno por muestra, para hacer contrastes de hipótesis. Y esa idea no está exenta de problemas. Veamos un ejemplo.

El fichero IntervalosVsContrastesNoSolapan.csv contiene dos columnas, cada una con datos de una muestra. Vamos a usar esos datos para contrastar la diferencia de medias entre ambas muestras. Empieza descargándote los datos al directorio de trabajo y leyéndolos con este comando.

datos = read.table("IntervalosVsContrastesNoSolapan.csv", header = TRUE)

Vamos a aprovechar que la función t.test nos permite obtener fácilmente intervalos de confianza para la media de cada muestra por separado.

Para la primera:

test1 = t.test(datos$muestra1, conf.level = 0.95)
test1$conf.int
## [1] 133.3888 135.6112
## attr(,"conf.level")
## [1] 0.95

Y para la segunda:

test2 = t.test(datos$muestra2, conf.level = 0.95)
test2$conf.int
## [1] 138.3714 141.2286
## attr(,"conf.level")
## [1] 0.95

Representando gráficamente ambos intervalos:

No solapan. La conclusión parece ser que podemos rechazar la igualdad de medias al 95%.

Si hacemos una prueba \(t\) para la diferencia de medias usando conjuntamente las dos muestras tenemos:

t.test(datos$muestra1, datos$muestra2, alternative = "two.sided")
## 
##  Welch Two Sample t-test
## 
## data:  datos$muestra1 and datos$muestra2
## t = -5.8105, df = 186.69, p-value = 2.645e-08
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -7.099427 -3.500573
## sample estimates:
## mean of x mean of y 
##     134.5     139.8

Ese p-valor tan pequeño y el intervalo de confianza que no contiene al 0 confirman nuestra idea inicial.

Ahora vamos a trabajar con un segundo fichero, llamado IntervalosVsContrastesNoSolapan.csv. Hagamos las mismas operaciones. Descargamos y leemos el fichero:

datos = read.table("IntervalosVsContrastesSolapan.csv", header = TRUE)

Calculamos intervalos de confianza para la media de cada muestra por separado.

test1 = t.test(datos$muestra1, conf.level = 0.95)
test1$conf.int
## [1] 131.4046 137.5954
## attr(,"conf.level")
## [1] 0.95
test2 = t.test(datos$muestra2, conf.level = 0.95)
test2$conf.int
## [1] 136.3871 143.2129
## attr(,"conf.level")
## [1] 0.95

Ahora, como muestra la figura, hemos obtenido intervalos de confianza que sí solapan.

¿Rechazamos la hipótesis alternativa de medias diferentes?

No tan rápido...

Vamos a hacer, como antes, una prueba \(t\) para la diferencia de medias usando ambas muestras:

t.test(datos$muestra1, datos$muestra2, alternative = "two.sided")
## 
##  Welch Two Sample t-test
## 
## data:  datos$muestra1 and datos$muestra2
## t = -2.2824, df = 196.14, p-value = 0.02354
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -9.8794255 -0.7205745
## sample estimates:
## mean of x mean of y 
##     134.5     139.8

Un p-valor menor que 0.05 y, en consecuencia, un intervalo de confianza que no contiene al origen. Así que rechazamos la hipótesis nula de igualdad de media. ¿En qué quedamos?

El problema es, por decirlo de una manera sencilla, que los intervalos de confianza tratan a cada variable por separado. Y al hacerlo, sobrevaloran la probabilidad de que las dos medias se parezcan.

Puede venir bien pensarlo de esta manera: para que, a partir de dos muestras de estas poblaciones, concluyamos que las dos medias se parecen mucho, la media poblacional de la primera debe haber caído en la parte más alta de su intervalo de confianza (o más arriba aún) y la media poblacional de la segunda debe caer en la parte baja de su intervalo (o más abajo aún).

El contraste de diferencia de medias tiene esto en cuenta al calcular las probabilidades. Pero al fijarnos en intervalos de confianza separados, asumimos que esos dos sucesos, de por sí poco probables, ocurren simultáneamente. Y eso causa nuestra sobreestimación de la probabilidad. Si además hacemos comparaciones múltiples, el efecto se multiplica, como sabemos.

Intervalos de confianza por bootstrap.

  • Muchos métodos de la Estadística clásica (intervalos de confianza, contrastes de hipótesis) asumen que las variables que estudiamos se distribuyen de forma normal. Entre otras cosas, eso significa que los intervalos de confianza para la media son simétricos respecto a la media muestral. Pero a menudo encontramos muestras muy asimétricas, que no justifican la simetría del intervalo.

  • El aumento de la capacidad de cómputo ha propiciado el desarrollo de métodos no paramétricos de construcción de intervalos de confianza basados en el remuestreo, como el bootstrap. Aquí vamos a usar ese método para calcular un intervalo de confianza de los datos contenidos en el fichero (basado de Crawley2011 , pág. 47): skewdata.csv

plot of chunk unnamed-chunk-32

Empezamos leyendo esos datos:

x = read.table(file = "skewdata.csv", header = TRUE)[, 1]
  • Ahora vamos a explorar tamaños muestrales entre \(n = 5\) y \(n = 40\).
  • Para cada tamaño construiremos \(10000\) remuestreos aleatorios de esa muestra.
  • Para cada remuestreo calculamos la media muestral, obteniendo una lista de 10000 medias muestrales.
  • Y representamos el intervalo del primer al tercer cuartil de esa lista con 10000 medias correspondientes a muestras de tamaño \(n\).
  • En la gráfica representamos esos intervalos (bootstrap intervals) en naranja, la media en azul y en rojo representamos los intervalos clásicos usando la \(t\) de Student. Fíjate en que para muestras grandes no hay apenas diferencia. Pero en muestras pequeñas el intervalo bootstrap refleja mucho mejor la asimetría de los datos.

plot of chunk bootstrap

El código R que hemos usado es muy sencillo, con apenas un bucle for. Casi todo lo demás es para los gráficos.

plot(c(0, 40), c(5,11), type="n", xlab="Sample size", ylab="Confidence interval") # Creamos la "caja" del gráfico.

for (k in seq(5, 40, 1)){ # Este bucle recorre los tamaños muestrales
  a =  numeric(10000) # el vector a almacenará las medias muestrales
  for (i in 1:10000){ # este es el bucle de remuestreo (bootstrap)
    a[i] = mean(sample(x, k, replace=T)) # generamos un remuestreo con reemp. y calculamos su media
    }
  points(c(k,k), quantile(a, c(.025,.975)), type="b", col = "orange", lwd= 3) # dibujo del intervalo bootstrap de este tamaño muestral
}

# el siguiente bloque de código genera una banda con los intervalos clásicos correspondientes a esas muestras.
xv = seq(5, 40, 0.1) 
yv = mean(x) - qt(0.975, xv) * sqrt(var(x) / xv)
lines(xv, yv, lty=2, col= "red", lwd = 2)
yv = mean(x) + qt(.975, xv) * sqrt(var(x) / xv)
lines(xv, yv, lty=2, col= "red", lwd = 2)

abline(h = mean(x), col="blue", lwd=2) # añadimos una línea horizontal en la media

De vuelta al laboratorio...

...diseño de experimentos.

Potencia y tamaño muestral.

Empezando por lo más básico, querríamos volver a la noción de potencia de un contraste. Recordemos que en un contraste se pueden dar esencialmente estas cuatro situaciones:

\(H_0\) es cierta \(H_0\) falsa
Rechazamos \(H_0\) error tipo I (\(\alpha\)) correcto
Rechazamos \(H_a\) correcto error tipo II (\(\beta\))

Hemos indicado entre paréntesis las probabilidades \(\alpha\) y \(\beta\) de cometer cada tipo de error. Entonces la potencia es \[ 1 - \beta\] Es decir, la probabilidad de no cometer un error de tipo II. Dicho de forma más sencilla:

La potencia es la probabilidad de rechazar \(H_0\) cuando \(H_0\) es, de hecho, falsa.

Vamos a fijarnos en un caso muy sencillo. Aunque las situaciones prácticas suelen ser mucho más complicadas los mismos principios generales se aplican a muchas de ellas.

Cuando usamos una muestra para contrastar si la media \(\mu\) de una variable normal es o no es igual a un cierto valor \(\mu_0\) la potencia del contraste es aproximadamente proporcional a:

\[\dfrac{\delta\,\sqrt{n}\,\alpha}{\sigma}\] donde:

  1. \(n\) es el tamaño de la muestra. A más muestra, más potencia.
  2. \(\delta\) es la diferencia entre \(\mu\) y \(\mu_0\). Es decir, el efecto que esperamos detectar. Cuanto mayor sea el efecto, mayor potencia.
  3. \(\alpha\) es el nivel de significación del contraste, que solemos fijar en \(0.05\). Cuanto más pequeño sea, menor potencia. No podemos tener lo mejor de los dos mundos.
  4. \(\sigma\) es la desviación típica, que indica la dispersión de la población. La solemos estimar mediante \(s\), a menudo procedente de estudios piloto o previos. Y a mayor dispersión, menor potencia.

Las ecuaciones de potencia suelen tener un aspecto similar a ese. Y se usan a menudo para determinar el tamaño muestral necesario para poder detectar un efecto de tamaño dado, con niveles de significación y potencia dados.

Ejemplo: vamos a averiguar el tamaño muestral \(n\) necesario para realizar un contraste unilateral de hipótesis con un nivel de confianza \(nc = 0.99\) (es decir, \(\alpha = 0.01\)), y una potencia \(1 − \beta = 0.80\), que sea capaz de detectar una diferencia en las medias (efecto) al menos igual a \(\delta = 0.1\). La cuasidesviación típica procedente de un estudio piloto es \(s = 0.5\). Suponiendo que las medias muestrales son normales queremos que sea: \[ 1-\beta=P\left(Z> z_{\alpha}-\dfrac{\delta}{\dfrac{s}{\sqrt{n}}}\right). \] Sustituyendo los datos se obtiene: \[ \dfrac{\sqrt{n}}{5} = z_{0.01} - z_{0.80} \] donde los \(z_p\) son cuantiles de la normal. Es decir, \(n \approx 250.90\).

Nuestro cálculo está basado en una estimación muy simplista de las probabilidades. En R existe una función que realiza un cálculo mucho más ajustado.

power.t.test(delta = 0.1, sd = 0.5, sig.level = 0.05, 
             power = 0.80, type="one.sample", alternative="one.sided")
## 
##      One-sample t test power calculation 
## 
##               n = 155.9257
##           delta = 0.1
##              sd = 0.5
##       sig.level = 0.05
##           power = 0.8
##     alternative = one.sided

Esta función puede usarse para determinar uno de los valores a partir de los restantes.

Por ejemplo, si nuestra muestra sólo es de tamaño \(n = 100\),

power.t.test(delta = 0.1, sd = 0.5, sig.level = 0.05, 
             n = 100, type="one.sample", alternative="one.sided")
## 
##      One-sample t test power calculation 
## 
##               n = 100
##           delta = 0.1
##              sd = 0.5
##       sig.level = 0.05
##           power = 0.6336178
##     alternative = one.sided

este cálculo nos informa de que la potencia ha caído al 63%.

En R hay varias de estas funciones para hacer cálculos de potencias y tamaños muestrales. Pero a menudo están dispersas por librerías (ver, por ejemplo, la librería pwr), y hay que entender bien su uso.

A veces es más sencillo usar alguno de los servicios web que proporcionan una interfaz simplificada para este tipo de cálculos.

Por ejemplo hay uno disponible en esta dirección (ver la pestaña Sample Size):

http://www.obg.cuhk.edu.hk/ResearchSupport/StatTools/ResourceIndex_Categories.php

¿Qué método estadístico usar?

  • Hemos dejado casi para el final otra de las preguntas frecuentes que los experimentadores dirigen a un estadístico. En general la respuesta, nos tememos, no será breve. El estadístico querrá detalles sobre el diseño experimental y tomarse su tiempo de reflexión para contestar.

  • Desafortunadamente, esa pregunta llega muchas veces cuando ya se han recogido los datos. Y a menudo se repite este mantra: la estadística no puede arreglar un mal diseño.

  • Que tiene su cara positiva en este otro: un experimento bien diseñado es más fácil de analizar y emplea menos recursos.

  • Siempre es una buena idea aprender más sobre diseño experimental. Y para eso "las dos orillas" tenemos que hablar más.

Muchas gracias por la atención.