La versión original en inglés puede ser consultada en este enlace

Se van a explorar variables aleatorias más a fondo. Imagine que en realidad se tiene el peso de todas las hembras de control y se pueden cargar en R. En estadística se refiere a este caso como la población. Estos son todos los ratones de control disponibles de los que se seleccionan 24. En la práctica no es usual que se se tenga acceso a la población, en este caso se tiene un conjunto de datos especial que se usa para ilustrar conceptos. Se cargan los datos de la población:

dir <- system.file(package = "dagdata")
filename <- file.path(dir,"extdata/femaleControlsPopulation.csv") 
poblacion <- read.csv(filename) %>% unlist

Luego se toma una muestra de 12 ratones tres veces y se mira cómo cambia el promedio de las muestras.

set.seed(1) # Se coloca esta semilla para que se obtenga el mismo resultado
muestra <- sample(poblacion,12)
mean(muestra) # [1] 23.47667
muestra <- sample(poblacion,12)
mean(muestra) # [1] 25.12583
muestra <- sample(poblacion,12)
mean(muestra) # [1] 23.72333

Los promedios varían, si se sigue haciendo esto repetidamente se puede entender el concepto de variable aleatoria. ### Hipótesis nula

Ahora se vuelve a la diferencia promedio de difobservaciones. Como científicos, se debe ser escéptico. ¿Cómo se sabe que este difobservaciones se debe a la dieta? ¿Qué sucede si se da a los 24 ratones la misma dieta? ¿Se verá una diferencia igual de grande? Los estadísticos se refieren a este escenario como la hipótesis nula. El nombre “nulo” se usa para recordar que se esta actuando como escépticos: se da crédito a la posibilidad de que no haya diferencia.

Como se tiene acceso a la población, se puede observar tantos valores como se quiera de la diferencia de los promedios cuando la dieta no tiene ningún efecto. Se puede hacer esto muestreando aleatoriamente 24 ratones de control, dándo la misma dieta y luego registrando la diferencia en la media entre dos grupos divididos aleatoriamente de 12 y 12. Aquí está este proceso escrito en código R:

set.seed(1)
control <- sample(poblacion,12) # 12 ratones de muestra
tratamiento <- sample(poblacion,12) # otros 12 ratones de control
print(mean(tratamiento) - mean(control)) # [1] 1.649167

Ahora se hace 10.000 veces. Primero se usa un “bucle-for”, una operación que permite automatizar esto.

set.seed(1)
n <- 10000
nulos <- vector("numeric",n)
for (i in 1:n) {
  control <- sample(poblacion,12)
  tratamiento <- sample(poblacion,12)
  nulos[i] <- mean(tratamiento) - mean(control)
}

Ahora se hace con la funcion replicate, los resultados son los mismos pero usar replicate es más limpio y eficiente

set.seed(1)
n <- 10000
N <- 12
nulos <- replicate(n,{
  control <- sample(poblacion,N)
  tratamiento <- sample(poblacion,N)
  mean(tratamiento) - mean(control)
})

Los valores en “nulo” forman lo que se llama la distribución nula. Se definirá esto de manera más formal a continuación.

Se verifica la proporción de estas diez mil diferencias entre las medias que son mayores que la diferencia entre las observaciones difobservaciones .

mean(nulos >= difobservaciones) # [1] 0.0125

Solo una pequeña proporción de las 10.000 simulaciones. Como escépticos, ¿qué se concluye? Cuando no hay un efecto de la dieta, se ve una diferencia tan grande como la que se observa solo 1,5% de las veces. Esto es lo que se conoce como valor p.

Distribuciones

Se ha explicado lo que se quiere decir con nulo en el contexto de la hipótesis nula, pero ¿qué es exactamente una distribución? La forma más sencilla de pensar en una distribución es como una descripción compacta de muchos números. Por ejemplo, se supone que se ha medido la altura de todos los hombres de una población. Se necesita describir estos números a alguien que no tiene idea de cuáles son estas alturas, como un extraterrestre que nunca ha visitado la Tierra. Se supone que todas estas alturas están contenidas en el siguiente conjunto de datos:

set.seed(1)
data(father.son, package="UsingR")
x <- father.son$fheight

