En el análisis estadístico de una muestra, una de las cantidades preferidas para representar la variabilidad de nuestros datos es la varianza muestral \(s^2\). La fórmula utilizada es la siguiente: \[s^2 = \frac{\sum_i^n (x_i - \bar{x})^2}{n-1},\] donde \(x_i\) es el valor de la variable \(X\) en la unidad \(i\) de la muestra, \(\bar{x}\) es la media muestral de \(X\) y \(n\) es el tamaño de la muestra.

El denominador puede resultar un poco desconcertante. ¿Por qué se divide entre \(n -1\)? Parece más natural, y más consistente con la interpretación que hacemos de la varianza, dividir entre \(n\); esto es, \[s^2_n = \frac{\sum_i^n (x_i - \bar{x})^2}{n}.\]

En este documento se da respuesta a dicha pregunta recurriendo a técnicas de simulación. El punto central a tomar en cuenta es el objetivo de tomar una muestra aleatoria. La teoría subyacente en el análisis estadístico de una muestra asume que utilizamos los datos observados en esta para conocer características de una población no observada. Este objetivo está implícito cuando nos referimos a las fórmulas utilizadas para calcular la media o la proporción en una muestra con el término estimadores.

Entonces, la fórmula utilizada para calcular la varianza muestral de la variable \(X\) debe juzgarse por sus propiedades como estimador de la varianza de \(X\) en la población de la que proviene la muestra. Veremos que, si bien \(s_n^2\) (dividiendo entre \(n\)) hace un excelente trabajo describiendo la dispersión de los datos en la muestra, como estimador carece de una propiedad deseable que sí está presente en \(s^2\). Mientras que \(s^2\) es un estimador insesgado, \(s^2_n\) no lo es.

¿Qué significa esto? Digamos que se quiere conocer un parámetro poblacional \(\theta\), y para eso se piensa tomar una muestra aleatoria de la población y aplicar el estimador \(\hat{\theta}\). Un estimador es insesgado cuando su valor esperado es igual al parámetro; esto es: \[\mathrm{E}(\hat{\theta}) = \theta.\]

La propiedad de ausencia de sesgo es explotada en diversos procedimientos inferenciales. Por ejemplo, los estimadores de la media y de la proporción son insesgados. Este hecho, junto con otras características del estimador reveladas por el Teorema Central del Límite (TCL), es usado en el procedimiento para para obtener intervalos de confianza a partir de las estimaciones obtenidas en una muestra.

La comparación de los dos estimadores de la varianza nos puede ayudar a comprender mejor la idea de sesgo estadístico.

1 La varianza muestral como estadístico de dispersión

Comenzamos con una ilustración simple de la varianza muestral.

Digamos que medimos la variable \(Y\) en una muestra de cinco observaciones, y que los valores obtenidos son 1, 2, 3, 4 y 5. Entonces, la media de \(Y\) en la muestra es \(\bar{y} = (1 + 2 + 3 + 4 + 5)/5 = 3\). Ahora, para informarnos sobre la dispersión, calculamos la varianza muestral. El primer para es tomar la diferencia del \(Y\) con respecto a la media y elevar esa diferencia al cuadrado: \((y_i - 3)^2\).

Veamos que arroja cada observación.

Y <- c(1:5) # vector con  enteros de 1 o 5.
v.media <- mean(Y) # media muestral
v.media
## [1] 3
dif.media <- v.media - Y # Diferencia con respecto a la media.
dif.media.2 <- dif.media^2 # Cuadrado de la diferencia
t <- data.frame(Y, dif.media, dif.media.2)
t
##   Y dif.media dif.media.2
## 1 1         2           4
## 2 2         1           1
## 3 3         0           0
## 4 4        -1           1
## 5 5        -2           4

Una alternativa para calcular la varianza es \(s_n^2\); es decir, sumamos los cuadrados de la cuarta columna y dividimos el resultado entre el tamaño de la muestra, que es 5.

## Suma de cuadrados dividida entre n (s^2_n)
sum(t$dif.media.2)/5
## [1] 2

