2.7 Generalización del contraste sobre la matriz de covarianzas
En esta sección se va a generalizar el test obtenido para el contraste de la matriz de covarianzas. En primer lugar ofrecemos el enunciado de un teorema de generalización, cuya demostración es innecesaria, pues consiste en la constatación de los mismos argumentos de máxima verosimilitud ya empleados, y de los desarrollos subsiguientes.
Teorema 2.2 Sea \(\bm{x_1},\ldots,\bm{x_n}\) una muestra aleatoria simple de \(N_d(\bm{\mu},\bm\Sigma)\). Si las hipótesis \(H_0\) y \(H_a\) conducen a los estimadores de máxima verosimilitud \(\hat{\bm\Sigma}\) y \(\bm{S}\), respectivamente, y si \(\overline{\bm{x}}\) es el estimador de máxima verosimilitud para \(\bm{\mu}\) bajo cualquiera de las dos hipótesis, entonces el estadístico de razón de verosimilitudes para contrastar \(H_0\) frente a \(H_a\) viene dado por
\[-2\log\lambda(\bm{x})=nd \left(a-\log g -1\right)\] siendo \(a\) y \(g\) las medias aritmética y geométrica, respectivamente, de los autovalores de \(\hat{\bm\Sigma}^{-1}\bm{S}\).
Una vez establecido este resultado general, trataremos diversas situaciones en las cuales se puede aplicar este teorema. En todos los casos supondremos que disponemos de una muestra aleatoria simple \(\bm{x_1},\ldots,\bm{x_n}\in N_d(\bm{\mu},\bm\Sigma)\) y que el vector de medias \(\bm{\mu}\) es desconocido.
2.7.1 Contraste de la hipótesis nula \(H_0: \bm\Sigma=k\bm{\Sigma_0},\ \ k\in (0,+\infty)\)
Supongamos que queremos contrastar una hipótesis nula compuesta sobre la matriz de covarianzas
\[H_0: \bm\Sigma= k \bm{\Sigma_0}\qquad k\in (0,+\infty)\] siendo \(\bm{\Sigma_0}\) una matriz de covarianzas fijada, frente a una alternativa en la que la matriz de covarianzas no está sujeta a restricciones. El vector de medias carece de restricciones tanto bajo la hipótesis nula como bajo la alternativa. Estamos en las condiciones del teorema anterior, por lo que sólo nos falta calcular el estimador de la matriz de covarianzas bajo la hipótesis nula, \(\hat{\bm\Sigma}=\hat{k}\bm{\Sigma_0}\). Se puede comprobar que \(\hat{k}=a_0\) es el estimador de máxima verosimilitud de \(k\), siendo \(a_0\) la media aritmética de los autovalores de \(\bm{\Sigma_0}^{-1}\bm{S}\).
Entonces, aplicando la expresión que figura en el teorema anterior, el estadístico de razón de verosimilitudes adopta la forma:
\[-2\log\lambda(\bm{x})=nd \left(a-\log g -1\right)\] siendo \(a\) y \(g\) las medias aritmética y geométrica, respectivamente, de los autovalores de la matriz \(\hat{\bm\Sigma}^{-1}\bm{S}=\frac{1}{a_0}\bm{\Sigma_0}^{-1}\bm{S}\). En consecuencia, \(a=1\) y \(g=\frac{1}{a_0}g_0\), siendo \(g_0\) la media geométrica de los autovalores de \(\bm{\Sigma_0}^{-1}\bm{S}\). Sustituyendo los valores de \(a\) y \(g\) obtenemos
\[-2\log\lambda(\bm{x})=nd \left(1-\log\left(\frac{1}{a_0}g_0\right)-1\right) =nd\log\frac{a_0}{g_0}.\]
Por último, no estando disponible la distribución exacta de este estadístico, la aproximamos por una \(\chi_m^2\), siendo el número de grados de libertad \(m=\frac{1}{2}d(d+1)-1=\frac{1}{2}(d-1)(d+2)\).
2.7.1.1 Test de esfericidad \(H_0: \bm\Sigma=\sigma^2\bm{I_d}\)
Hay un caso particular de este tipo de contraste que tiene un interés especial. Es el test de esfericidad, que consiste en contrastar la hipótesis nula de que las variables son incorreladas y tienen la misma varianza. La hipótesis nula es de la forma
\[H_0: \bm\Sigma=\sigma^2\bm{I_d},\] siendo \(\sigma^2\) la varianza común desconocida. Nótese que la incorrelación equivale a independencia cuando se trata de variables normales. Es inmediato que estamos ante un caso particular del test anterior. Para verlo basta con tomar \(\bm{\Sigma_0}=\bm{I_d}\). Por tanto, el estadístico de contraste sería:
\[-2\log\lambda(\bm{x})=nd\log\frac{a_0}{g_0}\sim \chi_m^2\] siendo \(a_0\) y \(g_0\) las medias aritmética y geométrica, respectivamente, de los autovalores de la matriz \(\bm{\Sigma_0}^{-1}\bm{S}=\bm{S}\), y \(m=\frac{1}{2}(d-1)(d+2)\).
Ejemplo. Vamos a contrastar la esfericidad de datos generados de distintas
normales multivariantes. Para ello utilizaremos la función rmvnorm de
la librería mvtnorm. En primer lugar generamos una muestra de tamaño
\(n=100\) de una normal \(N_3(\bm{0},\bm{I_3})\).
> set.seed(123456)
> library(mvtnorm)
> d <- 3
> mu <- c(0, 0, 0)
> n <- 100
> Sigma <- diag(d)
> x <- rmvnorm(n, mu, Sigma)Una vez generados los datos, vamos a realizar el contraste de esfericidad. Para ello calculamos \(a_0\) y \(g_0\) las medias aritmética y geométrica, respectivamente, de los autovalores de la matriz \(\bm{S}\).
> S <- cov(x) * (n - 1)/n
> lambda <- eigen(S)$values
> a0 <- mean(lambda)
> g0 <- prod(lambda)^(1/d)
> estad <- n * d * log(a0/g0)
> estad
[1] 4.483253Bajo \(H_0\), la distribución asintótica del estadístico es \(\chi_m^2\), siendo \(m=(d-1)(d+2)/2\). En este caso, para un nivel de significación \(\alpha=0.05\), no rechazamos la hipótesis nula de esfericidad, ya que el valor del estadístico es menor que \(\chi_{m,\alpha}^2\). De hecho el valor crítico del estadístico es considerablemente grande.
Trabajaremos ahora con datos generados de una normal multivariante \(N_3(\bm{0},4\bm{I_3})\).
Planteamos el contraste de esfericidad y, de nuevo, no rechazamos la hipótesis nula.
> S <- cov(x) * (n - 1)/n
> lambda <- eigen(S)$values
> a0 <- mean(lambda)
> g0 <- prod(lambda)^(1/d)
> estad <- n * d * log(a0/g0)
> estad
[1] 3.251188Como último ejemplo, vamos a generar datos de un vector normal que no verifique la hipótesis nula de esfericidad. Para ello elegiremos una matriz de covarianzas diagonal pero en la que no todos las variables tengan la misma varianza.
Calculamos el estadístico de contraste y obtenemos los siguientes resultados.
> S <- cov(x) * (n - 1)/n
> lambda <- eigen(S)$values
> a0 <- mean(lambda)
> g0 <- prod(lambda)^(1/d)
> estad <- n * d * log(a0/g0)
> estad
[1] 69.92047La esfericidad permite una mayor o menor varianza, pero esta varianza tendría que ser la misma en las tres variables y no debe existir correlación. Esto se pone en evidencia mediante los autovalores de la matriz de covarianzas, que deben ser iguales. Como se puede comprobar, los autovalores para este último ejemplo muestran una gran desigualdad. Esto hace que las medias aritmética y geométrica sean muy diferentes. El nivel crítico es muy pequeño, por lo que hay pruebas muy significativas de ausencia de esfericidad.
2.7.2 Contraste de independencia de dos subvectores \(H_0: \bm{\Sigma_{12}}=\bm{0}\)
Supongamos que el vector de interés \(\bm{x}\) se divide en dos subvectores con \(d_1\) y \(d_2\) variables (\(d_1+d_2=d\)). Por ejemplo, si \(\bm{x_1}\) denota el subvector de las primeras \(d_1\) componentes de \(\bm{x}\) y \(\bm{x_2}\) denota el subvector con las restantes \(d_2\) componentes, entonces podemos escribir:
\[\bm{x}=\left( \begin{array}{c} \bm{x_1}\\ \bm{x_2} \end{array} \right),\ \ \ \bm{\mu}=\left(\begin{array}{c}\bm{\mu_1}\\ \bm{\mu_2}\end{array}\right),\ \ \ \bm{\Sigma}=\left(\begin{array}{cc}\bm{\Sigma_{11}}&\bm{\Sigma_{12}}\\ \bm{\Sigma_{21}}&\bm{\Sigma_{22}}\end{array}\right).\]
Queremos contrastar la hipótesis de que los dos conjuntos de variables son independientes entre sí, lo cual, en esta situación donde se supone normalidad, equivale a incorrelación, esto es, a que \(\bm{\Sigma_{12}}=\bm{0}\). Bajo la hipótesis nula \(H_0: \bm{\Sigma_{12}}=\bm{0}\), la verosimilitud se descompone en el producto de dos factores correspondientes a las verosimilitudes que provienen de cada conjunto de variables. De este modo, bajo la hipótesis nula, los estimadores de \(\bm{\mu_1}\) y \(\bm{\Sigma_{11}}\) por un lado, y de \(\bm{\mu_2}\) y \(\bm{\Sigma_{22}}\) por otro, se obtienen maximizando la verosimilitud de cada conjunto de variables por separado. Así, suponiendo que el vector de medias \(\bm{\mu}=(\bm{\mu_1},\bm{\mu_2})^\prime{}\) es desconocido, su estimador de máxima verosimilitud será \(\hat{\bm{\mu}}=(\bm{\overline{x}_1},\bm{\overline{x}_2})^\prime{}=\overline{\bm{x}}\), mientras que
\[\hat{\bm\Sigma}=\left(\begin{array}{cc} \bm{S_{11}} & \bm{0} \\ \bm{0} & \bm{S_{22}} \end{array}\right).\]
Bajo la alternativa, no hay restricciones ni para el vector de medias ni para la matriz de covarianzas, de modo que los estimadores de máxima verosimilitud serán \(\overline{\bm{x}}\) y \(\bm{S}\), respectivamente.
Entonces el estadístico de razón de verosimilitudes adopta la forma que figura en la expresión del teorema anterior, donde
\[\hat{\bm\Sigma}^{-1}\bm{S}=\left(\begin{array}{cc} \bm{S_{11}} & \bm{0} \\ \bm{0} & \bm{S_{22}} \end{array}\right)^{-1} \left(\begin{array}{cc} \bm{S_{11}} & \bm{S_{12}} \\ \bm{S_{21}} & \bm{S_{22}} \end{array}\right) =\left(\begin{array}{cc} \bm{I_{d_1}} & \bm{S_{11}}^{-1}\bm{S_{12}} \\ \bm{S_{22}}^{-1}\bm{S_{21}} & \bm{I_{d_2}} \end{array}\right).\]
En tal caso, \(a=\frac{1}{d}\mbox{traza}\left(\hat{\bm\Sigma}^{-1}\bm{S}\right)=\frac{1}{d}d=1\). Por otro lado,
\[ \begin{aligned} g^d&=\left|\hat{\bm\Sigma}^{-1}\bm{S}\right|\\ &=\frac{|\bm{S}|}{|\bm{S_{11}}|\cdot |\bm{S_{22}}|}\\ &=\frac{|\bm{S_{22}}-\bm{S_{21}}\bm{S_{11}}^{-1}\bm{S_{12}}|}{|\bm{S_{22}}|}\\ &=\left|\bm{I_{d_2}}-\bm{S_{22}}^{-1}\bm{S_{21}}\bm{S_{11}}^{-1}\bm{S_{12}}\right|\\ &=\left|\bm{I_{d_2}}-\bm{R_{22}}^{-1}\bm{R_{21}}\bm{R_{11}}^{-1}\bm{R_{12}}\right|. \end{aligned}\]
En el último paso hemos sustituido las matrices de covarianzas por matrices de correlaciones, que se construyen a partir de las anteriores así: \(\bm{R}=\bm{D^{-1/2}}\bm{S}\bm{D^{-1/2}}\), siendo \(\bm{D}\) una matriz diagonal que contiene las varianzas. Entonces, el estadístico de contraste será:
\[-2\log \lambda(\bm{x})=-nd\log g=-n\log\left|\bm{I_{d_2}}-\bm{S_{22}}^{-1}\bm{S_{21}}\bm{S_{11}}^{-1}\bm{S_{12}}\right|.\]
Ahora podríamos apelar a la distribución asintótica del estadístico \(-2\log \lambda(\bm{x})\), como hemos hecho en el resto de contrastes. Sin embargo, se puede obtener la distribución exacta \(\lambda(\bm{x})^{2/n}\) bajo \(H_0\). Se tiene que, si \(H_0: \bm{\Sigma_{12}}=\bm{0}\) es cierta, \(\lambda(\bm{x})^{2/n}\) sigue una distribución \(\Lambda\) de Wilks5 de parámetros \(d_2\), \(d_1\), \(n-1-d_1\). Es decir,
\[\lambda(\bm{x})^{2/n}=\frac{|\bm{S}|}{|\bm{S_{11}}|\cdot |\bm{S_{22}}|}=\left|\bm{I_{d_2}}-\bm{S_{22}}^{-1}\bm{S_{21}}\bm{S_{11}}^{-1}\bm{S_{12}}\right| \in \Lambda\left(d_2,d_1,n-1-d_1\right).\]
Como corresponde a un test de razón de verosimilitudes, se rechazará la hipótesis nula \(H_0: \bm{\Sigma_{12}}=\bm{0}\) cuando \(\lambda(\bm{x})\) sea pequeño, o equivalentemente, cuando \(|\bm{I_{d_2}}-\bm{S_{22}}^{-1}\bm{S_{21}}\bm{S_{11}}^{-1}\bm{S_{12}}|\) sea pequeño. La distribución \(\Lambda\) de Wilks se puede transformar (ver Anexo A) en una \(F\) de Snédecor6. Puesto que las funciones para transformar una \(\Lambda\) de Wilks en una \(F\) de Snédecor son decrecientes, se rechazará la hipótesis nula \(H_0: \bm{\Sigma_{12}}=\bm{0}\) cuando su transformación en una \(F\) de Snédecor sea grande.
Ejemplo. En el primer capítulo habíamos presentado un estudio sobre diabetes en el que se analizaban diferentes variables de interés sobre pacientes que no presentan la enfermedad. Las variables de mayor relevancia para el estudio eran:
\(x_1\)= Intolerancia a la glucosa
\(x_2\)= Respuesta de la insulina a la glucosa oral
\(x_3\)= Resistencia de la insulina
Se recogían igualmente otras dos variables de menor relevancia:
\(y_1\)= Peso relativo
\(y_2\)= Glucemia basal
Suponiendo que el vector \(\bm{v}=(y_1,y_2,x_1,x_2,x_3)^\prime{}\) sigue una distribución normal multivariante, nos interesa contrastar la independencia de los subvectores \(\bm{y}=(y_1,y_2)^\prime{}\) y \(\bm{x}=(x_1,x_2,x_3)^\prime\). Es decir, planteamos la hipótesis nula \(H_0: \bm{\Sigma_{yx}}=\bm{0}\). Como habíamos visto al comentar los datos, podemos particionar la matriz de covarianzas muestral \(\bm{S}\) y definir a partir de ella las submatrices \(\bm{S_{11}}\), \(\bm{S_{22}}\), \(\bm{S_{12}}\) y \(\bm{S_{21}}\).
> diab <- read.table("data/diabetes.txt", header = TRUE)
> n <- dim(diab)[1]
> S <- cov(diab) * (n - 1)/n
> S
y1 y2 x1 x2 x3
y1 0.01583006 0.2113327 0.7700567 -0.2091966 2.141484
y2 0.21133270 69.0250473 25.6587902 -23.4352552 -20.388469
x1 0.77005671 25.6587902 1082.3610586 388.1077505 106.027410
x2 -0.20919660 -23.4352552 388.1077505 2330.1025520 1117.797732
x3 2.14148393 -20.3884688 106.0274102 1117.7977316 2089.952741Calculamos el estadístico del contraste \(\lambda(\bm{v})^{2/n}\). Lo podemos hacer de dos formas:
Recordamos que si \(H_0\) es cierta, \(\lambda(\bm{v})^{2/n}\) sigue una distribución \(\Lambda\) de Wilks de parámetros \(d_2\), \(d_1\), \(n-1-d_1\). En nuestro ejemplo, \(d_2=3\), \(d_1=2\) y \(n=46\), por lo tanto, \(\lambda(\bm{v})^{2/n}\) sigue una distribución \(\Lambda(3,2,43)\). Se verifica además (ver Anexo A) que
\[\frac{41(1-\sqrt{\Lambda(3,2,43)})}{3\sqrt{\Lambda(3,2,43)}}\stackrel{d}{=} F_{6,82}\]
Teniendo en cuenta que rechazamos \(H_0\) cuando \(\lambda(\bm{x})\) es pequeño y que la transformación anterior es una función decreciente, rechararemos la hipótesis de independencia (para una significación \(\alpha\)) si
\[\frac{41(1-\sqrt{\Lambda(3,2,43)})}{3\sqrt{\Lambda(3,2,43)}}>f_{6,82,\alpha}\] siendo \(f_{6,82,\alpha}\) el cuantil \(1-\alpha\) de la distribución \(F_{6,82}\).
> d1 <- 2
> d2 <- 3
> # Parámetros de la distribución Lambda de Wilks (d,mH,mE)
> d <- d2
> mH <- d1
> mE <- n - 1 - d1
> # Transformación del estadístico a F de Snedecor
> estadF <- (mE - d + 1)/d * (1 - sqrt(estad))/sqrt(estad)
> estadF
[1] 2.394993> # Grados de libertad F de Snedecor
> glF1 <- 2 * d
> glF2 <- 2 * (mE - d + 1)
> qf(0.95, glF1, glF2)
[1] 2.211303Es decir, rechazamos la hipótesis nula de independencia.
2.7.2.1 Caso particular. Coeficiente de correlación múltiple
Consideremos que uno de los conjuntos de variables tenga un único elemento, por ejemplo, \(d_1=1\) y \(d_2=d-1\). En esta situación, \(\bm{R_{11}}=1\), y si denotamos \(\bm{\alpha}=\bm{R_{21}}\), éste será un vector (\(d-1\))–dimensional. De este modo, el estadístico de contraste resulta:
\[\lambda(\bm{x})^{2/n} = \left|{\bm{I_{d_2}}}-\bm{R_{22}}^{-1}\bm{\alpha}\bm{\alpha}^\prime{}\right| \stackrel{(a)}{=} 1-\bm{\alpha}^\prime{}\bm{R_{22}}^{-1}\bm{\alpha} = 1-\bm{R_{12}}\bm{R_{22}}^{-1}\bm{R_{21}} = 1-R^2\in \Lambda\left(d_2,1,n-2\right)\] siendo \(R\) el coeficiente de correlación múltiple entre la primera variable y las restantes. Usando que en general
\[\frac{1-\Lambda(d,1,m)}{\Lambda(d,1,m)}=\frac{d}{m-d+1}F_{d,m-d+1},\] tenemos que
\[\frac{R^2}{1-R^2}\in \frac{d-1}{n-d}F_{d-1,n-d}\]
Finalmente, el coeficiente de correlación múltiple se considerará significativo cuando el estadístico anterior \(R^2/(1-R^2)\) sea grande, comparado con la distribución \(F\) de Snédecor.
2.7.3 Contraste de independencia de todas las variables \(H_0:\bm\Sigma\) es diagonal
La hipótesis nula consiste en suponer que las variables son incorrelacionadas, pero, a diferencia del test de esfericidad, no exigimos que tengan la misma varianza. De nuevo, la incorrelación equivale a independencia en un contexto de normalidad. Así, bajo la hipótesis nula, se maximiza la verosimilitud separadamente para cada variable, dando lugar a los estimadores de la media y la varianza de dichas variables y en consecuencia a los estimadores del vector de medias y matriz de covarianzas:
\[\overline{\bm{x}}\qquad \mbox{y} \qquad \hat{\bm\Sigma}=\textnormal{diag}(s_{1}^2, \ldots, s_d^2)\]
Bajo la alternativa, no hay ninguna clase de restricciones sobre los parámetros \(\bm{\mu}\) y \(\bm\Sigma\) que, por tanto, admiten como estimadores de máxima verosimilitud \(\overline{\bm{x}}\) y \(\bm{S}\).
Entonces, aplicando el teorema general, el estadístico de contraste será:
\[-2\log \lambda(\bm{x})=nd\left(a-\log g-1\right)=-n\log |\bm{R}|\] siendo \(a\) y \(g\) las medias aritmética y geométrica, respectivamente, de los autovalores de la matriz \(\hat{\bm\Sigma}^{-1}\bm{S}\), y \(\bm{R}\) la matriz de correlaciones.
El último paso de la expresión anterior se debe a que \(\bm{R}=\hat{\bm{\Sigma}}^{\bm{-1/2}}\bm{S}\hat{\bm\Sigma}^{\bm{-1/2}}\) y esta matriz, aún siendo distinta de \(\hat{\bm\Sigma}^{-1}\bm{S}\), tiene los mismos autovalores que élla. Además, como la diagonal de una matriz de correlaciones está formada por unos, la traza vale \(d\) y, en consecuencia, la media de los autovalores vale uno, \(a=1\).
Por último, aproximamos la distribución del estadístico así:
\[-n\log|\bm{R}|\sim \chi_{\frac{1}{2}d(d-1)}^2\] donde el número de grados de libertad resulta de la diferencia del número de parámetros independientes bajo la hipótesis nula y bajo la alternativa: \(d+\frac{1}{2}d(d+1)-(d+d)=\frac{1}{2}d(d-1)\).
Ejemplo. Para ilustrar los contrastes de esfericidad habíamos generado, entre otros, datos de un vector normal multivariante con matriz de covarianzas diagonal y tal que no todas las variables tenían la misma varianza. En concreto, habíamos ejecutado el siguiente código para generar la muestra.
> library(mvtnorm)
> d <- 3
> mu <- c(0, 0, 0)
> Sigma <- diag(c(3, 5, 1))
> n <- 100
> x <- rmvnorm(n, mu, Sigma)Si queremos contrastar \(H_0:\bm\Sigma\) es diagonal, debemos calcular la matriz de correlaciones muestral y su determinante. Como hemos visto, bajo \(H_0\) la distribución asintótica del estadístico es \(\chi_m^2\), siendo \(m=d(d-1)/2\). En este caso, para un nivel de significación \(\alpha=0.05\), no rechazamos la hipótesis nula, ya que el valor del estadístico es menor que \(\chi_{m,\alpha}^2\). El valor crítico del estadístico nos lleva a aceptar sin lugar a dudas la hipotesis nula.
La distribución \(\Lambda\) de Wilks sirve para extender al caso multivariante la distribución \(F\) de Snédecor. Así, si la \(F\) de Snédecor se obtiene mediante el cociente de dos variables de tipo ji-cuadrado, la \(\Lambda\) de Wilks surge como el cociente de los determinantes de dos matrices de covarianzas.↩︎
En una \(\Lambda(d,m_H,m_E)\), si \(m_H=2\) (\(d\geq 2\)), \[\frac{1-\sqrt{\Lambda(d,2,m_E)}}{\sqrt{\Lambda(d,2,m_E)}}\stackrel{d}{=} \frac{d}{m_E-d+1} F_{2d,2(m_E-d+1)}.\]↩︎