Actividad realizada como practica academica tomada y traducida de https://www.r-bloggers.com/2019/08/where-do-p-values-come-from-fundamental-concepts-and-simulation-approach/

Dado que los datos se interpretan como una muestra aleatoria, cualquier cantidad calculada a partir de los datos (por ejemplo, la media aritmética de los valores observados, x¯) puede tratarse como una variable aleatoria. Es decir, cada muestra aleatoria de la población producirá un valor diferente de la misma cantidad (por ejemplo, diferentes medias muestrales). Además, significa que estas cantidades seguirán una distribución de probabilidad. Por ejemplo, si tomamos varias muestras aleatorias de la misma población y calculamos la media de cada una de estas muestras, algunos valores de la media de la muestra serán más comunes que otros. Esta distribución es lo que llamamos distribución muestral y la cantidad derivada de las muestras se llama estadística. En las estadísticas frecuentistas, cualquier enunciado de probabilidad (por ejemplo, valores p, intervalos de confianza) siempre se refiere a una estadística y su distribución muestral asociada.

El otro concepto importante detrás del valor p es el de la prueba de hipótesis. Una hipótesis nula es un enunciado sobre diferencias entre grupos o fenómenos. La hipótesis nula representa una posición escéptica predeterminada (es decir, sin diferencia). El objetivo de la prueba significativa de hipótesis nula (NHST) es desafiar la hipótesis nula dados 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 (hay muchos recursos en línea que explican esta idea, por ejemplo, esta). 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.

\[\color{blue}{Ejemplo~para~la~prueba~de~media~de~una~muestra~(prueba~t~de~una~muestra)}\] Imaginemos que tenemos un pequeño conjunto de datos como:

sample = c(1.52, 5.24, -0.23, 2.47, 2.63)
sample
## [1]  1.52  5.24 -0.23  2.47  2.63

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, μ> 0). Como se describió anteriormente, establecemos como hipótesis nula (H0) la opción más conservadora, que en el caso sería μ = 0, mientras que la hipótesis contestataria o alternativa (H1) representa el caso positivo (μ> 0). Formalmente, esto se escribe como:

\[H_0:\mu_{0}\geq\mu_{0}\] \[H_0:\mu_{0}<\mu_{0}\] Donde μ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 t, (x¯ – μ0) / (s / n −− √), donde x¯ es la media de la muestra, s es la desviación estándar de la muestra yn 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 (σ) de la muestra, mientras fijamos la media (μ) 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 5 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 (x¯ – μ0) / (s / n −− √) 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.

\[\color{blue}{Intermezzo~técnico}\] Para poner en práctica estas ideas con R, necesitamos cargar un par de bibliotecas:

  1. 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.
  2. distr6 para crear y manipular distribuciones de probabilidad.
library(furrr)
## Loading required package: future
library(distr6)
## 
## Attaching package: 'distr6'
## The following object is masked from 'package:stats':
## 
##     qqplot
## The following object is masked from 'package:base':
## 
##     truncate
for(lib in c("furrr", "distr6")) {
  library(lib, character.only = TRUE)
}

Para activar el cálculo paralelo, necesitamos el siguiente código:

plan(multiprocess)
## Warning: Strategy 'multiprocess' is deprecated in future (>= 1.20.0). Instead,
## explicitly specify either 'multisession' or 'multicore'. In the current R
## session, 'multiprocess' equals 'multisession'.

Finalmente, para asegurarnos de que siempre obtengamos la misma salida, necesitamos inicializar el generador de números aleatorios:

set.seed(2019)

\[\color{blue}{Estimación~de~parámetros~de~la~población}\] Existen métodos de diferencias 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:

sigma_hat = sd(sample)
sigma_hat
## [1] 1.986663

El otro parámetro de la distribución Normal lo fija la hipótesis nula:

mu = 0

\[\color{blue}{Construyendo~la~distribución~muestral}\] Primero, construyamos la distribución normal que describe la población:

population = Normal$new(mean = mu, sd = sigma_hat)
population
## Norm(mean = 0, sd = 1.98666303131658)

Luego, necesitamos una función que calcule el estadístico (x¯ – μ0) / (s / n −− √) 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
sampling_distribution = future_map_dbl(1:N, ~ calc_statistic(population$rand(5)))
## Deprecated. Use $properties$kurtosis instead.
## Deprecated. Use $properties$skewness instead.
## Deprecated. Use $properties$support instead.
## Deprecated. Use $properties$symmetry instead.
## Deprecated. Use $traits$type instead.
## Deprecated. Use $traits$valueSupport instead.
## Deprecated. Use $traits$variateForm instead.
## 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".
## 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".
head(sampling_distribution)
## [1]  0.4046007 -0.5544431 -0.7818772 -0.5914122  0.2107333  1.0442065

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

plot(density(sampling_distribution,), xlab = "t statistic", main = "", xlim = c(-10,10)) # Histogram of sampling_distribution
curve(dt(x,4), -10, 10, add = T, col = 2, n = 1e3) # Theoretical distribution t(n-1)
abline(v = calc_statistic(sample), col = 3) # Observed statistic

Podemos ver que, de hecho, la distribución muestral es la distribución t de Student con 4 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.

\[\color{blue}{Calcular~el~p-valor}\] 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 sampling distribution es tan simple como calcular la fracción de casos en los que se cumple esta condición:

sum(sampling_distribution >= calc_statistic(sample))/N
## [1] 0.03014

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

1 - pt(calc_statistic(sample), 4)
## [1] 0.02946129

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(sample, alternative = "greater")$p.val
## [1] 0.02946129

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 de alrededor del 3,7% 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, pero los utilizaré para publicaciones futuras.

\[\color{blue}{¿Por~qué~no~utilizar~la~media~muestral~como~estadística?}\] Usar la media muestral (x¯) como estadística en lugar de (x¯ – μ0) / (s / n −− √) 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
sampling_distribution = future_map_dbl(1:N, ~ mean(population$rand(5)))
## Deprecated. Use $properties$kurtosis instead.
## Deprecated. Use $properties$skewness instead.
## Deprecated. Use $properties$support instead.
## Deprecated. Use $properties$symmetry instead.
## Deprecated. Use $traits$type instead.
## Deprecated. Use $traits$valueSupport instead.
## Deprecated. Use $traits$variateForm instead.
## 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".
## 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".
head(sampling_distribution)
## [1]  0.3465767 -0.5515707 -1.2287648  0.4815007 -1.2585042  0.2225078

que se parece a:

plot(density(sampling_distribution, bw = "SJ"), main = "", xlab = "sample mean") # Histogram of sampling_distribution
abline(v = mean(sample), col = 2) # Observed statistic

El valor p se calcula de la misma manera que antes:

sum(sampling_distribution >= mean(sample))/N
## [1] 0.00438

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 5.