La varianza generalizada

Introducción

La varianza generalizada (el determinante de la matriz de covarianzas) es una medida de la dispersión de una muestra de observaciones multivariantes. Vamos a comentar cuál es su distribución aproximada para datos normales. Es un resultado curioso porque es muy sencillo pero la aproximación requiere un conjunto de datos muy grande, especialmente si la dimensión de los datos es media o alta.

La varianza generalizada es una medida global de dispersión de una nube de puntos, es decir, una medida numérica única que resume la dispersión de los datos. Si disponemos de un conjunto de datos \(x_1,\ldots,x_n\) normales de dimensión \(p\) y llamamos \(S\) a la matriz de covarianzas muestral, la varianza generalizada de \(x_1,\ldots,x_n\) es el determinante de \(S\), que denotamos por \(|S|\).

Es ilustrativo considerar el caso \(p=2\). En esta situación, cada dato está formado por dos variables. Supongamos que las varianzas muestrales de las variables son \(s_1^2\), \(s_2^2\) y la correlación muestral entre ellas es \(r\). Tenemos que \[|S|= \left|\begin{array}{cc} s_1^2 & rs_1s_2\\ rs_1s_2 & s_2^2 \end{array}\right| = s_1^2s_2^2 (1-r^2).\] Si \(r\approx 0\), los puntos se sitúan aproximadamente en el interior de una elipse cuyos semiejes miden \(3s_1\) y \(3s_2\) (debido a que las variables son normales y aproximadamente incorreladas). Por tanto \(|S|^{1/2}\approx s_1s_2\) es proporcional al área ocupada por la nube de puntos. Para \(r\approx 1\), los puntos están aproximadamente alineados con lo que el área que ocupa la nube de puntos es muy pequeña. En este caso, \(|S|\approx 0\). En general, diagonalizando \(S\) se puede ver que para cualquier valor de \(r\) se mantiene esta interpretación geométrica de \(|S|^{1/2}\) en términos del área ocupada por los puntos.

Un inconveniente de \(|S|^{1/2}\) es que sus unidades son las del producto de las variables, lo que dificulta la comparación de dispersiones de conjuntos de datos con diferente dimensión. Para corregir este problema Peña y Rodríguez (2003) han propuesto la varianza efectiva \(|S|^{1/p}\) como medida de dispersión global. La varianza efectiva coincide con la media geométrica de los autovalores de \(S\).

La distribución aproximada de la varianza generalizada

Se puede demostrar el siguiente resultado, sorprendentemente simple (véase Anderson (2003), An Introduction to Multivariate Statistical Analysis, tercera edición, teorema 7.5.4.):

Supongamos que tenemos \(n\) v.a.i.i.d. con distribución normal multivariante \(\mbox{N}_p(\mu,\Sigma)\). Entonces, si \(p\) está fijo y \(n\to \infty\), \(\sqrt{n/p}\, (|S|/|\Sigma| - 1)\) converge en distribución a una normal de media cero y varianza 2.

Por tanto, la distribución aproximada del cociente entre la varianza generalizada muestral y la poblacional es muy sencilla: \[\frac{|S|}{|\Sigma|}\cong \mbox{N}\Big(1,\, \frac{2p}{n}\Big).\]

A partir de aquí resulta fácil obtener intervalos de confianza y contrastes aproximados para \(|\Sigma|\).

Desgraciadamente, la convergencia es extremadamente lenta para dimensiones moderadas. En el siguiente ejemplo, he generado 1000 muestras de 1000 datos de dimensión 20 cada una. (La media poblacional es cero y la matriz de covarianzas poblacional corresponde a 20 variables de varianza 1 y equicorreladas con \(\rho=0.5\).) Para cada muestra he calculado \(\sqrt{n/p}\, (|S|/|\Sigma| - 1)\). En un principio pensé que me había equivocado al observar que la distribución de los valores simulados de las varianzas generalizadas (línea continua gruesa) era tan diferente de \(\mbox{N}(0,2)\) (línea discontinua más fina):

library(MASS)
library(tidyverse)

# Parámetros
R <- 1000
n <- 1000
dimen <- 20
mu <- rep(0, dimen)
rho <- 0.5
sigma <- (1-rho) * diag(rep(1, dimen)) + rho * matrix(1, dimen, dimen)

# Genera muestras y calcula las varianzas generalizadas transformadas para 
# que la distribución sea N(0,2) aproximadamente
set.seed(100)
determinante <- NULL
for (i in 1:R){
  datos <- mvrnorm(n, mu, sigma)
  S <- cov(datos)
  determinante[i] <- sqrt(n/dimen)*(det(S)/det(sigma) - 1)
}
df <- data.frame(determinante)

# Compara la distribución asintótica con la de los valores generados
ggplot(df) +
  geom_density(aes(x = determinante, y=..density..), size=1.1) +
  geom_function(fun = dnorm, args = list(sd = sqrt(2)), linetype=2) +
  labs(x = "Varianza generalizada", y = "Densidad")

Sin embargo, si en vez de muestras de 1000 observaciones generamos muestras de 50000 observaciones, las distribuciones sí coinciden:

Esto significa que en dimensión 20 hay que disponer de bastantes miles de datos para que el resultado asintótico dé una aproximación adecuada.

La diferencia entre la distribución simulada y la aproximada se puede atribuir a la maldición de la dimensionalidad. De hecho, en dimensiones bajas la convergencia es mucho más rápida. Por ejemplo, si \(p=2\) y \(n=1000\), el resultado es:

La típica pregunta de los alumnos acerca de cuánto de grande tiene que ser el tamaño muestral para aplicar alguna aproximación no tiene respuesta general: depende de los valores de los parámetros (usualmente desconocidos) y, como en este ejemplo, de la dimensión de los datos. Podríamos decir que el tamaño importa, pero no sabemos en general cómo de grande debe ser… hay que ir caso a caso.

La distribución exacta

Como curiosidad, la distribución exacta de la varianza generalizada de datos normales es conocida y coincide con la de un producto de \(p\) variables \(\chi^2\) independientes. Según un resultado clásico debido a Wilks, 1934, se tiene \[(n-1)^p\frac{|S|}{|\Sigma|} \equiv \prod_{j=1}^p Y_{n-j+1},\] donde \(Y_n,\ldots,Y_{n-p+1}\) son independientes y \(Y_j\equiv \chi^2_{n-j+1}\).