Un método para resumir estos números es simplemente enumerarlos todos para que los vea el extraterrestre. Aquí hay 10 alturas (en pulgadas) seleccionadas al azar de las 1078 que hay en el conjunto de datos:

round(sample(x,10), 1) #  [1] 70.6 67.4 70.1 68.5 62.7 67.1 66.5 67.6 67.2 68.5

Función de distribución acumulada

Esto se denomina función de distribución acumulativa (CDF). Cuando la CDF se deriva de los datos, en lugar de teóricamente, también se llama CDF empírica (ECDF). El ECDF para los datos de altura se ve así y usualmente se conoce como gráfica en ojiva:

\[ F(a)\equiv{Pr(x\leq{a})} \]

plot(ecdf(x), main="", xlab = "Altura de los padres en pulgadas", ylab = TeX("$F(a)\\equiv Pr(x\\leq a)$"))

Histogramas

Aunque el concepto de CDF empírico se discute ampliamente en los libros de texto de estadística, la grafíca de ojiva no es muy popular en la práctica. La razón es que los histogramas dan la misma información y son más fáciles de interpretar. Los histogramas muestran la frecuencia de valores en intervalos:

\[ \mbox{Pr}(a\leq{x}\leq{b}) = F(b)-F(a) \]

Trazar estas alturas como barras es lo que se conoce como histograma. Es una gráfica útil cuando se está interesado en los intervalos, tal o cual porcentaje está entre 70 pulgadas y 71 pulgadas, etc., en lugar del porcentaje menor que una altura en particular.

También es fácil distinguir diferentes tipos (familias) de distribuciones al observar histogramas. Aquí hay un histograma de alturas:

bins <- seq(floor(min(x)), ceiling(max(x)))
hist(x, breaks = bins, xlab =  "Altura de los padres en pulgadas", main = "Altura de los padres")

Mostrar esta gráfica al extraterrestre es mucho más informativo que mostrar números. Con esta gráfica simple, se puede aproximar el número de individuos en cualquier intervalo dado. Por ejemplo, hay alrededor de 70 individuos de más de seis pies (72 pulgadas) de altura.

Gráficos combinados

Se pueden usar también gráficos combinados cuando se quiere mostrar de diferentes maneras y en una sola presentación el comportamiento de los datos, por ejemplo se muestra el histograma, la ojiva y la curva normal.

par(mar = c(5, 4, 4, 4) + 0.3)  # Se deja espacio para el eje z
hist(x, col = 'lightblue', main="Altura de los padres", xlab="Altura en pulgadas", ylab="Frecuencia")
# Se crea un vector de 100 elementos para graficar la ojiva y la campana normal
bins = seq(min(mean(x)-3*sd(x),min(x)), max(mean(x)+3*sd(x),max(x)), len=100)
# Se crea el vector del acumulado de alturas expresado en porcentaje
x.cut = cut(x, bins, right=FALSE) 
x.freq = table(x.cut) 
x.cumfreq = round(100 * cumsum(x.freq) / length(x),2)

# Se grafica la ojiva
par(new = TRUE)
plot(bins, c(0,x.cumfreq), type="l", axes=F, bty="n", xlab="", ylab ="", col="red", lwd=2)
axis(side=4, at = pretty(range(x.cumfreq)))
mtext("%", side=4, line=3)

# Se grafica la curva normal (campana)
par(new=TRUE)
de <- dnorm(bins, mean=mean(x), sd=sd(x))
plot(bins, de, type="l", lwd=2, axes=FALSE, bty = "n", xlab = "", ylab = "",col="black")

Distribuciones de probabilidad

Resumir listas de números es un uso poderoso de la distribución de probabilidad. Un uso aún más importante es describir los posibles resultados de un variable aleatoria. A diferencia de una lista fija de números, en realidad no observamos todos los resultados posibles de las variables aleatorias, por lo que en lugar de describir proporciones, describimos probabilidades. Por ejemplo, si elegimos una altura aleatoria de nuestra lista, entonces la probabilidad de que caiga entre \(a\) y \(b\) se denota con:

\[ \mbox{Pr}(a \leq X \leq b) = F(b) - F(a) \]

