blog en español.

¿De donde vienen los valores p? conceptos fundamentales y un acercamiento por simulación.

el valor-p es la probabilidad de las colas calculado de la distribución muestral de un estadístico basado en muestras. Esta distribución muestral depende del tamaño de la muestra el estadistico a utilizar y la población a la cual se le toma la muestra.

¿Qué es P-valor?

El valor p es un concepto dado por la estadística frecuentista que se enseña en el colegio o las universidades. En este tipo de estadística se asume hipotéticamente una población infinita donde se considera que los datos constituyen una muestra aleatoria. La idea de este planteamiento es crear proposiciones acerca de los fenómenos de la población así como respecto a la muestra encontrar una distribución apropiada. Como los datos son tratados como una muestra aleatoria, cualquier cantidad calculada de los datos debe ser tratada como una variable aleatoria, estas cantidades deben seguir una distribución de probabilidad. Por ejemplo si se toman multiples muestras de una misma población y se calcula la media de cada una de estas muestras algunos valores vas a ser más comunes que otros. Los valores que se puedan calcular de las muestras son llamados estadísticos.

Otro concepto importante es la prueba de hipótesis del valor p. Una hipótesis nula es una proposición acerca de las diferencias entre fenomenos o grupos, la hipotesis nula representa una posición escéptica, el objetivo de una prueva significativa de hipótesis nula es desafiar la hipótesis nula con los datos observados.

Este procedimiento requiere la construcción de un modelo estadístico para describir la población aleatoria de la cual los datos podrían haber sido muestreados, asumiendo que la hipótesis nula es verdadera. Luego, se calcula la probabilidad de que ocurran los datos observados asumiendo que se tomaron muestras de esa población. Cuanto más improbables son los datos, menos probable se vuelve la hipótesis nula. Sin embargo, comparar los datos observados con todos los datos alternativos posibles es demasiado difícil, por lo que el problema a menudo se simplifica al comparar una estadística calculada a partir de los datos con todos los valores posibles de la misma estadística que podrían haberse obtenido de la población. En la práctica, esto significa que necesitamos construir una distribución muestral y evaluar cuán improbable es la estadística observada dentro de esa distribución muestral.

La complicación final es que las cantidades que usamos como estadística son a menudo medidas continuas de los datos (incluso cuando los datos están compuestos de valores discretos) pero las leyes básicas de la probabilidad nos dicen que la probabilidad de cualquier valor particular de una variable aleatoria continua es exactamente cero. Por lo tanto, calcular la probabilidad de un valor estadístico particular en una distribución muestral es inútil. Como alternativa, se puede calcular la probabilidad de todos los valores más extremos que el observado (es decir, el área sombreada en la figura anterior). Cuanto menor sea esta probabilidad, más improbable será el valor.¡Y esta probabilidad es el valor p!

Dado que todo esto fue muy abstracto, pasaremos por un ejemplo clásico: la prueba t de una muestra. Ilustraré con un enfoque computacional cómo la distribución muestral surge de la población y cómo se relaciona con la estadística muestral. Elegí R como el lenguaje de programación para este artículo ya que es un lenguaje común para las estadísticas, pero el mismo enfoque podría haberlo hecho en Matlab/Octave, Python o Julia.

Ejemplo para la prueba de media de una muestra (prueba t de una muestra)

Imagine que tenemos un pequeño conjunto de datos.

muestra=c(-0.88,0.58,2.83,0.57,1.49,1.69,0.43,0.58,-0.88,0.77,2.50,
          0.45,1.24,0.28,1.95,2.19,0.18,-0.81,1.49,1.91,0.99,0.96,0.91,
          0.72,-0.31,0.28,1.21)

La mayoría de los valores son positivos y es posible que deseemos probar si el efecto subyacente es realmente positivo. Asumiremos que estos datos podrían tratarse como una muestra aleatoria de una población normal con una media positiva (es decir \(\mu>0\)). Como se describió anteriormente, establecemos como hipótesis nula (\(H_0\)) la opción más conservadora, que en el ejemplo sería \(\mu = 0\), mientras que la hipótesis alternativa o alternativa (\(H_1\)) representa el caso positivo (\(\mu> 0\)). Formalmente, esto se escribe como:

\[H_0:\mu=\mu_0\] \[H_1:\mu>\mu_0\]

Donde \(\mu_0=0\) en este caso. Suponiendo que el procedimiento de muestreo asegura la independencia estadística y dados los supuestos anteriores, podríamos usar el estadístico:

\[\frac{(\bar{x}-\mu_0)\sqrt{n}}{s} \]

donde \(\bar{x}\) es la media, s es la desviación estándar y n es el tamaño de la muestra. La razón para usar esta estadística es que se sabe que la distribución muestral es una distribución t de Student. La prueba de esto la dio William Gosset en 1908 (firmado con el seudónimo de Student). Si desea comprobar lo que se necesita para realizar dicha prueba, eche un vistazo al documento original.

