Introducción

Ahora se demostrará cómo obtener un valor p en la práctica. Se comienza cargando datos experimentales y luego hay unos pasos utilizados para formar una estadística t y calcular un valor p. Se puede realizar esta tarea con solo unas pocas líneas de código (en al final de esta sección se pueden ver). Sin embargo, para comprender los conceptos, se construirá una estadística t desde “cero”.

Leer y preparar datos

Se comienza leyendo los datos.

dir <- system.file(package = "dagdata")
filename <- file.path(dir,"extdata/femaleMiceWeights.csv") 
tratamiento <- read.csv(filename)

Un primer paso importante es identificar qué filas están asociadas con el tratamiento y el control, y calcular la diferencia de las medias.

control <- filter(tratamiento,Diet=="chow") %>% select(Bodyweight) %>% unlist
tratamiento <- filter(tratamiento,Diet=="hf") %>% select(Bodyweight) %>% unlist
diff <- mean(tratamiento) - mean(control)
print(diff) # [1] 3.020833

Se solicita que se informe un valor p, ¿qué se debe hacer? Se aprendio que la diferencia de medias diff, denominado tamaño del efecto observado, es una variable aleatoria. También se aprendió que bajo la hipótesis nula, la media de la distribución de “diff” debería ser 0. ¿Qué pasa con el error estándar? También se aprendió que el error estándar de esta variable aleatoria es la desviación estándar de la población dividida por la raíz cuadrada del tamaño de la muestra:

\[ SE(\bar{X}) = \dfrac{\sigma_x}{\sqrt{N}}\]

Se usa la desviación estándar de la muestra como una estimación de la desviación estándar de la población. En R, simplemente se usa la función sd y el error estándar SE es:

sd(control)/sqrt(length(control)) # [1] 0.8725323
## [1] 0.8725323

Este es el error estándar del promedio de la muestra, pero en realidad se quiere el error estándar de la diferencia de las medias muestrales “diff”. Se vio cómo la teoría estadística dice que la varianza de la diferencia de dos variables aleatorias es la suma de sus varianzas, entonces se calcula la varianza y se saca la raíz cuadrada:

se <- sqrt( 
  var(tratamiento)/length(tratamiento) + var(control)/length(control) 
) # [1] 1.469867

La teoría estadística dice que si dividimos una variable aleatoria por su error estándar, se obtiene una nueva variable aleatoria con un error estándar de 1.

tstat <- diff/se

Esta relación es lo que se llama estadístico t. Es la razón de dos variables aleatorias y, por lo tanto, una variable aleatoria. Una vez que se conoce la distribución de esta variable aleatoria, se puede calcular fácilmente un valor p.

Como se explicó en la sección anterior, el teorema del límite central dice que para tamaños de muestra grandes, tanto los promedios de la muestra mean(tratamiento) como la mean(control) son normales. La teoría estadística dice que la diferencia de dos variables aleatorias distribuidas normalmente es también normal, por lo que el teorema del límite central dice que tstat es aproximadamente normal con media 0 (la hipótesis nula) y desviación estándar 1 (dividimos por su SE).

Entonces, ahora para calcular un valor p, todo lo que hay que hacer es preguntar: ¿con qué frecuencia una variable aleatoria distribuida normalmente excede diff? R tiene una función incorporada, pnorm, para responder esta pregunta específica. pnorm(a) devuelve la probabilidad de que una variable aleatoria que sigue la distribución normal estándar caiga por debajo de a. Para obtener la probabilidad de que sea mayor que “a”, simplemente se usa 1-pnorm(a). Si se quiere saber la probabilidad de ver algo tan extremo como diff: más pequeño (más negativo) que -abs(diff) o más grande que abs(diff), se llama a estas dos regiones “colas” y se calcula su tamaño:

righttail <- 1 - pnorm(abs(tstat)) 
lefttail <- pnorm(-abs(tstat))
pval <- lefttail + righttail
print(pval) # [1] 0.0398622

En este caso, el valor p es menor que 0.05 y si se usa el corte convencional de 0.05, se llamaría a la diferencia estadísticamente significativa.