\(X\) ahora está en mayúsculas para distinguirlo como una variable aleatoria y que la ecuación anterior define la distribución de probabilidad de la variable aleatoria. Conocer esta distribución es increíblemente útil en ciencia. Por ejemplo, en el caso anterior, si se conoce la distribución de la diferencia en la media de los pesos de los ratones cuando la hipótesis nula es verdadera, conocida como distribución nula, podemos calcular la probabilidad de observar un valor tan grande como se hizo, denominado p-valor. En una sección anterior se ejecutó lo que se llama una simulación de Monte Carlo (se darán más detalles sobre la simulación de Monte Carlo en una sección posterior) y se obtuvieron 10,000 resultados de la variable aleatoria bajo la hipótesis nula. Se repite el ciclo anterior, pero esta vez se agrega un punto a la figura cada vez que se ejecuta el experimento. Si se ejecuta este código, puede ver la distribución nula que se forma a medida que los valores observados se apilan uno encima del otro.

n <- 100
nullplot(-5,5,1,30, xlab="Diferencias observadas (gramos)", ylab="Frecuencia") 
totals <- vector("numeric",11)
set.seed(1)
for (i in 1:n) {
  control <- sample(poblacion,12)
  treatment <- sample(poblacion,12)
  nulldiff <- mean(treatment) - mean(control)
  j <- pmax(pmin(round(nulldiff)+6,11),1)
  totals[j] <- totals[j]+1
  text(j-6,totals[j],pch=15,round(nulldiff,1))
  ##if(i < 15) Sys.sleep(1) ## Agregar esta línea para ver los valores aparecer lentamente
}

La figura anterior equivale a un histograma. A partir de un histograma del vector “nulos” que se calculó anteriormente, se puede ver que valores tan grandes como difobservaciones son relativamente raros:

hist(nulos, freq=TRUE, main="Histograma de nulos")
abline(v=difobservaciones, col="red", lwd=2)
Distribución Nula con la diferencia observada marcada con una línea roja vertical.

Distribución Nula con la diferencia observada marcada con una línea roja vertical.

Un punto importante a tener en cuenta aquí es que mientras se define \(\mbox{Pr}(a)\) contando casos, se aprende que, en algunas circunstancias, las matemáticas dan fórmulas para \(\mbox{Pr}(a)\) que horran la molestia de calcular como se hizo aquí. Un ejemplo de este poderoso enfoque utiliza la aproximación de distribución normal.

Distribución normal

La distribución de probabilidad que vemos arriba se aproxima a una que es muy común en la naturaleza: la curva de campana, también conocida como distribución normal o distribución gaussiana. Cuando el histograma de una lista de números se aproxima a la distribución normal, se puede usar una fórmula matemática para aproximar la proporción de valores o resultados (área bajo la curva) en cualquier intervalo dado: \[ \mbox{Pr}(a<x<b)=\int_{a}^{b} \dfrac{1}{\sqrt{2\pi\sigma^2}}exp\left(\dfrac{-{(x-\mu)}^2}{2\sigma^2}\right)\, dx \]

Si bien la fórmula puede parecer intimidante, no hay que preocuparse, nunca habrá que escribirla, ya que está almacenada en una forma más conveniente (como pnorm en R que establece a en \(-\infty\), y toma b como argumento).

Aquí se hace referencia a \(\mu\) y \(\sigma\) como la media y la desviación estándar de la población, que se calculan así:

\[ \mu=\dfrac{1}{M}\sum_{i=1}^{M}x_{i} \]

\[ \sigma^2=\dfrac{1}{M}\sum_{i=1}^{M}(x_{i}-\mu)^2 \] Estandarización: si se encuentra que los datos se distribuyen de manera normal se pueden transformar a unidades estándar, se resta la media y se divide por la desviación estandar, para un determinado punto \(X_i\) de la distribución. \[ Z_i=\dfrac{X_i-\bar{X}}{\sigma_X} \]

Si esta aproximación normal es válida para esta lista, entonces la media poblacional y la varianza de esta lista se pueden usar en la fórmula anterior. Un ejemplo de esto sería cuando se notó anteriormente que menos del 1.5% de los valores en la distribución nula estaban por encima de difobservaciones. Se puede calcular la proporción de valores por debajo de un valor x con pnorm(x, mu, sigma) sin conocer todos los valores.

