En este ejercicio nos aproximaremos al supuesto de normalidad del modelo de regresión lineal. Antes que nada debemos de recordar que el modelo pretende modelar la media de la población; esta es la clave no se trata de la modelación de la población sino de la media de la población. Esto tomando en cuenta el tamaño de la muestra, y la idea intuitiva de que cuando la Muestra es grande, los datos se comportan de manera normal. Veamos pues una curva de una normal:
hist(rnorm(5000, 1, 1), freq = 0, main = "y~N(1,1)", ylim = c(0, 0.45))
lines(seq(-4, 4, by = 0.01), dnorm(seq(-4, 4, by = 0.01), 1, 1), col = 2)
Ahora utilizaremos datos provenientes de otra función, digamos una gamma, con parámetros 1,1
hist(rgamma(5000, 1, 1), freq = 0, main = expression(y %~% Gamma(1, 1)), ylim = c(0,
1), xlim = c(0, 6))
lines(seq(0, 5, by = 0.001), dgamma(seq(0, 5, by = 0.001), 1, 1), col = 2)
Es claro que la distribución no es ni se parece a una normal, ahora supongamos que a partir de esta distribución sacamos cuatro muestras de tamaños diferentes digamos 10,50,100 y 1000;si todo falla pues digamos que una muestra de mil es los suficientemente grande:
par(mfrow = c(2, 2), mar = c(2, 2, 2, 2))
hist(rgamma(10, 1, 1), freq = 0, main = "n=10", xlab = "", ylab = "")
hist(rgamma(50, 1, 1), freq = 0, main = "n=50", xlab = "", ylab = "")
hist(rgamma(100, 1, 1), freq = 0, main = "n=100", xlab = "", ylab = "")
hist(rgamma(1000, 1, 1), freq = 0, main = "n=1000", xlab = "", ylab = "")
Es evidente que entre más grande sea la muestra el histograma se aproxima más a la distribución original, y no a una normal. Esta es una interpretación erronea, sin embargo y regresando a la idea del modelo de regresión lineal, en realidad lo que nos interesa es la media, tomemos entonces en lugar de una muestra \( n \) muestras, digamos de tamaño pequeño \( m=10 \) es decir que tomaremos \( n={10,50,10,1000} \) muestras de tamaño 10. Y veamos que sucede:
n10 <- mean(rgamma(10, 1, 1))
while (length(n10) < 10) {
n10 <- c(n10, mean(rgamma(10, 1, 1)))
}
n50 <- mean(rgamma(10, 1, 1))
while (length(n50) < 50) {
n50 <- c(n50, mean(rgamma(10, 1, 1)))
}
n100 <- mean(rgamma(10, 1, 1))
while (length(n100) < 100) {
n100 <- c(n100, mean(rgamma(10, 1, 1)))
}
n1000 <- mean(rgamma(10, 1, 1))
while (length(n1000) < 1000) {
n1000 <- c(n1000, mean(rgamma(10, 1, 1)))
}
par(mfrow = c(2, 2), mar = c(2, 2, 5, 2))
hist(n10, freq = 0, xlab = "", ylab = "", main = "n10")
hist(n50, freq = 0, xlab = "", ylab = "", main = "n50")
hist(n100, freq = 0, xlab = "", ylab = "", main = "n100")
hist(n1000, freq = 0, xlab = "", ylab = "", main = "n1000")
Notemos que en este caso conforme n aumenta la media tiende a distribuirse de manera normal, aun cuando la distribución de la variable de interes no es normal. Esto es el teorema central del límite.
Recordemos que dados los parámetros de la distribución la \( E[x]=\frac{k}{\lambda}=\frac{1} {1}=1 \) veamos esto dado nuestros tamaños de muestra, por otra parte podemos aplicar pruebas de normalidad en este caso aplicaremos la Kolmorogov-Smirnoff y la de Shapiro-Wilks, en la tabla se muestran los p-valores de las pruebas:
| n | \( \bar{x} \) | p - KS | p - SW |
|---|---|---|---|
| 10 | 1.0012 | 5.1699 × 10-6 | 0.1607 |
| 50 | 1.0308 | 8.8818 × 10-16 | 0.7471 |
| 100 | 0.9667 | 0 | 0.0174 |
| 1000 | 1.009 | 0 | 1.7436 × 10-9 |
Como no fijamos la generación de números aleatorios, cada que se ejecute el código podemos esperar resultados diferentes cada vez que se ejecute el script, sin enmbargo esto se permitió ya que la tendencia debe de mantenerse en todos los casos, i.e. que entre mayor sea el tamaño de muestra la distribución de la media claramente se aproxima más a la normal, y el parámetro del valor esperado también.