En lugar de pasar por la demostración matemática, usaré un procedimiento basado en simulación que encuentro más intuitivo y más fácil de seguir, y también es muy general, mientras produce los mismos resultados. Efectivamente, construiremos una aproximación de Monte Carlo a la distribución muestral del estadístico t. El procedimiento de cálculo es el siguiente:

1.) Estime todos los parámetros de la población de la muestra, excepto el parámetro que ya está fijado por la hipótesis nula. En este ejemplo, esto significa que estimamos la desviación estándar de la población normal \(\sigma\) de la muestra, mientras fijamos la media \(\sigma\) en 0.

2.) Genere N muestras aleatorias de la población de tamaño n utilizando el método de Monte Carlo. En este caso, n es igual a 100 y N puede ser tan grande como desee (cuanto más grande, mejor).

3.) Para cada muestra generada en el paso 2, calcule la estadística. En este caso, calculamos \(\frac{(\bar{x}-\mu_0)\sqrt{n}}{s}\) para cada muestra.

4.) El estadístico N calculado en el paso 3 es una muestra grande de valores de la distribución muestral. El valor p se puede calcular como la fracción de esta muestra que excede (o es menor, dependiendo de la cola que se esté probando) que el valor de la estadística en los datos reales.

Intermezzo técnico

Para poner en práctica estas ideas con R, necesitamos cargar un par de bibliotecas:

  • furrr que facilita la ejecución paralela de código en múltiples procesos basados en el futuro del paquete, porque la vida es demasiado corta para esperar los cálculos.

  • distr6 para crear y manipular distribuciones de probabilidad.

for(lib in c("furrr", "distr6")) {
  library(lib, character.only = TRUE)
}
plan(multiprocess)

Estimación de parámetros de la población

Existen diferentes métodos para estimar los parámetros de un modelo estadístico o distribución de una muestra. Para el caso de la desviación estándar de una distribución normal, el enfoque canónico es usar la desviación estándar con la corrección de Bessel, también conocida como desviación estándar de la muestra (es decir, usando n - 1 en lugar de n). Resulta que este es el comportamiento predeterminado de la función sd de R, por lo que el estimador es simplemente:

\(\hat{\sigma} = sd(muestra)\)

El otro parámetro de la distribución muestral es dado por la hipótesis nula

Construyendo la distribución muestral

Primero, construyamos la distribucion normal que describe la población:

dsta=sd(muestra)  #desviación estandar de la muestra.
mu0=0
set.seed(2019) 
poblacion = rnorm(10000,mu0,dsta) #Población generada con la desvaición muestral
                                #y con media dada por la hipótesis nula.

Luego, necesitamos una función que calcule el estadístico \(\frac{(\bar{x}-\mu_0)\sqrt{n}}{s}\) para una muestra dada

calc_statistic = function(x, mu0 = 0) {
  (mean(x) - mu0)/(sd(x)/sqrt(length(x)))
}

Luego, generamos N muestras de la población y calculamos la estadística para cada una de ellas. Cuanto mayor sea N, más preciso será el cálculo del valor p. Resulta que calcular valores p precisos es computacionalmente difícil, por lo que, para asegurarnos de que nuestros cálculos sean precisos, generemos 50000 muestras:

N = 5e4
set.seed(2019)
distr.mues = future_map_dbl(1:N, ~ calc_statistic(sample(poblacion,size = 27,replace =FALSE))) #50 000 muestras de tamaño 27.

Comparemos la con la distribución teórica y la estadística de muestra observada:

plot(density(distr.mues,), xlab = "Estadística t", main = "Teórica vs estadística muestral", xlim = c(-10,10)) # Histograma de distr.mues
curve(dt(x,26), -10, 10, add = T, col = 2, n = 1e3) #Distribución teórica t(n-1)
abline(v = calc_statistic(muestra), col = 3) # Estadística observada

Podemos ver que, de hecho, la distribución muestral es la distribución t de Student con 26 grados de libertad. Además, el valor observado para el estadístico se encuentra en un área de baja probabilidad en la distribución muestral, lo que sugiere que es poco probable que este resultado ocurra bajo la hipótesis nula.

Calcule el valor p.

En este ejemplo, el valor p se define como la probabilidad de que una muestra del modelo de hipótesis nula conduzca a una estadística igual o mayor que la observada. Calcular este valor a partir de la distribución muestral es tan simple como calcular la fracción de casos en los que se cumple esta condición:

sum(distr.mues >= calc_statistic(muestra))/N
## [1] 0

Podemos comparar este valor con el obtenido de la distribución t:

1 - pt(calc_statistic(muestra), 26)
## [1] 4.407106e-05

y con los resultados de llamar a la función t.test en R (que debería ser el mismo que para la distribución t a menos que haya aprendido la teoría incorrecta …):