Ahora hay un problema. El teorema del límite central funciona para muestras grandes, pero ¿12 es lo suficientemente grande? Una regla general para el teorema del límite central es que 30 es un tamaño de muestra suficientemente grande (pero esto es solo una regla heurística). El valor p que calculamos es solo una aproximación válida si los supuestos se cumplen, lo que no parece ser el caso aquí. Sin embargo, existe otra opción además de usar el teorema del límite central.

La distribución t en la práctica

Como se dijo anteriormente, la teoría estadística ofrece otro resultado útil. Si la distribución de la población es normal, entonces podemos calcular la distribución exacta del estadístico t sin necesidad del teorema del límite central. Este es un gran “si” dado que, con muestras pequeñas, es difícil verificar si la población es normal. Pero para algo como el peso, se puede suponer que la distribución de la población probablemente se aproxima bien a la normal y que se puede usar esta aproximación. Además, se puede mirar un gráfico qq para las muestras. Esto muestra que la aproximación es al menos cercana:

library(rafalib)
mypar(1,2)
qqnorm(tratamiento, main = 'Normal Q-Q Plot tratamiento')
qqline(tratamiento,col=2)
qqnorm(control, main = 'Normal Q-Q Plot control')
qqline(control,col=2)
Gráfica de Quantile-quantile para cada muestra contra la distribución normal teórica.

Gráfica de Quantile-quantile para cada muestra contra la distribución normal teórica.

Si se usa esta aproximación, entonces la teoría estadística dice que la distribución de la variable aleatoria tstat sigue una distribución t. Esta es una distribución mucho más complicada que la normal. La distribución t tiene un parámetro de ubicación como el normal y otro parámetro llamado grados de libertad. R tiene una buena función que calcula todo

t.test(tratamiento, control)
## 
##  Welch Two Sample t-test
## 
## data:  tratamiento and control
## t = 2.0552, df = 20.236, p-value = 0.053
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.04296563  6.08463229
## sample estimates:
## mean of x mean of y 
##  26.83417  23.81333

Si se quiere ver solo el valor p entonces se puede usar el extractor $:

result <- t.test(tratamiento,control)
result$p.value
## [1] 0.05299888

El valor p es un poco mayor ahora. Esto es de esperar porque la aproximación de teorema del límite central consideró el denominador de tstat prácticamente fijo (con muestras grandes prácticamente lo es), mientras que la aproximación de distribución t tiene en cuenta que el denominador (el error estándar de la diferencia) es un valor aleatorio. variable. Cuanto menor sea el tamaño de la muestra, más varía el denominador.

Puede resultar confuso que una aproximación haya dado un valor p y otra haya dado otro, porque se espera que haya solo una respuesta. Sin embargo, esto no es infrecuente en el análisis de datos. Se usaron diferentes supuestos, diferentes aproximaciones y, por lo tanto, se obtuvieron diferentes resultados.

Más adelante, en la sección de cálculo de potencia, se describirán los errores de tipo I y tipo II. Como vista previa, se señala que la prueba basada en la aproximación de teorema del límite central tiene más probabilidades de rechazar incorrectamente la hipótesis nula (un falso positivo), mientras que la distribución t tiene más probabilidades de aceptar incorrectamente la hipótesis nula (falso negativo).

Ejecución de la prueba t en la práctica

Ahora que se han repasado los conceptos, se puede mostrar el código relativamente simple que se usaría para calcular una prueba t:

library(dplyr)
dir <- system.file(package = "dagdata")
filename <- file.path(dir,"extdata/mice_pheno.csv") 
datos <- read.csv(filename) %>% na.omit
datos <- read.csv("mice_pheno.csv")
control <- datos %>% filter(Diet=="chow") %>% select(Bodyweight) # Filtrado con dplyr
tratamiento <- datos$Bodyweight[datos$Diet=="hf"] # Filtrado usando el selector $
t.test(tratamiento,control)
## 
##  Welch Two Sample t-test
## 
## data:  tratamiento and control
## t = 7.1932, df = 735.02, p-value = 1.563e-12
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  2.231533 3.906857
## sample estimates:
## mean of x mean of y 
##  30.48201  27.41281

Los argumentos de t.test pueden ser de tipo * data.frame * y, por lo tanto, no es necesario hacer unlist para convertirlos en objetos numéricos.

Ejercicios TLC en la práctica Capítulo de inferencia Intervalos de confianza