La interpretación del resultado es bastante directa: en promedio, el cuadrado de la distancia de las observaciones con respecto a la media es de 2 unidades. Más aún si se toma la raíz cuadrada de la varianza obtenemos la desviación estándar, y podemos decir que en promedio la diferencia entre cada observación y la media muestral es de \(\sqrt{2}\).1

Ahora, si obtenemos la varianza muestral de \(Y\) con la función apropiada en R, var(), obtenemos una cantidad distinta.

## Función de R para obtener varianza muestral
var(t$Y)
## [1] 2.5

Esto es porque var() utiliza \(s^2\); es decir, no divide entre \(n\) sino entre \(n-1\):

## Suma de cuadrados dividida entre n (s^2_n)
sum(t$dif.media.2)/4
## [1] 2.5

Ahora, 2.5 no es realmente el promedio del cuadrado de las diferencias con respecto a la media. De hecho, es mayor. Si mi interés es la dispersión muestral de \(Y\), \(s_n^2\) hace un mejor trabajo que \(s^2\).

Pero en realidad ese no es mi objetivo.

2 Estimador insesgado: la media

Digamos que se desea conocer las características que tiene una variable \(X\) en una población determinada. Para ello, tomamos de la población una muestra aleatoria de tamaño \(n = 30\) y en la muestra se calcula la media y la varianza de \(X\).

Vamos a suponer en la población \(X\) tiene una distribución normal con media \(\mu = 0\) y varianza \(\sigma^2 = 1\); esto es \(X \sim N(0,1)\)2. Ver la figura 2.1. A las características poblacionales como \(\mu\) y \(\sigma^2\) se les llama parámetros.

Distribución poblacional de la variable X ~ N(0,1)

Figure 2.1: Distribución poblacional de la variable X ~ N(0,1)

Vamos a tomar una muestra aleatorio de 30 observaciones de esta población con la función rnorm()3 y observemos los valores de \(X\) obtenidos.

## Tomar una muestra de n = 30 de X ~ N(0,1)
muestra <- rnorm(30)
muestra # valores en la muestra
##  [1]  0.19834463  0.66260061 -0.50851016  0.17244492 -1.32946300  0.69496585
##  [7] -1.81344249 -0.11094755  1.20840026  0.77206142  0.88858300  1.05775439
## [13] -1.54208811 -1.34931209 -1.11109491 -0.29455671  1.06689745  0.28615393
## [19]  0.88246466  0.05498059 -0.73921222 -0.03669539  0.02791200 -0.55752255
## [25]  0.55838756  0.15806131  0.89309819 -0.07353653 -1.74486883  0.83848808

Ahora tomamos la media y la varianza muestrales:

mean(muestra) # media muestral
## [1] -0.02632172
var(muestra)  # varianza muestral
## [1] 0.8114038

Se puede ver que los valores muestrales rondan los valores verdaderos de los respectivos parámetros, pero no son exactamente iguales. Esto es lo que cabe esperar de la naturaleza aleatoria del muestreo.

Para ver esto más claro, tomemos 10 muestras aleatorias y en cada una de ellas vamos a obtener la media.

Para ello, primero vamos a generar una función que en un sólo paso tome una muestra de tamaño \(n\) y obtenga la media muestral.

# Función para generar muestra de tamaño n y obtener media
mi.fun.media <- function(n){ 
  muestra <- rnorm(n)
  mean(muestra)
}

mi.fun.media(30) # media en muestra de n = 30
## [1] -0.3947447

Ahora, con la función replicate() vamos a repetir el procedimiento anterior 10 veces: tomar 10 muestras aleatorias de \(n = 30\) y en cada una calcular la media.

replicate(10, mi.fun.media(30))
##  [1] -0.06633709 -0.48680387  0.01222077 -0.04580272  0.02634239 -0.08121338
##  [7]  0.01554403 -0.04087080 -0.11009388 -0.13555126