t.test(muestra, alternative = "greater")$p.val
## [1] 4.407106e-05

Podemos ver que la estimación de Monte Carlo es muy similar al valor analítico. Tenga en cuenta que todavía hay un pequeño error numérico a pesar de la gran cantidad de muestras de Monte Carlo. Hay métodos más avanzados para refinar una estimación del valor p de Monte Carlo, que se vera en publicaciones futuras.

¿Por qué no utilizar la media muestral como estadística?

Usar la media muestral \(\bar{x}\) como estadística en lugar de \(\frac{(\bar{x}-\mu_0)\sqrt{n}}{s}\) habría sido más intuitivo, ya que estamos haciendo una pregunta sobre la media, no una función complicada de la misma. Con el enfoque basado en simulación, esto es trivial (o al menos no más difícil que el anterior). Creemos la distribución muestral para la media muestral:

N = 5e4
set.seed(2019)
distr.mues2 = future_map_dbl(1:N, ~ mean(sample(poblacion,size = 27,
                           replace =FALSE))) #50 000 muestras de tamaño 27
## Warning: UNRELIABLE VALUE: Future ('<none>') unexpectedly generated random
## numbers without specifying argument 'seed'. There is a risk that those random
## numbers are not statistically sound and the overall results might be invalid.
## To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe random
## numbers are produced via the L'Ecuyer-CMRG method. To disable this check, use
## 'seed=NULL', or set option 'future.rng.onMisuse' to "ignore".
## Warning: UNRELIABLE VALUE: Future ('<none>') unexpectedly generated random
## numbers without specifying argument 'seed'. There is a risk that those random
## numbers are not statistically sound and the overall results might be invalid.
## To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe random
## numbers are produced via the L'Ecuyer-CMRG method. To disable this check, use
## 'seed=NULL', or set option 'future.rng.onMisuse' to "ignore".
## Warning: UNRELIABLE VALUE: Future ('<none>') unexpectedly generated random
## numbers without specifying argument 'seed'. There is a risk that those random
## numbers are not statistically sound and the overall results might be invalid.
## To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe random
## numbers are produced via the L'Ecuyer-CMRG method. To disable this check, use
## 'seed=NULL', or set option 'future.rng.onMisuse' to "ignore".
## Warning: UNRELIABLE VALUE: Future ('<none>') unexpectedly generated random
## numbers without specifying argument 'seed'. There is a risk that those random
## numbers are not statistically sound and the overall results might be invalid.
## To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe random
## numbers are produced via the L'Ecuyer-CMRG method. To disable this check, use
## 'seed=NULL', or set option 'future.rng.onMisuse' to "ignore".

que se ve como:

#Histograma de la distribución muetsral para la media muestral.
plot(density(distr.mues2, bw = "SJ"), main = "", xlab = "sample mean",xlim = c(-0.86,0.9)) 
 abline(v = mean(muestra), col = 2) #Estadística observada.

y el valor p es calculado como en el proceso anterior

sum(distr.mues2 >= mean(muestra))/N
## [1] 0

Observe que el valor p es diferente al anterior. Esto no es incorrecto, dado que dada la misma muestra y modelo de población, el valor p dependerá de la estadística que se calcule. Quizás me pregunte cuál sería el valor analítico en este caso. Hasta donde yo sé, no hay ninguno, la razón es que hasta ahora no ha sido posible derivar la distribución muestral para la media muestral. Los únicos resultados conocidos son para una muestra hipotética de tamaño infinito (cuando la distribución muestral se vuelve Normal, independientemente de la población). Si bien estos resultados pueden usarse razonablemente para muestras grandes, claramente no se aplica a un tamaño de muestra de 27.

Observaciones finales

Como se ha mostrado en este artículo, para calcular el valor p de un estadístico necesitamos conocer su distribución muestral, que depende del modelo asumido para la población. Tradicionalmente, se tenía que derivar la distribución muestral matemáticamente (si era posible) para cada combinación de modelo estadístico y poblacional. Esto restringe en gran medida el número de escenarios que pueden ser manejados por las pruebas estadísticas tradicionales, lo que lleva a la temida pregunta: “¿Qué prueba estadística debo usar para mis datos?”. Hay docenas de árboles de decisión que se supone que debe seguir el practicante de estadística.

Hoy en día, puede utilizar un enfoque de simulación como el que se presenta en este artículo. Con este enfoque, siempre puede calcular el valor p para cualquier modelo estadístico y poblacional, incluso si nadie ha probado esa combinación en particular antes. Restringirse a los métodos analíticos tenía sentido a principios del siglo XX cuando se desarrollaron muchas de estas pruebas, ¡pero en el siglo XXI se pueden usar computadoras! En la próxima publicación, mostraré cómo aplicar este método a datos del mundo real que no estaban (y no podían) distribuirse normalmente. ¡Manténganse al tanto!