La aproximación normal funciona muy bien aquí, por ejemplo para calcular la proporción de resultados nulos que están por encima de la diferencia observada, la línea roja en la gráfica se tiene.

1 - pnorm(difobservaciones,mean(nulos),sd(nulos)) # [1] 0.01247381

Más adelante, se aprenderá que hay una explicación matemática para esto. Una característica muy útil de esta aproximación es que solo se necesita conocer \(\mu\) y \(\sigma\) para describir la distribución completa. A partir de esto, se puede calcular la proporción de valores en cualquier intervalo.

Por ejemplo, para la distribución normal estándar se tiene:

pnorm(1)-pnorm(-1) # [1] 0.6826895
puntosx <- c(-1, seq(-1, 1, length=100), 1)
puntosy <- c(0, dnorm(seq(-1, 1, length=100)), 0)

de <- dnorm(seq(-3, 3, length=600), mean=0, sd=1)
plot(seq(-3, 3, length=600), de, type="l", lwd=2, axes=FALSE, bty = "n", xlab = "", ylab = "",col="black")
polygon(puntosx, puntosy, col='gray')

axis(1, at=seq(-3, 3, by=1), labels = FALSE)
lablist = as.vector(seq(-3, 3,1))
text(x = c(-3, -2, -1, 0, 1, 2, 3), y = 0, labels = c(-3, -2, -1, 0, 1, 2, 3), adj = c(0, 3), xpd=T)

El 68,29% de los datos están a una distancia menor o igual a 1 desviación estándar:

pnorm(2)-pnorm(-2) # [1] 0.9544997
puntosx <- c(-2, seq(-2, 2, length=200), 2)
puntosy <- c(0, dnorm(seq(-2, 2, length=200)), 0)

de <- dnorm(seq(-3, 3, length=600), mean=0, sd=1)
plot(seq(-3, 3, length=600), de, type="l", lwd=2, axes=FALSE, bty = "n", xlab = "", ylab = "",col="black")
polygon(puntosx, puntosy, col='gray')

axis(1, at=seq(-3, 3, by=1), labels = FALSE)
lablist = as.vector(seq(-3, 3,1))
text(x = c(-3, -2, -1, 0, 1, 2, 3), y = 0, labels = c(-3, -2, -1, 0, 1, 2, 3), adj = c(0, 3), xpd=T)

El 95,45% de los datos están a una distancia menor o igual a 2 desviaciones estándar:

pnorm(3)-pnorm(-3) # [1] 0.9973002
puntosx <- c(-3, seq(-3, 3, length=300), 3)
puntosy <- c(0, dnorm(seq(-3, 3, length=300)), 0)

de <- dnorm(seq(-3, 3, length=600), mean=0, sd=1)
plot(seq(-3, 3, length=600), de, type="l", lwd=2, axes=FALSE, bty = "n", xlab = "", ylab = "",col="black")
polygon(puntosx, puntosy, col='gray')

axis(1, at=seq(-3, 3, by=1), labels = FALSE)
lablist = as.vector(seq(-3, 3,1))
text(x = c(-3, -2, -1, 0, 1, 2, 3), y = 0, labels = c(-3, -2, -1, 0, 1, 2, 3), adj = c(0, 3), xpd=T)

El 99,73% de los datos están a una distancia menor o igual a 3 desviaciones estándar:

pnorm(6)-pnorm(-6) # [1] 1

El 99,9999% de los datos están a una distancia menor o igual a 6 desviaciones estándar:

Resumen

Entonces, calcular un valor p para la diferencia en la dieta de los ratones fue bastante fácil. Pero, ¿por qué no se ha terminado? Para hacer el cálculo, se hizo el equivalente a comprar todos los ratones disponibles en The Jackson Laboratory y realizar el experimento repetidamente para definir la distribución nula. Sin embargo, esto no es algo que se pueda hacer en la práctica. La inferencia estadística es la teoría matemática que permite aproximar esto solo con los datos de una muestra, es decir, los 24 ratones originales.

Introducción capítulo Capítulo de inferencia Ejercicios variables aleatorias