Una vez más, nótese que las 10 medias rondan el valor de la media poblacional, pero no son iguales. Sabemos que el estimador de la media es insesgado: \(\mathrm{E}(\bar{x}) = \mu.\) ¿Qué significa esto? Ya vimos que al tomar la media de una variable en una muestra aleatoria podríamos obtener valores distintos. Pensemos entonces, en el conjunto de todos los valores que podríamos obtener al tomar un muestra aleatoria de tamaño \(n = 30\) y calcular la media. Ese conjunto, con su distribución de probabilidad, es la distribución muestral. Podemos hacernos idea muy precisa de sus características si simulamos este proceso un número de veces considerablemente alto. A continuación, se toman 100,000 muestras de tamaño \(n = 30\) y en cada una de ellas se calcula el promedio. Los resultados se guardan en el vector x.b.

x.b <- replicate(100000, mi.fun.media(30))

Para saber si es un estimador insesgado, tomamos la media de las cien mil estimaciones.4

round(mean(x.b), digits = 3)
## [1] 0

Al tomar una muestra aleatoria de una población y calcular el promedio de la variable \(X\), es posible obtener distintos valores. Sin embargo, la media de esos valores posibles es el verdadero valor del parámetro, que es la media de \(X\) en la población. Formalmente, escribimos \[\mathrm{E}(\bar{x}) = \mu.\]

Para visualizar la distribución muestral de manera global, tomamos un histograma de los 100 mil promedios simulados. La línea vertical azul destaca la ubicación de la media de las simulaciones.

Distribución muestral del estimador de la media

Figure 2.2: Distribución muestral del estimador de la media

3 La varianza

Volvamos a las dos fórmulas para la varianza muestra \[s^2_n = \frac{\sum_{i=1}^n(x_i - \bar{x})^2}{n} \qquad \text{y} \qquad s^2 = \frac{\sum_{i=1}^n(x_i - \bar{x})^2}{n-1}.\]

Como decíamos, \(s^2_n\) hace un buen trabajo describiendo la dispersión de la muestra. Sin embargo, vistos como estimadores, resulta que \(s_2\) es insesgado, mientras que \(s_n^2\) es un estimador con sesgo.

Para ver esto, vamos a simular la distribución muestral de ambos estimadores siguiendo un procedimiento igual al de la sección anterior.

Para ver esto, vamos a seguir un procedimiento similar al anterior. Primero vamos a tomar 100,000 mil muestras aleatorias de una población en la que la distribución de \(X\) es normal, con media \(\mu = 0\) y \(\sigma^2 = 1\). En cada una de ellas, vamos a calcular la varianza con \(s^2 = \sum_{i=1}^n(x_i - \bar{x})^2/(n-1)\).5

## Función para muestra de n = 30 y tomar la varianza
mi.fun.var1 <- function(n){ 
  muestra <- rnorm(n)
  sum((muestra - mean(muestra))^2)/(n-1) # Estimador s^2
}
s2 <- replicate(100000, mi.fun.var1(30))

Ahora tomamos el promedio de las 100 mil varianzas estimadas- Recuérdese que, si el estimador es insesgado, este promedio tendría que ser igual a 1, y este es el caso.

round(mean(s2), digits = 3)
## [1] 1.001

