Se va a revisar qué tan bien el teorema del límite central se aproxima a los promedios de muestra de los datos que se tienen. Se aprovecha todo el conjunto de datos de población para comparar los resultados que obtenemos al muestrear realmente la distribución con lo que predice el Teorema del límite central.

Para los siguientes ejercicios se usará el siguiente conjunto de datos.

dir <- system.file(package = "dagdata")
filename <- file.path(dir,"extdata/mice_pheno.csv") 
datos <- read.csv(filename) %>% na.omit # Ya se había leído previamente
head(datos)
##   Sex Diet Bodyweight
## 1   F   hf      31.94
## 2   F   hf      32.48
## 3   F   hf      22.82
## 4   F   hf      19.92
## 5   F   hf      32.22
## 6   F   hf      27.50

Se van a seleccionar solo las hembras, dado que las hembras y los machos tienen diferentes pesos. Se realizará de dos maneras diferentes para ver las opciones.

# Primero se hara uso de filter y select de la librería dplyr
hembrasControl <- datos %>% filter(Sex == "F" & Diet == "chow") %>% select(Bodyweight) %>% unlist
hembrasTratamiento <- datos$Bodyweight[datos$Sex == "F" & datos$Diet == "hf"]

Con los datos obtenidos se pueden calcular los parametros de población usando las funciones mean y popvar.

mu_hTratamiento <- mean(hembrasTratamiento)
mu_hControl <- mean(hembrasControl)

var_hTratamiento <- mean((hembrasTratamiento - mean(hembrasTratamiento))^2) # [1] 25.70358
popvar(hembrasTratamiento) # [1] 25.70358 # Función popvar de la librería rafalib
identical(var(hembrasTratamiento), var_hTratamiento) # [1] FALSE
identical(popvar(hembrasTratamiento), var_hTratamiento) # [1] TRUE
var_hControl <- popvar(hembrasControl)

Se muestra la diferencia entre las dos medias calculadas.

print(mu_hTratamiento - mu_hControl) # [1] 2.375517

Como se puede ver se usa popvar y popsd para calcular la varianza y la distribución estándar de la población, dado que las funciones var y sd no calculan la varianza \(\sigma^2\) y desviación estándar \(\sigma\) de la población, sino que calculan la varianza \(s^2\) y desviación estándar \(s\) de una muestra.

Entonces, para ser matemáticamente correctos, no usamos sd ovar. En su lugar, usamos las funciones popvar ypopsd en rafalib:

library(rafalib)
sd_hf <- popsd(hembrasTratamiento)
sd_control <- popsd(hembrasControl)

En la práctica no es común que se calculen estos parámetros de población, son valores que no se ven de manera directa. En general, son estimados a partir de muestras.

N <- 12
hf <- sample(hembrasTratamiento, 12)
control <- sample(hembrasControl, 12)

Como se dijo, el Teorema del límite central dice que para \(N\) grandes, el conjunto de datos sigue aproximadamente una distribución normal con la media poblacional como promedio y con un error estandar igual a la desviación estándar poblacional dividida por \(\sqrt{N}\). Se mencionó que una regla general es que \(N\) debe mayor que 30. Sin embargo, esa es solo una regla empírica ya que la precisión de la aproximación depende de la distribución de la población. Aquí se puede verificar la aproximación y se hace para varios valores de \(N\).

Ahora se usará sapply yreplicate en lugar de bucles for, lo que hace que el código sea más limpio (no hay que preasignar un vector, R se encarga de esto):

Ns <- c(3,12,25,50)
B <- 10000 # número de simulaciones
res <-  sapply(Ns, function(n) {
  replicate(B,mean(sample(hembrasTratamiento,n))-mean(sample(hembrasControl,n)))
})

Se puede usar qq-plots para ver qué tan bien funcionan las aproximaciones TCE para los datos simulados. Si la distribución normal es una buena aproximación, los puntos deberían caer en línea recta en comparación con los cuantiles normales. Cuanto más se desvía, peor es la aproximación. En el título se pone el promedio y la DE de la distribución observada, se ve cómo la DE disminuye con \(\sqrt{N}\) como se predijo.

mypar(2,2)
for (i in seq(along=Ns)) {
  titleavg <- signif(mean(res[,i]),3)
  titlesd <- signif(popsd(res[,i]),3)
  title <- paste0("N=",Ns[i],", media=",titleavg,", SD=",titlesd)
  qqnorm(res[,i],main=title)
  qqline(res[,i],col=2)
}
Gráfica de cuantiles teóricos contra cuantiles de las muestras.

Gráfica de cuantiles teóricos contra cuantiles de las muestras.

Se ve un ajuste bastante bueno incluso para \(N=3\). ¿Por qué? Debido a que la población en sí está relativamente cerca de la distribución normal, los promedios de las muestras también están cerca de lo normal (la suma de las normales también es normal). En la práctica, en realidad se calcula una relación: se divide por la desviación estándar estimada. Aquí es donde el tamaño de la muestra comienza a importar más.

Ns <- c(3,12,25,50) # tamaño de las muestras
B <- 10000 # número de simulacionoes
##  función para calcular una estadistica t
computetstat <- function(n) {
  y <- sample(hembrasTratamiento,n) # muestras de hembras de tratamiento de tamaños 3, 12, 25 y 50
  x <- sample(hembrasControl,n) # muestras de hembras de control de tamaños 3, 12, 25 y 50
  (mean(y)-mean(x))/sqrt(var(y)/n+var(x)/n) # Calculo de t
}
res <-  sapply(Ns,function(n) {
  replicate(B,computetstat(n))
})
mypar(2,2)
for (i in seq(along=Ns)) {
  qqnorm(res[,i],main=Ns[i])
  qqline(res[,i],col=2)
}
Grafica de cuantiles teóricos contra cuantiles de cuatro muestras de diferentes tamaños.

Grafica de cuantiles teóricos contra cuantiles de cuatro muestras de diferentes tamaños.

Se puede ver que para $ N = 3 $, el TLC no proporciona una aproximación que se pueda usar Para \(N=12\), hay una ligera desviación en los valores más altos, aunque la aproximación parece útil. Para 25 y 50, la aproximación es más acertada.

Hay que tener muy presente que esta simulación sólo demuestra que \(N=12\) es lo suficientemente grande en este caso particular, no en general. Como se mencionó anteriormente, no se puede realizar esta simulación en la mayoría de las situaciones. Solo se usa la simulación en esta ocasión para ilustrar los conceptos detrás del TLC y sus limitaciones. En secciones futuras, se describirán los enfoques que realmente se usan en la práctica.

Ejercicios de TLC y distribución t Capítulo de inferencia Ejercicios TLC en la práctica