Vamos ahora a ubicar el valor esperado del estimador \(s^2\) en la distribución muestral (figura 3.1.

Distribución muestral de $s^2$

Figure 3.1: Distribución muestral de \(s^2\)

Como una nota aparte, recuérdese que el estimador de la varianza no está contemplado por el (TCL), por lo que el hecho de que su distribución no tenga forma normal no tendría que sorprendernos. El punto aquí es que \(s^2\) es insesgado, lo que es una propiedad independiente de la forma específica de su función de densidad.

Ahora vamos a aplicar el mismo procedimiento para \(s^2_n = \sum_{i=1}^n(x_i - \bar{x})^2/n\). Primero simulamos la distribución muestral tomando 100 mil muestras y en cada una obteniendo la varianza con \(s^2_n\).

## Función para muestra de n = 30 y tomar la varianza
mi.fun.var2 <- function(n){ 
  muestra <- rnorm(n)
  sum((muestra - mean(muestra))^2)/n # Estimador s^2_n
}
s2n <- replicate(100000, mi.fun.var2(30))

A continuación, tomamos el promedio de las 100 mil varianzas.

round(mean(s2n), digits = 3)
## [1] 0.967

El promedio de las varianzas calculadas con \(s^2_n\) no es igual a la varianza poblacional. Para ver que no se trata de un mero asunto de aproximación, repite el proceso unas cuantas veces. La cantidad obtenida siempre va a ser de alrededor de 0.967.

Entonces, independientemente del error en cada muestra particular, el estimador \(s^2_n\) subestima de manera el verdadero valor del parámetro \(\sigma^2\).

Este fenómeno es visualizado en la figura 3.2.

Distribución muestral de $s^2_n$

Figure 3.2: Distribución muestral de \(s^2_n\)

Esta es la razón por la que, a pesar de las ventajas intuitivas de \(s_n^2\) para describir la dispersión de la variable analizada en la muestra, \(s^2\) es preferido: tiene propiedades que se consideran deseables para fines inferenciales.

4 Explicando el sesgo

Lo anterior ilustra cómo la diferencia entre los estimadores \(s^2\) y \(s_n^2\) cuando tomamos el valor esperado. En esta sección daremos una explicación analítica del sesgo en \(s_n^2\), recurriendo a su formulación matemática.6

Para entender en qué consiste el sesgo, o, para ponerlo en términos coloquiales, qué “hace mal” \(s_n^2\), tomemos el valor esperado. Un primer punto a considerar es el siguiente: \[\mathrm{E}\left(\frac{\sum_{i=1}^n \left(x_i - \bar{x} \right)^2}{n} \right) = \frac{1}{n} \mathrm{E}\left(\sum_{i=1}^n \left(x_i - \bar{x} \right)^2\right);\] esto es, como el tamaño de la muestra \(n\) es constante, el valor esperado de \(s_n^2\) se obtiene multiplicando el valor esperado de su numerador por \(1/n\). Así,

\[\begin{align*} \mathrm{E}\left(\sum_{i=1}^n \left(x_i - \bar{x} \right)^2\right) & = \mathrm{E}\left( \sum_{i=1}^n x_i^2 - 2\bar{x}\sum_{i=1}^n x_i + \sum_{i=1}^n\bar{x}^2 \right) \\ & = \sum_{i=1}^n \mathrm{E}(x_i^2) - n \mathrm{E}(\bar{x}^2)\\ & = n\sigma^2 + n \mu^2 - \sigma^2 - n \mu^2\\ & = n\sigma^2 - \sigma^2\\ & = (n-1)\sigma^2. \end{align*}\]

Así, el valor estimado de \(s^2_n\) es \[\begin{align} \frac{1}{n} \mathrm{E}\left(\sum_{i=1}^n \left(x_i - \bar{x} \right)^2\right) = \frac{n-1}{n} \sigma^2. \tag{4.1} \end{align}\]

Entonces, el valor esperado de \(s^2_n\) subestima sistemáticamente el verdadero valor del parámetro \(\sigma^2\) al multiplicarlo por una cantidad menor a 1.7

Podemos ir un paso más allá en la interpretación del sesgo si descomponemos el resultado anterior: \[\begin{align} \frac{n-1}{n} \sigma^2 = \frac{n\sigma^2}{n} - \frac{\sigma^2}{n} = \sigma^2 - \frac{\sigma^2}{n} = \sigma^2 - \left(\frac{\sigma}{\sqrt{n}}\right)^2. \tag{4.2} \end{align}\]

Por el TCL, sabemos que \(\sigma/\sqrt{n}\) es el error estándar del estimador de la media. Esto nos habla del origen del sesgo: lo que \(s^2_n\) “hace mal” como estimador. El término \((x_i - \bar{x})^2\) toma la diferencia de cada observación con respecto a la media muestral. Sin embargo, la verdadera varianza es con relación a la media poblacional; esto es, \((x_i - \mu)^2\). Entonces, el numerador \(\sum_i (x_i - \bar{x})^2\) no toma en cuenta la diferencia entre \(\bar{x}\) y \(\mu\). La ecuación (4.2) muestra el impacto de no tomar en cuenta esta diferencia, pues al verdadero valor del parámetro se le descuenta el cuadrado del error estándar: esto es, la diferencia esperada entre la media poblacional \(\mu\) y la media muestral \(\bar{x}\).8

El estimador \(s^2\) tienen el mismo numerador que \(s^2_n\), pero corrige el sesgo con una ligera modificación en el denominador: \[\mathrm{E}\left(\frac{\sum_{i=1}^n \left(x_i - \bar{x} \right)^2}{n-1} \right) = \frac{1}{n-1} \mathrm{E}\left(\sum_{i=1}^n \left(x_i - \bar{x} \right)^2\right) = \frac{(n-1)\sigma^2}{n-1} = \sigma^2.\]

5 En conclusión

La discusión sobre el estimador de la varianza buscó cubrir los siguientes objetivos:

6 Para saber más…

El estimador de la varianza es explicado en prácticamente cualquier manual de estadística. La sustitución de \(n\) por \(n-1\) (corrección de Bessel) es discutida en la entrada de .

En la demostración para \(\mathrm{E}\left(\sum_{i=1}^n \left(x_i - \bar{x} \right)^2\right) = (n - 1)\sigma^2\) en la sección 4 se recurre a diversas propiedades tanto del valor esperado como de la varianza. En se exponen dichas propiedades de manera explícita (ecuacianes 3 a 7), lo que permite seguir más de cerca cada uno de los pasos.

7 Repaso

  1. En el ejemplo se asume la normal estándar, pero esto es sólo un caso especial escogido por su simplicidad. Para un sabor más “realista”, asumamos lo siguiente. Se toma al azar una muestra de 30 personas y se les aplica la prueba IQ. En la población, el IQ tiene una distribución normal con media \(\mu = 100\) y desviación estándar \(\sigma = 15\). Con la ecuación (4.1), calcula el valor que arrojaría el estimador sesgado. Simula las distribuciones muestrales de \(s^2\) y \(s^2_n\) y compara los promedios correspondientes con la cantidad previamente obtenida.

  2. De la misma forma, la distribución de la variable poblacional \(X\) no tiene que ser normal. Para ver esto, asume en \(X\) una distribución uniforme con valores entre 0 y 1. Simula las distribuciones muestrales de \(s^2\) y \(s^2_n\) y compara sus promedios con la varianza poblacional verdadera.

  3. Aunque el estimador no corregido \(s^2_n\) es sesgado, es consistente. Para ver de qué va esta propiedad, asume para la población \(X \sim N(0,1)\) y simula la distribución muestral de \(s^2_n\) varias veces incrementando sucesivamente el tamaño de la muestra. Describe qué pasa con el valor esperado y la forma de la distribución.

¿Cuál de las siguientes alternativas es preferible?

Argumenta tu respuesta.

  1. Nota para puristas: eso no es exactamente lo que da la desviación estándar, pero es una interpretación aceptada por su carácter intuitivo.↩︎

  2. Por supuesto, quien toma una muestra no conoce con precisión la distribución poblacional de la variable de interés. Para eso toma una muestra. Lo que hacemos aquí es explotar las posibilidades de la programación para simular estos objetos no observados y así entender mejor la teoría subyacente.↩︎

  3. Por defecto, la función rnorm() asume los argumentos mean = 0 y sd = 1.↩︎

  4. A fin de concentrarnos en el objetivo fundamental del ejercicio, presentamos las cantidades de interés redondeadas a tres digitos.↩︎

  5. Podríamos usar directamente la función var() en vez de sum((muestra - mean(muestra))^2)/(n-1), pero la idea es tener claro exactamente que estamos haciendo.↩︎

  6. En cierta medida, la sección es una respuesta a esta pregunta, planteada en un foro de discusión en línea y que encuentro bastante simpática: “¿Cómo es exactamente que los estadísticos acordaron usar \(n-1\) como el estimador insesgado de la varianza poblacional sin simulación?”.↩︎

  7. En nuestro ejemplo, \(\sigma^2 = 1\) y \(n = 30\). Compara lo que se obtiene al sustituir esos valores en la ecuación (4.1) con el resultado de la simulación.↩︎

  8. Claramente, este cuadrado es la varianza del estimador, pero seguimos la costumbre de expresar la variabilidad del estimador en términos de errores estándar.↩︎