4  Alternativas ante violaciones de supuestos

4.1 Matriz de diseño \(X\) con rango menor a \(p\)

Si \(X_{n \times p}\) tiene rango \(r<p\), no existe una única solución \(\widehat{\boldsymbol{\beta}}\) de las ecuaciones normales. En este caso tenemos tres formas de encontrar una solución \(\widehat{\boldsymbol{\beta}}\) y la proyección ortogonal \(\widehat{\boldsymbol{Y}}\):

1. Reducir el modelo para que \(X\) tenga rango completo.
2. Encontrar una inversa generalizada de \(\left(X^{\top} X\right)\).
3. Imponer condiciones para lograr identificabilidad de la solución a usar.

4.1.1 Reducir el modelo para que \(X\) tenga rango completo

La primera opción consiste en reducir el modelo para considerar las columnas linealmente independientes de \(X\). Si \(\operatorname{ran}(X)=r\), sea \(X_1\) la matrix con \(r\) columnas independientes de \(X\) y sea \(X_2\) la matriz con las columnas restantes de \(X\). Entonces \(X_2=X_1 F\) porque las columnas de \(X_2\) se pueden expresar como la combinación lineal de las columnas de \(X_1\).

\[ X=\left[X_1,\ X_2\right]=\left[X_1,\ X_1 F\,\right] =X_1\left[I_{r \times r},\ F\,\right] \]

Este es un caso especial de la factorización \(X= K L\), donde \(\operatorname{ran}\left(K_{n \times r}\right)=r\) y \(\operatorname{ran}\left(L_{r \times p}\right)=r\). Ahora, \[ \mathbb{E}[\boldsymbol{Y}]=X \boldsymbol{\beta}=K L \boldsymbol{\beta}=K \boldsymbol{\alpha}. \] Como \(K\) tiene rango completo, el estimador de mínimos cuadrados de \(\boldsymbol{\alpha}\) es \(\widehat{\boldsymbol{\alpha}}=\left(K^{\top} K\right)^{-1} K^{\top} \boldsymbol{Y}\) y la proyección ortogonal es \[ \widehat{\boldsymbol{Y}}=K \widehat{\boldsymbol{\alpha}}=K\left(K^{\top} K\right)^{-1} K^{\top} \boldsymbol{Y}. \] Entonces, \(P=K\left(K^{\top} K\right)^{-1} K^{\top}\), i.e. \(P=X_1\left(X_1^{\top} X_1\right)^{-1} X_1^{\top}\).

Ejemplo: One-way ANOVA con 2 categorías.

Consideramos \(n\) datos de los cuales \(n_1\) son categoría \(A\) y \(n_2\), \(B\). \[ \left(\begin{array}{c} \boldsymbol{Y}_1 \\ \boldsymbol{Y}_2 \\ \end{array}\right)= \left(\begin{array}{c} Y_{11} \\ \vdots \\ Y_{1 n_1} \\ Y_{21} \\ \vdots \\ Y_{2 n_2} \end{array}\right)=\left(\begin{array}{ccc} 1 & 1 & 0 \\ \vdots & \vdots & \vdots \\ 1 & 1 & 0 \\ 1 & 0 & 1 \\ \vdots & \vdots & \vdots \\ 1 & 0 & 1 \end{array}\right)\left(\begin{array}{c} \mu \\ \beta_1 \\ \beta_2 \end{array}\right)+\left(\begin{array}{c} \varepsilon_{11} \\ \vdots \\ \varepsilon_{1 n_1} \\ \varepsilon_{21} \\ \vdots \\ \varepsilon_{2 n_2} \end{array}\right) \] Sea \(X_1\) la matriz con las primeras dos columnas de \(X\). Entonces, \[ X=X_1\left[I_{r \times r},\ F\,\right] = X_1\left(\begin{array}{rrr} 1 & 0 & 1 \\ 0 & 1 & -1 \end{array}\right) \] y \[ X \boldsymbol{\beta}=X_1 \boldsymbol{\alpha},\quad \text{ donde }\quad \boldsymbol{\alpha}=\left(\begin{array}{rrr} 1 & 0 & 1 \\ 0 & 1 & -1 \end{array}\right)\left(\begin{array}{c} \mu \\ \beta_1\\ \beta_2\end{array}\right)=\binom{\mu+\beta_2}{\beta_1-\beta_2}. \] Entonces, \[ \begin{aligned} \widehat{\boldsymbol{\alpha}}&=\left(X_1^{\top} X_1\right)^{-1} X_1^{\top} \boldsymbol{Y}\\[10pt] & =\left(\begin{array}{cc} n & n_1 \\ n_1 & n_1 \end{array}\right)^{-1}\binom{\sum_j Y_{1 j}+\sum_j Y_{2 j}}{\sum_j Y_{1 j}}\\[10pt] &=\frac{1}{(n-n_1)\,n_1}\left(\begin{array}{cc} n_1 & -n_1 \\ -n_1 & n \end{array}\right)\binom{\sum_j Y_{1 j}+\sum_j Y_{2 j}}{\sum_j Y_{1 j}} \\[10pt] &=\left(\begin{array}{cc} n_2^{-1} & -n_2^{-1} \\ -n_2^{-1} & n_1^{-1}+n_2^{-1} \end{array}\right)\binom{\sum_j Y_{1 j}+\sum_j Y_{2 j}}{\sum_j Y_{1 j}} \\[10pt] & =\binom{\overline{\boldsymbol{Y}}_2}{\overline{\boldsymbol{Y}}_1-\overline{\boldsymbol{Y}}_2} \end{aligned} \] y entonces \(\widehat{\boldsymbol{Y}}=X_1 \widehat{\boldsymbol{\alpha}}=\left(\overline{\boldsymbol{Y}}_1, \ldots, \overline{\boldsymbol{Y}}_1, \overline{\boldsymbol{Y}}_2, \ldots, \overline{\boldsymbol{Y}}_2 .\right)^{\top}\).

4.1.2 Encontrar una inversa generalizada

En matemáticas, y en particular, álgebra, una generalizada inversa (o, g-inversa) de \(x\) es un elemento \(y\) que tiene algunas propiedades de un elemento inverso pero no necesariamente todos ellos.

El propósito de construir un inverso generalizado de una matriz es obtener una matriz que pueda servir como inversa en algún sentido para una clase más amplia de matrices.

Definición 4.1 (Inversa Generalizada) Una matriz \(A^{\mathrm{-}} \in \mathbb{R}^{n \times m}\) es una inversa generalizada de \(A \in \mathbb{R}^{m \times n}\) si \[ A A^{\mathrm{-}} A=A \]

Siempre (ver Teorema 4.2) se puede definir una inversa generalizada de una matriz arbitraria, pero cuando una matriz tiene un inverso regular, este inverso debe coincidir con la inversa común. La matriz \(A^{-1}\) ha sido denominada como inversa regular por algunos autores.

Consideramos el sistema \[ A x=y \] donde \(A\) es una \(m \times n\) matriz y \(y \in \mathcal{C}(A)\), el espacio de la columna de \(A\). Si \(m = n\) y \(A\) es no singular, entonces \(x = A^{-1} y\) será la solución del sistema. En particular tenemos que \[ A A^{-1} A=A. \]

Ahora supongamos \(A\) que es rectangular \((m \neq n)\), o cuadrada y singular. Entonces necesitamos un candidato \(G\) de dimensión \(n \times m\) tal que para todos \(y \in \mathcal{C}(A)\), \[ A G y=y, \] para que \(x=G y\) sea una solución del sistema lineal \(A x = y\).

Equivalentemente, necesitamos una matriz \(G\) de orden \(n \times m\) de tal manera que \[ A G A=A \] (ya que si \(AGA=A\) entonces \(AGy=AGAx=Ax\), pero \(Ax=y\)  y \(\ \therefore \ \ AGy=y\)).

Los tipos importantes de inversa generalizada incluyen:

Si la matriz \(A\) es \(m \times n\)

• Inversos un parcial (inverso derecho o inverso izquierdo)
  1. Inverso derecho: Si \(\operatorname{ran}(A)=m\), entonces existe una \(n \times m\) matriz \(A_{\mathrm{R}}^{-1}\) llamada la derecha inversa de \(A\) tal que \(A A_{\mathrm{R}}^{-1}=I_m\).
  2. Inverso izquierdo: Si \(\operatorname{ran}(A)=n\), entonces existe una \(n \times m\) matriz \(A_{\mathrm{L}}^{-1}\) llamada la izquierda inversa de \(A\) tal que \(A_{\mathrm{L}}^{-1} A=I_n\).
• Inversa de Bott-Duffin. Introducida en Bott y Duffin (1953).
• Inversa de Drazin. Introducida en Drazin (1958).
• Inversa de Moore-Penrose (también conocida como “pseudoinversa” o “p-inversa”).

Algunos inversos generalizados se definen y clasifican en función de las condiciones Penrose:

1. \(A A^{-} A=A\)
2. \(A^{-} A A^{-}=A^{-}\)
3. \(\left(A A^{-}\right)^*=A A^{-}\)
4. \(\left(A^{-} A\right)^*=A^{-} A\),

donde \({ }^*\) denota la transpuesta Hermitiana (transpuesta conjugada).

• Si \(A^{-}\) satisface la primera condición, entonces es una inversa generalizada de \(A\). Si cumple las dos primeras condiciones, entonces es una inversa generalizada reflexiva de \(A\).

• Si cumple las cuatro condiciones, entonces es la inversa de Moore-Penrose \(A\), que se denota como \(A^{+}\). El lector interesado, puede referirse a Thome, N. (2019). Esta inversa recibe su nombre por los matemáticos E. H. Moore y Roger Penrose, quienes descubrieron y estudiaron sus propiedades de manera independiente.

Como comentamos arriba, cuando \(A\) es no singular, cualquier inverso generalizado \(A^{-}=A^{-1}\) y por lo tanto es único.

Para una matriz \(A\) singular, algunos inversos generalizados, como las inversas de Drazin y de Moore-Penrose, son únicas; mientras que otras no están necesariamente definidas de manera única.

Teorema 4.1 (Unicidad de la inversa) Si \(A\) es no singular, la única inversa generalizada de \(A\) es \(A^{-1}\).

Demostración

Por un lado, \(A A^{-1} A=I A=A\), entonces \(A^{-1}\) es una inversa generalizada.

Ahora, si \(A A^{-} A=A\), entonces \(A A^{-}=A A^{-}I = A A^{-} A A^{-1}=A A^{-1}=I\).  Así \(A^{-}\) es la inversa de \(A\).

Ejemplo: Una inversa generalizada

Un método para obtener una inversa generalizada de \(B\) es el siguiente:

1. Si \(\operatorname{ran}(B)=r<\operatorname{mín}(m,n)\), eliminamos \(m-r\) filas y \(n-r\) columnas para dejar una \(r \times r\) matriz que sea no singular.

2. Invertimos esta matriz cuadrada \(r \times r\).

3. Obtenemos \(B^{-}\) insertando bloques de renglones y columnas a la inversa obtenida para que esta nueva matriz tenga \(n\) renglones y \(m\) columnas.

Por ejemplo si \[ B=\left(\begin{array}{ll} B_{11} & B_{12} \\ B_{21} & B_{22} \end{array}\right) \] y \(B_{11}\) es una matriz \(r \times r\) no singular, entonces \[ B^{-}=\left(\begin{array}{cc} B_{11}^{-1} & 0 \\ 0 & 0 \end{array}\right) \]

Omitimos la demostración pero ilustramos con un ejemplo numérico

Ejemplo Numérico para la inversa generalizada

\[ B=\left[\begin{array}{lll} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{array}\right] \]

Como \(\operatorname{ran}(B)=2\) y el bloque superior izquierdo es invertible, podemos encontrar entonces fácilmente la inversa generalizada de \(B\) \[ B^{-}=\left[\begin{array}{ccc} -\frac{5}{3} & \frac{2}{3} & 0 \\ \frac{4}{3} & -\frac{1}{3} & 0 \\ 0 & 0 & 0 \end{array}\right] \] Verificamos que cumpla con la propiedad de inverso generalizado: \[ B B^{-} B=\left[\begin{array}{lll} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{array}\right]\left[\begin{array}{ccc} -\frac{5}{3} & \frac{2}{3} & 0 \\ \frac{4}{3} & -\frac{1}{3} & 0 \\ 0 & 0 & 0 \end{array}\right]\left[\begin{array}{lll} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{array}\right]=\left[\begin{array}{lll} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{array}\right]=B \]

Lema 4.1 (Una inversa generalizada para matrices simétricas) Sea \(A\) una matriz simétrica, entonces se puede definir una matriz inversa generalizada de \(A\) a partir de los elementos de la descomposición espectral de \(A\).

Demostración

Como \(A\) simétrica, entonces existe \(T\) ortogonal tal que \(T^{\top} A T=\Lambda\) y \(A=T\Lambda T^{\top}\).

Sea \(\Gamma=\operatorname{diag}(\gamma_1,\ldots,\gamma_p),\) donde \[ \gamma_i= \begin{cases}1 / \lambda_i, & \text { si }\quad \lambda_i \neq 0 \\ 0, & \text { si }\quad \lambda_i=0\end{cases} \] y \(G = T\ \Gamma \ T^{\top}\).

Mostramos ahora que \(G\) es una inversa generalizada de \(A\) . Como \(T\) es ortogonal, \(T^{\top} T=I\) y \[ \begin{aligned} A G A & =T \Lambda T^{\top} T\ \Gamma\ T^{\top} T \Lambda T^{\top} \\ & =T \Lambda \ \Gamma \ \Lambda T^{\top} \\ & =T \Lambda T^{\top} \\ & =A \qquad \qquad \because \quad A \text{ es simétrica} \end{aligned} \]

Es trivial demostrar que esta inversa generalizada es reflexiva.

El resultado anterior se puede generalizar para cualquier matriz real rectangular usando la descomposición SVD (?thm-svd):

Teorema 4.2 (Existencia de la inversa generalizada) Una inversa generalizada de una matriz \(A\) que considera la descomposición en valores singulares de \(A\) (?thm-svd): \[ \underset{m \times n}{A}=\underset{m\times m}{U}\ \ \underset{m \times n}{\Sigma}\ \ \underset{n\times n}{V^{\top} } \] es \[ \underset{n\times m}{A^{-}}=\underset{n\times n}{V}\ \ \underset{n\times m}{\Sigma^{+}}\ \ \underset{m\times m}{U^{\top}}\ , \] donde \(\Sigma^{+}\) tiene en su diagonal los valores \(\sigma_i\) cuando \(\sigma_i>0\) y 0 cuando \(\sigma_i=1\).

Ejemplo: Una inversa generalizada

Sea \[ A = \begin{pmatrix} 1 & 2 \\ 3 & 4 \\ 5 & 6 \end{pmatrix} \]

La descomposición en valores singulares de \(A\) es:

\[ A = U \Sigma V^{\top} \approx \begin{pmatrix} -0.2298 & 0.8835 \\ -0.5247 & 0.2405 \\ -0.8196 & -0.4016 \end{pmatrix} \begin{pmatrix} 9.51 & 0 \\ 0 & 0.71 \\ 0 & 0 \end{pmatrix} \begin{pmatrix} -0.619 & -0.785 \\ -0.785 & 0.619 \end{pmatrix}. \] Ahora, \[ \Sigma^+ = \begin{pmatrix} 1/9.52 & 0 & 0\\ 0 & 1/0.51 & 0\\ \end{pmatrix} = \begin{pmatrix} 0.105 & 0 & 0 \\ 0 & 1.944 & 0 \\ \end{pmatrix} \] Así que la inversa generalizada de \(A\) es: \[ \begin{align} A^- &= V \Sigma^+ U^{\top} \\ &=\begin{pmatrix} 0.619 & -0.785 \\ 0.785 & 0.619 \end{pmatrix} \begin{pmatrix} 0.105 & 0 & 0\\ 0 & 1.408 & 0\\ \end{pmatrix} \begin{pmatrix} -0.229 & -0.5247 & -0.8196 \\ 0.883 & 0.2405 & -0.4016 \\ 0.408 & -0.816 & 0.4082 \\ \end{pmatrix}\\ &\approx \begin{pmatrix} -1.333 & -0.406 & 0.560 \\ 1.045 & 0.247 & -0.552 \end{pmatrix} \end{align} \]

Lema 4.2 Si \(G\) y \(H\) son dos inversas generalizadas de \(\left(X^{\top} X\right)\), entonces

1. \(X G X^{\top} X = X H X^{\top} X = X\),
2. \(X G X^{\top} = X H X^{\top}\).

La demostración se deja como ejercicio.

En la regresión, si la matriz \(n\times p\) de diseño \(X\) tiene un rango \(r<p\), las ecuaciones normales \(X^{\top} X \boldsymbol{\beta}=X^{\top} \boldsymbol{Y}\) no tienen una solución única. Estableciendo \(B=X^{\top} X\) y \(\boldsymbol{c}=X^{\top} \boldsymbol{Y}\), tenemos \[ \begin{aligned} \boldsymbol{c} & =B \boldsymbol{\beta} \\ \Leftrightarrow \qquad B\left(B^{-} \boldsymbol{c}\right)&=B B^{-} B \boldsymbol{\beta} \end{aligned} \] y entonces \(B^{-} \boldsymbol{c}\) es una solución de \(B \boldsymbol{\beta}=\boldsymbol{c}\). (De hecho, se puede demostrar que cada solución de \(B \boldsymbol{\beta}=\boldsymbol{c}\) puede expresarse en la forma \(B^{-} \boldsymbol{c}\) para algún \(B^{-}\)).

Como hemos visto, hay varias formas de calcular un \(B^{-}\) En particular, Thome (2019) describe algunos métodos para calcular \(B^{+}\) para la matriz simétrica \(B\). Algunos son:

1. Método basado en la pseudoinversa,
2. La descomposición en valores singulares,
3. Método basado en el algoritmo de Gauss-Jordan,
4. Método de Greville.

Teorema 4.3 (Proyección bajo inversa generalizada) \[ X\left(X^{\top} X\right)^{-} X^{\top} \] es el operador de proyección perpendicular en \(\mathcal{C}(X)\).

Demostración

Necesitamos demostrar las siguientes dos condiciones para \[ M=X\left(X^{\top} X\right)^{-} X^{\top} \]

1. \(v \in \mathcal{C}(X)\) implica \(M v=v \quad\) (proyección),
2. \(w \in \mathcal{C}^{\perp}(X)\) implica \(M w=0 \quad\) (perpendicularidad).

En relación a ii. sabemos que en particular, una proyección es ortogonal si \(M^{\top}(I-M)=0\)

(i) Sea \(v \in \mathcal{C}(X)\), entonces podemos escribir a \(v\) como \(v=X b\) y entonces tenemos \[ X\left(X^{\top} X\right)^{-} X^{\top} v = X\left(X^{\top} X\right)^{-} X^{\top} X b = X b = v. \]

(ii)

Queremos mostrar que \(M^{T}(I-M)=0\). \[ \left[X\left(X^{\top} X\right)^{-} X^{\top}\right]\left[I-X\left(X^{\top} X\right)^{-} X^{\top}\right] = X\left(X^{\top} X\right)^{-} X^{\top}-X\left(X^{\top} X\right)^{-} X^{\top}=0. \]

Entonces el estimador de mínimos cuadrados, considerando la inversa generalizada \((X^{\top}X)^{-}\) es \[ \begin{aligned} \widehat{\boldsymbol{\beta}} & =(X^{\top}X)^{-} X^{\top} \boldsymbol{Y} \\ \end{aligned} \] Este valor es una solución de las ecuaciones normales.

En particular, \(X^{+} \mathbf{Y}\) es la única solución que minimiza \(\widehat{\boldsymbol{\beta}}^{\top} \widehat{\boldsymbol{\beta}}\) (Peters and Wilkinson [1970]).

Finalmente notamos que \(\widehat{\boldsymbol{\theta}}=X \widehat{\boldsymbol{\beta}}=X\left(X^{\top} X\right)^{-} X^{\top} \boldsymbol{Y}\) y \(P= X\left(X^{\top} X\right)^{-} X^{\top}\) es matriz de proyección de \(\mathbb{R}^n\) a \(\mathcal{C}(\mathbf{X})\).

4.1.3 Imponer condiciones para lograr identificabilidad de la solución a usar

En ocasiones cuando los parámetros \(\beta\)’s tienen un significado importante y no queremos remover columnas de \(X\), queremos producir una solución única imponiendo algunas restricciones en \(\beta\).

Imponemos \(p-r\) restricciones en \(\boldsymbol{\beta}\) para hacer que \(\boldsymbol{\beta}\) sea determinada en forma única (identificable), i.e. tal que para cualquier \(\boldsymbol{\theta} \in \mathcal{C}(X)\), haya una único \(\boldsymbol{\beta}\) que satisface \[ X \boldsymbol{\beta}=\boldsymbol{\theta} \quad \text { y } \quad H \boldsymbol{\beta}=\mathbf{0} . \]

Estas ecuaciones se puede escribir como \[ \binom{\boldsymbol{\theta}}{\mathbf{0}}=\binom{X}{H} \boldsymbol{\beta} =: G \boldsymbol{\beta}. \tag{4.1}\]

Teorema 4.4 La expresión Ecuación 4.1 tiene una solución única ssi \(G\) tiene rango \(p\) y los renglones de \(H\) son linealmente independientes de los renglones de \(X\).

Teorema 4.5 La expresión Ecuación 4.1 tiene una solución única ssi \(G\) tiene rango \(p\) y \(H\) tiene rango \(p-r\).

Para estimar \(\boldsymbol{\beta}\), resolvemos \(\widehat{\boldsymbol{Y}}=X \widehat{\boldsymbol{\beta}}\)   y   \(H \widehat{\boldsymbol{\beta}}=\mathbf{0}\), es decir, resolvemos las ecuaciones normales aumentadas: \[ X^{\top} X \widehat{\boldsymbol{\beta}}= X^{\top} \boldsymbol{Y} \qquad \text{ y }\qquad H^{\top} H \widehat{\boldsymbol{\beta}}=\mathbf{0} \] esto es, \[ \binom{X^{\top} X}{H^{\top}H} \widehat{\boldsymbol{\beta}}=\left(G^{\top}G\right) \widehat{\boldsymbol{\beta}} =\binom{ X^{\top}\boldsymbol{Y}}{\boldsymbol{0}} \]

Entonces, \[ \widehat{\boldsymbol{\beta}}=\left(G^{\top} G\right)^{-1} X^{\top} \boldsymbol{Y} \quad \text { y } \quad \widehat{\boldsymbol{Y}}=X \widehat{\boldsymbol{\beta}}=P \boldsymbol{Y} \] donde \(P=X\left(G^{\top}G\right)^{-1} X^{\top}\).

Ejemplo: One-way ANOVA con 2 categorías (continuación).

Imponemos la condición \(\beta_1+\beta_2=0\), esto es, \[ H \boldsymbol{\beta} \equiv (0,1,1)\left(\begin{array}{c} \mu \\ \beta_1 \\ \beta_2 \end{array}\right)=0. \] Cuando \(n_1=n_2\), entonces \[ \widehat{\boldsymbol{\beta}}=\left(\begin{array}{c} \overline{Y}_{. .} \\ \frac{1}{2}\left(\overline{Y}_{1 .}-\overline{Y}_{2 .}\right) \\ \frac{1}{2}\left(\overline{Y}_{2 .}-\overline{Y}_{1 .}\right) \end{array}\right) \] satisface las ecuaciones normales y la restricción \(\beta_1+\beta_2=0\). Entonces, como antes, tenemos \(\widehat{\boldsymbol{Y}} = X \widehat{\boldsymbol{\beta}}=\left(\overline{Y}_1 ., \ldots, \overline{Y}_1, \overline{Y}_2, \ldots, \overline{Y}_2\right)^{\top}.\)

4.2 Multicolinearidad

Como comentamos en el capítulo anterior, aunque con la multicolinearidad no se viola el supuesto que el rango de la matriz de diseño sea completa, la cercana dependencia lineal de algunos matrices de datos de al manos un par de covariables puede ocasionar una disminución importante en la calidad e interpretación de los resultados.

Algunos posible métodos para disminuir o eliminar la multicolinearidad son:

1. Aumentar el tamaño de la muestra. Ampliar la muestra introducirá más variación en la serie de datos, lo que reduce el efecto del error de muestreo y ayuda a aumentar la precisión al estimar varias propiedades de los datos.

2. Eliminar algunas de las características altamente correlacionadas.

   • Método manual: Como el Factor de Inflación de Varianza (VIF) mide cuánto aumenta la varianza de un coeficiente de regresión debido a la colinealidad con otras variables. Un VIF alto indica una alta colinealidad. En este método, calculamos el VIF para cada covariable y se considera para eliminar si su VIF excede un umbral (por ejemplo, 10).

   • Eliminación de características usando Descomposición PCA. PCA es una técnica de reducción de dimensionalidad que transforma las características originales en un nuevo conjunto de características no correlacionadas llamadas componentes principales.

3. Reemplazar los regresores altamente correlacionados con una combinación lineal de ellos. Por ejemplo en el caso de la regresión de la esperanza de vida, en lugar de considerar el tamaño de la población y la extensión del estado, si estuvieran fuertemente correlacionadas, se puede optar por usar la nueva variable °°densidad poblacional**.

4. Usar métodos de regularización como RIDGE o LASSO. Estos métodos buscan ajustar un modelo con las covariables más informativas, así que las covariables redundantes, por ser altamente correlacionadas con las incluidas, tenderán a ser eliminadas.

4.2.1 Contracción de \(X^{\top}X\) a una matriz de menor rango

Este enfoque para la estabilización de las estimaciones de regresión se basa en la contracción de \(X^{\top} X\) a una matriz de menor rango. Uno de los métodos para hacer esto consiste en expandir \(X^{\top} X\) con su descomposición espectral y luego reter solo los primeros \(r\) de estos \(p\) vectores propios para representar \(X^{\top} X\), \((r<p)\). El estimador de mínimos cuadrados se obtiene entonces con una inversa generalizada.

Los \(r\) primeros vectores singulares pueden entonces considerarse como nuevas variables predictoras independientes. Al modelo derivado de este procedimiento se le conoce como la “regresión de componentes principales”, pero es differente a lo que ahora se conoce ahora más ampliamente como “Regresión por componentes principales” que revisamos en la siguiente sección.

Consideremos el modelo (con variables estandarizadas) \(W\): \[ y = \alpha \boldsymbol{1}_n + W \boldsymbol{\beta} + \boldsymbol{\epsilon} \] y hemos comentado que los estimadores de mínimos cuadrados para los parámetros son: \[ \widehat{\alpha}=\overline{y} \quad \text{y} \quad \widehat{\boldsymbol{\beta}} = \left(W^{\top}W\right)^{-1}W^{\top}\boldsymbol{y}=\left(\sum_{j=1}^k \frac{1}{\lambda_j} \boldsymbol{v}_j \boldsymbol{v}_j^{\top} \right)W^{\top} \boldsymbol{y}. \] Ahora, si los que dan problemas son los términos en los que las \(\lambda_j\)’s son pequeñas, una opción que puede ser natural es la de eliminarlos. Supongamos que hay dos \(\lambda_j\)’s bastante más chicas que las demás (las últimas dos, ya que el vector de eigenvalores está ordenado en forma descendente), entonces bajo este enfoque, el estimador de \(\boldsymbol{\beta}\) es : \[ \widehat{\boldsymbol{\beta}}_{(-)} = \left(\sum_{j=1}^{k-2} \frac{1}{\lambda_j} \boldsymbol{v}_j \boldsymbol{v}_j^{\top} \right)W^{\top} \boldsymbol{y} = (W^{\top}W)^{-}_*W^{\top}\boldsymbol{y}, \] con la inversa generalizada definida como en Lema 4.2, pero haciendo los dos menores eigenvalores iguales a ser cero (por eso el asterisco como subíndice). Con esto, estamos considerando que \(W\) tiene rango \(k-2<k\) pero no eliminamos ninguna covariable para obtener \(\widehat{\boldsymbol{\beta}}_{(-)}\).

Ejemplo: Medidas corporales

Ejemplificamos este estimador con los datos de mediciones físicas. Los datos son bajados de https://opendata.usm.my/items/564b9a77-89a1-453b-b6ed-2697f32e3d73 y contienen los siguientes campos:

• Gender (Male and Female (M=1 & F= 2) (391 Males & 324 Females)
  
• Age (1 year and above)
  
• HeadCircumference (in inches)
  
• ShoulderWidth (in inches)
• ChestWidth (in inches)
  
• Belly (in inches)
• Waist (in inches)
• Hips (in inches)
• ArmLength (in inches)
• ShoulderToWaist (in inches)
• WaistToKnee (in inches)
• LegLength (in inches)
• TotalHeight - from head to toe (in inches)
• Class Label (Not defined)

library(knitr)
dat <- read.csv("data/bjv6c9pmp4-1/Body_Measurements_original_CSV.csv", head=TRUE)
#Hacemos m'as cortos los nombres
names(dat)<-c("Gender","Age","Head_Circ", "Should_WD", "Chest_WD", "Belly",
              "Waist","Hips", "Arm_LEN","Should_Waist", "Waist_Knee", "Leg_LEN",
              "T_Height")
kable(head(dat))

Gender	Age	Head_Circ	Should_WD	Chest_WD	Belly	Waist	Hips	Arm_LEN	Should_Waist	Waist_Knee	Leg_LEN	T_Height
1	30	22	18	20	18	14	22	22	25	25	22	52
1	28	19	22	17	18	21	25	28	23	25	20	56
2	27	21	18	16	14	10	15	21	18	14	18	53
1	29	20	20	18	11	19	14	24	21	20	21	45
2	28	16	14	18	13	11	30	25	22	32	13	47
2	22	17	19	18	14	16	18	20	24	21	19	60

library(ggplot2)
library(GGally)
ggpairs(dat,columns = 2:13, title = "Matrix Plot")

library(corrplot)
# Calcular la matriz de correlación
cor_matrix <- cor(subset(dat, select = -c(Gender)))
# Generar la gráfica de matriz de correlación
corrplot(cor_matrix, method = "circle", type = 'lower',title = "")

library(regclass)
dat1 <- subset(dat, select = -c(Gender))
mod <- lm(T_Height~., data=dat1)
par(oma=c(4,2,2,0))
plot(VIF(mod),t="h",lwd=2,axes = FALSE,xlab="",main="FIV")
points(VIF(mod),pch=19,cex=2)
grid()

# Add Y-axis
axis(2)

# Add X-axis
axis(1,
     at=seq(1,11),
     labels = names(dat1)[-12],
     col = "blue",       # Axis line color
     col.ticks = "blue", # Ticks color
     col.axis = "blue",  # Labels color
     las=3)    

Estamos interesados en obtener un modelo para la altura total (T_Height) y notamos que hay una correlación positiva entre varias medidas corporales, como Age-Shoud_WD-Should-Waist-Waist_Knee, Waist-Hips.

X    <- as.matrix(dat1[,-12])
n    <- dim(X)[1]
k    <- dim(X)[2]
y    <- dat1[,12]

# Estandarizamos los datos:
Cen  <- (diag(n)-matrix(1,n,n)/n) %*% X                  # Centrar
W    <- Cen %*% diag(1/sqrt(diag(t(Cen) %*% Cen)/(n-1))) # Escalamos
#Las dos lineas anteriores equivalen a 
# W  <- scale(X)

# Descomposición espectral
WtW   <- t(W) %*% W
eig   <- eigen(WtW)
vec   <- (eig$vectors)[,1:(k-2)]  #Elimnamos los dos eigenvalores más pequeños
val   <- (eig$values)[1:(k-2)]
WTWm  <- vec %*% diag(1/val) %*% t(vec) #Descomposici'on Espectral

# Estimaci'on de beta_CP
bcp   <- WTWm %*% t(W) %*% y

vifp  <- diag(WTWm)   # En la sig prop se ve  por qué los vif's se calculan así
MSres <- sum((y-mean(y)-W%*%(as.vector(bcp)))^2) / (n-k-1)  # MS_Res

# Contratamos con valores de beta en la regresion regular
# los vectores contienen las ordenadas al origen
bcpa      <- c(mean(y), bcp)
dat1_est  <- data.frame(scale(dat1))
dat1_est$T_Height <- dat1$T_Height
reg       <- lm(T_Height~., data=dat1_est)

plot(reg$coefficients, bcpa)

kable(cbind(reg$coefficients, bcpa),col.names = c('Nombre variable', 'beta_ord', 'bca'))

Nombre variable	beta_ord	bca
(Intercept)	48.1187151	48.1187151
Age	1.1423529	0.8430151
Head_Circ	0.8569511	0.7690213
Should_WD	0.7322031	0.6972929
Chest_WD	1.1269183	0.9569469
Belly	1.3488450	1.4553820
Waist	0.7875668	0.1926268
Hips	-0.2410071	0.3614580
Arm_LEN	1.2727307	1.2734246
Should_Waist	0.7099306	1.6272107
Waist_Knee	2.4470683	1.9976160
Leg_LEN	3.8472007	3.7548091

print(sum(reg$residuals^2/(n-k-1)))  #estimación de varianza común

#> [1] 55.99043

print(MSres)                         #estimación de la varianza resultante de reducir rango

#> [1] 56.5971

Proposición 4.1 (FIV como elementos en la diagonal de la inversa generalizada) Tenemos que \(\operatorname{FIV}_j= \left[\operatorname{diag}(W^{\top}W)^{-}_*\right]_j\).

Demostración

Tenemos \[ \left(W^{\top} W\right)^{-}_* = \sum_{j=1}^{k-2} \frac{1}{\lambda_j}\boldsymbol{v}_j \boldsymbol{v}_j^{\top} = [\boldsymbol{v}_1,\cdots,\boldsymbol{v}_{k-2}] \left( \begin{array}{ccc} 1/\lambda_1 & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & 1/\lambda_{k-2} \end{array} \right) \left[ \begin{array}{c} \boldsymbol{v}_1^{\top} \\ \vdots \\ \boldsymbol{v}_{k-2}^{\top} \end{array} \right] \] y por otro lado \[ \begin{align} W^{\top}W & = \; \sum_{j=1}^k \lambda_j\boldsymbol{v}_j \boldsymbol{v}_j^{\top} = \sum_{j=1}^{k-2} \lambda_j \boldsymbol{v}_j \boldsymbol{v}_j^{\top} + \sum_{j=k-1}^k \lambda_j \boldsymbol{v}_j \boldsymbol{v}_j^{\top} \\ & = \; [\boldsymbol{v}_1,\cdots,\boldsymbol{v}_{k-2}] \left( \begin{array}{ccc} \lambda_1 & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & \lambda_{k-2} \end{array} \right) \left[ \begin{array}{c} \boldsymbol{v}_1^{\top} \\ \vdots \\ \boldsymbol{v}_{k-2}^{\top} \end{array} \right] + \sum_{j=k-1}^k \lambda_j \boldsymbol{v}_j \boldsymbol{v}_j^{\top} \end{align} \]

entonces, debido a que los vectores propios son ortonormales, tenemos que \[ (W^{\top}W)^{-}_*(W^{\top}W) = [\boldsymbol{v}_1,\cdots,\boldsymbol{v}_{k-2}] \left[ \begin{array}{c} \boldsymbol{v}_1^{\top} \\ \vdots \\ \boldsymbol{v}_{k-2}^{\top} \end{array} \right] \] y además \((W^{\top}W)^{+}_*(W^{\top}W)(W^{\top}W)^{+}_*\) es igual a \[ \begin{align} (W^{\top}W)^{+}_*(W^{\top}W)(W^{\top}W)^{+}_*&= [\boldsymbol{v}_1,\cdots,\boldsymbol{v}_{k-2}] \left[ \begin{array}{c} \boldsymbol{v}_1^{\top} \\ \vdots \\ \boldsymbol{v}_{k-2}^{\top} \end{array} \right] [\boldsymbol{v}_1,\cdots,\boldsymbol{v}_{k-2}] \left( \begin{array}{ccc} 1/\lambda_1 & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & 1/\lambda_{k-2} \end{array} \right) \left[ \begin{array}{c} \boldsymbol{v}_1^{\top} \\ \vdots \\ \boldsymbol{v}_{k-2}^{\top} \end{array} \right]\\[15pt] &=[\boldsymbol{v}_1,\cdots,\boldsymbol{v}_{k-2}]\left( \begin{array}{ccc} 1/\lambda_1 & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & 1/\lambda_{k-2} \end{array} \right) \left[ \begin{array}{c} \boldsymbol{v}_1^{\top} \\ \vdots \\ \boldsymbol{v}_{k-2}^{\top} \end{array} \right] \end{align} \] que es igual a \((W^{\top}W)^{+}_*\).

Entonces \[ \begin{align} \text{Var}(\widehat{\boldsymbol{\beta}}_{(-)}) & = \; \text{Var}( \; (W^{\top}W)^{-}_*W^{\top}y \; ) \\ & = \; \sigma^2 (W^{\top}W)^{-}_*W^{\top}W(W^{\top}W)^{-}_* = \sigma^2 (W^{\top}W)^{-}_* \end{align} \] así los factores de inflación de varianza estan dados por los elementos diagonales de \((W^{\top}W)^{-}_*\).

Puede verificarse que
\[ \text{Var}(\widehat{\boldsymbol{\beta}}) = \text{Var}(\widehat{\boldsymbol{\beta}}_{(-)}) + \sum_{j=k-1}^k \frac{1}{\lambda_j} \boldsymbol{v}_j \boldsymbol{v}_j^{\top}, \] así que \(\text{Var}(\widehat{\boldsymbol{\beta}}) \ge \text{Var}(\widehat{\boldsymbol{\beta}}_{(-)})\) (orden de Loewner).

Por otro lado, al igual que algunos métodos de regularización, se obtiene un estimador sesgado de \(\boldsymbol{\beta}\): \[ \begin{align} \mathbb{E}\left[\widehat{\boldsymbol{\beta}}_{(-)}\right] & = \; \mathbb{E}\left[\; (W^{\top}W)^{-}_*W^{\top}\boldsymbol{y} \; \right] = (W^{\top}W)^{-}_*W^{\top}\mathbb{E}[\boldsymbol{y}]\\ & = \; (W^{\top}W)^{-}_*W^{\top}(\alpha \boldsymbol{1}_n + W \boldsymbol{\beta}) = (W^{\top}W)^{-}_*W^{\top}W \boldsymbol{\beta} \\ & = \; [\boldsymbol{v}_1,\cdots,\boldsymbol{v}_{k-2}] \left[ \begin{array}{c} \boldsymbol{v}_1^{\top} \\ \vdots \\ \boldsymbol{v}_{k-2}^{\top} \end{array} \right] \boldsymbol{\beta} = (I - \boldsymbol{v}_{k-1}\boldsymbol{v}_{k-1}^{\top} - \boldsymbol{v}_k \boldsymbol{v}_k^{\top}) \boldsymbol{\beta} \\ & = \; \boldsymbol{\beta} - (\boldsymbol{v}_{k-1}\boldsymbol{v}_{k-1}^{\top} + \boldsymbol{v}_k \boldsymbol{v}_k^{\top}) \boldsymbol{\beta} = \boldsymbol{\beta} + \text{sesgo} \end{align} \]

4.2.2 Regresión de Componentes Principales

La regresión por componentes principales aborda el problema de la colinealidad desde la perspectiva de eliminar, del análisis, aquellas dimensiones de \(\mathcal{C}(X)\) que están causando el problema de colinealidad. Esto es similar, en concepto, a eliminar una variable independiente del modelo cuando no hay suficiente dispersión en dicha variable para aportar información significativa sobre \(\boldsymbol{Y}\). Sin embargo, en la regresión por componentes principales, la dimensión eliminada del análisis está definida por una combinación lineal de las variables, en lugar de por una sola variable independiente.

Este método consiste en ajustar un modelo de regresión lineal por mínimos cuadrados empleando como predictores las componentes generadas a partir de un anålisis de componentes principales (PCA). El objetivo es poder utilizar los principales componentes principales para reducir la dimensión aplicado a un problema de regresión.

Consideramos \(X\) la matriz de diseño sin la columna correspondiente a la ordenada al origen y la descomposición espectral \[ X^{\top} X =T \Lambda T^{\top}. \] Entonces \[ \begin{aligned} \mathbb{E}(\boldsymbol{Y}) & =\alpha \boldsymbol{1}_n+X \boldsymbol{\beta} \\ & =\alpha \boldsymbol{1}_n+X\, T\, T^{\top} \boldsymbol{\beta} \quad\left(\text { porque } T\, T^{\top}=I_p\right) \\ & =\alpha \boldsymbol{1}_n+(X T)\left(T^{\top} \boldsymbol{\beta}\right) \\ & =\alpha \boldsymbol{1}_n+Z \boldsymbol{\gamma}, \end{aligned} \] donde \(Z=X\,T=(T^{\top}X^{\top})^{\top}\) son los datos transformados según la rotación. En las columnas tiene los componentes principales. Por otro lado, \(\boldsymbol{\gamma}=T^{\top} \boldsymbol{\beta}\) es el vector de parámetros desconocidos. La estimación por mínimos cuadrados de \(\boldsymbol{\gamma}\) es:

\[ \begin{aligned} \widehat{\boldsymbol{\gamma}} & =\left(Z^{\top} Z\right)^{-1} Z^{\top} \mathbf{Y} \\ & =\Lambda^{-1} Z^{\top} \mathbf{Y} \end{aligned} \]

y el estimador de \(\boldsymbol{\beta}\) es: \[ \widehat{\boldsymbol{\beta}}=T \, \widehat{\boldsymbol{\gamma}}. \]

Los estimadores de regresión por componentes principales para \(\boldsymbol{\gamma}\) se obtienen eliminando \((p-r)\) columnas de \(Z\) para formar \(Z_{(r)}\), y luego utilizando mínimos cuadrados para obtener un estimador \(\widetilde{\boldsymbol{\gamma}}\) para \(\boldsymbol{\gamma}\) de: \[ \widetilde{\boldsymbol{\gamma}}=\left(Z_{(r)}^{\top}\, Z_{(r)}\right)^{-1} Z_{(r)}^{\top} \boldsymbol{Y}. \] El estimador por componentes principales de \(\boldsymbol{\beta}\) es entonces: \[ \begin{aligned} \widetilde{\boldsymbol{\beta}} & = T_{(r)} \boldsymbol{\gamma} \\ & =T_{(r)}\left(Z_{(r)}^{\top}\, Z_{(r)}\right)^{-1} Z_{(r)}^{\top} \boldsymbol{Y} \end{aligned} \tag{4.2}\]

donde \(T_{(r)}\) se forma eliminando las columnas correspondientes de \(T\), de modo que \(Z_{(r)}=X\, T_{(r)}\). Es decir, \(Z_{(r)}\) se formó eliminando las últimas \((p-r)\) columnas de \(Z\). Entonces, \[ \begin{aligned} \operatorname{Var}(\widetilde{\boldsymbol{\beta}}) & =\sigma^2 \operatorname{tr}\left\{\left(Z_{(r)}^{\top}\, Z_{(r)}\right)^{-1}\right\} \\ & =\sigma^2 \sum_{i=1}^r \frac{1}{\lambda_i}. \end{aligned} \]

Claramente, la mayor reducción en la varianza del estimador por componentes principales respecto a mínimos cuadrados se logra eliminando los \((p-r)\) componentes principales correspondientes a los valores propios más pequeños de \(X^{\top}X\).

La regresión por componentes principales puede considerarse como el uso de un estimador de mínimos cuadrados restringidos (Johnson et al., 1973). La ecuación Ecuación 4.2 es el estimador de mínimos cuadrados restringidos, sujeto a las \((p-r)\) restricciones: \[ T_{(-r)}\, \boldsymbol{\beta}=0, \] donde \(T_{(-r)}=\left(T_{r+1}, T_{r+2}, \ldots, T_k\right)\). Garnham (1979) muestró que el sesgo cuadrado debido al uso del estimador por componentes principales es: \[ \sum_{i=r+1}^p \gamma_i^2 \]

y, entonces el error cuadrático medio es: \[ \operatorname{MSE}(\widetilde{\boldsymbol{\boldsymbol{\beta}}})=\sigma^2 \sum_{i=1}^r \frac{1}{\lambda_i}+\sum_{i=r+1}^k \gamma_i^2. \]

Si consideraramos únicamente la varianza podriamos eliminar componentes si importar o no (Massy, 1965). En tal caso, el error cuadrático medio de \(\widehat{\boldsymbol{\beta}}\) puede ser grande y las predicciones de la respuesta poco precisas. Por lo tanto, es necesario considerar la correlación entre un componente principal y la respuesta.

Ejemplo

cov_escal<-data.frame(scale(dat1[,-12]))
kable(head(cov_escal))

Age	Head_Circ	Should_WD	Chest_WD	Belly	Waist	Hips	Arm_LEN	Should_Waist	Waist_Knee	Leg_LEN
1.2389990	0.3809359	0.7733769	1.0217364	-0.2170850	-0.6003351	0.2997646	0.5920035	1.3197154	1.6216904	-0.6098670
1.0699588	-0.4204846	1.6133282	0.4572586	-0.2170850	0.1977762	0.6428091	1.7079521	0.9479207	1.6216904	-0.8622015
0.9854386	0.1137957	0.7733769	0.2690993	-0.6118351	-1.0563988	-0.5006724	0.4060120	0.0184340	-0.4919011	-1.1145360
1.1544789	-0.1533444	1.1933526	0.6454178	-0.9078977	-0.0302556	-0.6150205	0.9639863	0.5761260	0.6609670	-0.7360343
1.0699588	-1.2219050	-0.0665744	0.6454178	-0.7105226	-0.9423828	1.2145498	1.1499778	0.7620234	2.9667032	-1.7453721
0.5628380	-0.9547649	0.9833648	0.6454178	-0.6118351	-0.3723033	-0.1576280	0.2200206	1.1338181	0.8531117	-0.9883687

medidasCuerpo.pca <- prcomp(W, center=TRUE, scale.=TRUE)
summary(medidasCuerpo.pca)

#> Importance of components:
#>                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
#> Standard deviation     2.1697 1.1765 1.00513 0.88059 0.79918 0.77046 0.71147
#> Proportion of Variance 0.4279 0.1258 0.09184 0.07049 0.05806 0.05396 0.04602
#> Cumulative Proportion  0.4279 0.5538 0.64563 0.71612 0.77418 0.82815 0.87416
#>                            PC8     PC9    PC10    PC11
#> Standard deviation     0.66047 0.58822 0.56350 0.53332
#> Proportion of Variance 0.03966 0.03145 0.02887 0.02586
#> Cumulative Proportion  0.91382 0.94528 0.97414 1.00000

res1 <- cor(medidasCuerpo.pca$x, method="pearson")
corrplot::corrplot(res1, method= "color", order = "hclust", tl.pos = 'n')

cumsum(medidasCuerpo.pca$sdev^2)/sum(medidasCuerpo.pca$sdev^2)

#>  [1] 0.4279445 0.5537823 0.6456263 0.7161206 0.7741836 0.8281475 0.8741645
#>  [8] 0.9138210 0.9452753 0.9741422 1.0000000

pcs <- as.data.frame(medidasCuerpo.pca$x)
pcr.data <- cbind(T_Height=dat1[,12], pcs)
omr.data <- cbind(T_Height=dat1[,12], cov_escal)

lmodel <- lm(T_Height ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 , data = pcr.data)
omodel <- lm(T_Height ~ ., data = omr.data)

summary(omodel)

#> 
#> Call:
#> lm(formula = T_Height ~ ., data = omr.data)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -40.517  -3.107   0.739   3.543  23.563 
#> 
#> Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   48.1187     0.2796 172.073  < 2e-16 ***
#> Age            1.1424     0.3901   2.928 0.003517 ** 
#> Head_Circ      0.8570     0.3176   2.698 0.007132 ** 
#> Should_WD      0.7322     0.3449   2.123 0.034085 *  
#> Chest_WD       1.1269     0.3473   3.245 0.001231 ** 
#> Belly          1.3488     0.3255   4.144 3.82e-05 ***
#> Waist          0.7876     0.4026   1.956 0.050855 .  
#> Hips          -0.2410     0.4295  -0.561 0.574877    
#> Arm_LEN        1.2727     0.3526   3.609 0.000329 ***
#> Should_Waist   0.7099     0.4414   1.608 0.108241    
#> Waist_Knee     2.4471     0.4362   5.610 2.91e-08 ***
#> Leg_LEN        3.8472     0.3663  10.504  < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 7.483 on 704 degrees of freedom
#> Multiple R-squared:  0.627,  Adjusted R-squared:  0.6211 
#> F-statistic: 107.6 on 11 and 704 DF,  p-value: < 2.2e-16

summary(lmodel)

#> 
#> Call:
#> lm(formula = T_Height ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7, 
#>     data = pcr.data)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -38.687  -3.256   0.800   3.728  23.748 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  48.1187     0.2804 171.611  < 2e-16 ***
#> PC1          -4.2534     0.1293 -32.889  < 2e-16 ***
#> PC2           0.5547     0.2385   2.326  0.02030 *  
#> PC3           1.5282     0.2792   5.474 6.10e-08 ***
#> PC4           0.9599     0.3186   3.012  0.00268 ** 
#> PC5           1.8357     0.3511   5.229 2.25e-07 ***
#> PC6           0.9185     0.3642   2.522  0.01188 *  
#> PC7          -1.2009     0.3944  -3.045  0.00241 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 7.503 on 708 degrees of freedom
#> Multiple R-squared:  0.6228, Adjusted R-squared:  0.6191 
#> F-statistic:   167 on 7 and 708 DF,  p-value: < 2.2e-16

Una crítica a la regresión por componentes principales es que la eliminación de algunos componentes no provocará, en general, la eliminación de ninguna de las variables originales, por lo que la ecuación final es tan compleja como la que se forma utilizando mínimos cuadrados ordinarios.

Otra desventaja de esta metodología para resolver el problema de multicolinealidad, es la disminución/pérdida de interpretabilidad. Los componentes principales obtenidos a través del ACP son combinaciones lineales de las variables originales, y su interpretación en términos de las covariables originales puede no ser directa. En consecuencia, los componentes transformados pueden carecer de un significado directo o relevancia teórica.

Estimación del modelo de regresión vía Componentes Principales es similar a otra técnica, PLS, (Partial Least Squares), en el sentido que hacen la regresión sobre un subespacio del espacio de columnas de \(X\).

Otros métodos se presentan en Jolliffe (1972).

4.3 Transformaciones de variables

La transformación de variables es una técnica utilizada en estadística y análisis de datos para modificar una variable con el objetivo de cumplir con ciertos supuestos o mejorar la interpretación de los datos. Esta transformación implica aplicar una función matemática a los valores de la variable original para obtener una nueva variable.

La transformación de variables implica aplicar una función matemática a los valores de una variable con el objetivo de modificar su distribución o relación con otras variables.

Algunas transformaciones comunes incluyen la transformación logarítmica, la transformación exponencial y la transformación de raíz cuadrada. Estas transformaciones se utilizan principalmente cuando los datos presentan asimetría o heterocedasticidad.

Por ejemplo, si se tiene una variable con una distribución sesgada hacia la derecha, se puede aplicar una transformación logarítmica para reducir la asimetría y hacer que los datos se aproximen más a una distribución normal (gráfica o distribución se asemeja a la campana de Gauss).

Existen diferentes tipos de transformaciones de variables que se utilizan según las características de los datos y los objetivos del análisis. Algunas de las transformaciones más comunes son:

• Logarítmica,
• Recíprocos,
• Raíz cuadrada,
• Arcoseno
• Transformaciones estabilizadoras de varianza
• Transformaciones Box-Cox

Ejemplo: Transformación a homocedasticidad.

Si se tiene idea de la relación funcional entre la media y la varianza de una variable aleatoria \(Y\), es posible encontar una transformación \(h(Y)\) tal que su varianza no dependa de la media:

• Supongamos \(\sigma^2 = \Omega(\mu)\). Queremos \(h(Y)\) tal que \(\operatorname{Var}(h(Y)) \equiv c\), \(c\) constante.
• Usando una aproximación de primer orden \(h(Y) \approx h(\mu)+h'(\mu)(Y-\mu)\), \[ \begin{align} \operatorname{Var}(h(Y)) &\approx [h'(\mu)]^2\operatorname{Var}(Y)\\ &=[h'(\mu)]^2\ \Omega\,(\mu). \end{align} \] Si queremos que la varianza de la variable transformada no dependa de \(\mu\), entonces \(h\) debe ser tal que \[ h'(\mu) = \left[\frac{c}{\Omega\,(\mu)}\right]^{1/2}. \]

Equivalentemente, \[ h(\mu) = \int \frac{d\mu}{[\Omega(\mu)]^{1/2}}. \]

Por ejemplo, si \(\sigma^2\) es proporcional a la media \(\mu\), como es el caso de datos de conteos Poisson, entonces \(h(\mu) = \int \frac{d\mu}{{a\mu}^{1/2}} = c \sqrt{\mu}\). Entonces la transformación raíz cuadrada estabiliza la varianza.

par(mfrow=c(2,1),oma=c(0, 0, 0, 0),mar=c(2,2,1,1))
nobs<-500
mu <- seq(1,100,length=nobs)
dat <- rpois(n=nobs,mu)
plot(mu,dat,pch=19)

plot(mu,sqrt(dat),pch=19)

4.4 Heterogeneidad de varianza y correlación

4.4.1 Transformaciones Box-Cox

En 1964, G. Box y D. Cox estudiaron la familia de transformaciones potencia \[ y^{(\lambda)} = \left\{ \begin{array}{ll} \frac{y^{\lambda}-1}{\lambda} & \text{si } \lambda \ne 0\\ \text{log}(y) & \text{si } \lambda = 0 \end{array} \right. \]

• Esta familia contiene las transformaciones, muy usadas, raíz cuadrada, log, e inversa.
• Está reescalada de forma que sea continua en \(\lambda=0\).
• Es particularmente útil en casos de variables con colas largas.

Consideremos los datos de datos \((\boldsymbol{x}_1,y_1),\ldots,(\boldsymbol{x}_n,y_n)\) quen en notación de vectores y matrices es \(X\) y \(\boldsymbol{y}\).

Suponemos que, para cada \(i\), \(y_i^{(\lambda)}\) es una función monótona de \(y_i\) y que, para alguna \(\lambda\) (desconocida) el vector de respuestas transformadas \(\boldsymbol{y}^{(\lambda)}=(y_1^{(\lambda)},\ldots, y_n^{(\lambda)})^{\top}\) sigue un modelo lineal \[ \boldsymbol{y}^{(\lambda)} = X\boldsymbol{\beta} +\boldsymbol{\epsilon}, \qquad \boldsymbol{\epsilon} \sim N_n(0, \ \sigma^2 I). \] Como la verosimilitud es la conjunta de los observables, ¿Cuál es la distribución de los datos observables \(\boldsymbol{y}\)?

Del resultado de cambio de variable: \[ u \sim f_U(u), \qquad v= h(u), \qquad f_V(v) = f_U\left(h^{-1}(v)\right) \; \left| \frac{d h^{-1}(v)}{dv} \right|, \] tenemos, \[ \begin{align} y^{(\lambda)} \sim N(\boldsymbol{y}^{(\lambda)}) , \qquad y = g^{-1}(\boldsymbol{y}^{(\lambda)}) , \qquad f_Y(\boldsymbol{y}) = N(g(\boldsymbol{y})) \left| \frac{d g(\boldsymbol{y})}{d\boldsymbol{y}} \right|\\[10pt] \text{con }\quad g(\boldsymbol{y})=\left(\frac{y_1^{\lambda}-1}{\lambda},\cdots, \frac{y_n^{\lambda}-1}{\lambda}\right). \end{align} \] En nuestro caso, el valor absoluto del Jacobiano es \[ \left| \frac{d g(\boldsymbol{y})}{d\boldsymbol{y}} \right| = \left|\frac{d \boldsymbol{y}^{(\lambda)})}{d\boldsymbol{y}} \right| = \text{abs}\left| \begin{array}{cccc} y_1^{\lambda-1} & & & \\ & y_2^{\lambda-1} & & \\ & & \ddots & \\ & & & y_n^{\lambda-1} \end{array} \right| = \prod_{i=1}^n |y_i^{\lambda - 1}| \] Nota: Se pueden usar diferentes \(\lambda\)’s para las diferentes \(y_i\)’s

La verosimilitud es entonces \[ f(\boldsymbol{y}) = \frac{1}{(2\pi \sigma^2)^{n/2}} \text{exp}\left\{ -\frac{1}{2 \sigma^2} (\boldsymbol{y}^{(\lambda)}-X\boldsymbol{\beta})^{\top}(\boldsymbol{y}^{(\lambda)}-X\boldsymbol{\beta}) \right\} \; \prod_{i=1}^n |y_i^{\lambda - 1}|. \] A partir de esta, podemos obtener la verosimilitud perfil de \(\lambda\). Primero derivando la logverosimilitud e igualando a cero. Así, \[ \begin{align} \widehat{\beta}(\lambda) = & \; (X^{\top}X)^{-1}X^{\top}\boldsymbol{y}^{(\lambda)}\\ \widehat{\sigma}^2 (\lambda) = & \; \frac{1}{n} SS_{\mathrm{Res}}(\lambda,\widehat{\beta}) \end{align} \] donde \[ SS_{\mathrm{Res}}(\lambda,\beta) = (y^{(\lambda)}-X\boldsymbol{\beta})^{\top}(y^{(\lambda)}-X\boldsymbol{\beta}). \]

Perfil de \(\lambda\)

Sustituyendo los estimadores en la logverosimilitud completa, obtenemos la perfil \[ \ell(\lambda) = -\frac{n}{2} \text{log} \frac{SS_{\mathrm{Res}}(\lambda,\widehat{\beta}(\lambda))}{n} + (\lambda -1) \sum_{i=1}^n \text{log}|y_i| \]

Un intervalo de aproximadamente \((1-\alpha)\times 100\%\) de confianza es el conjunto de \(\lambda^*\)’s tales que \[ \left\{ \lambda^* \left| -2\text{log}\frac{L(\lambda^*)}{L(\widehat{\lambda})} \le \chi^2(\alpha;\ v) \right. \right\} \]

Ejemplo: Experimento textil

El siguiente ejemplo aparee en el artículo original Box & Cox (1964) “A Textile Experiment using a Single Replicate of a 3a Design” que usa datos de un reporte técnico de la organización internacional de lana textil. Los datos reportan el número de ciclos para falla (\(y\)) obtenidos de experimentos con factores:

• \(x_1\): Longitud de la muestra (250, 300, 350 mm)
• \(x_2\): Amplitud del ciclo de carga (8, 9, 10 mm)
• \(x_3\): Carga (40, 45, 50 gm)

Los datos en tabla representan los niveles de \(x\) denotadas convencionalmente como \(-1, 0, 1\). La figura presenta los primeros renglones de la tabla original en el artículo

Figura 4.1: Datos ejemplo en Box & Cox (1964)

# Datos
#    originales

x1     <- c(250,300,350)
x2     <- c(8,9,10)
x3     <- c(40,45,50)
diseno <- as.matrix(expand.grid(Var3 = x1, Var2 = x2, Var1 = x3))[,c(3,2,1)]

#    codificados
a      <- c(-1,0,1)
diseno <- as.matrix(expand.grid(Var3 = a ,Var2 = a ,Var1 = a ))[,c(3,2,1)]
ciclos <- c( 674, 370, 292, 338, 266, 210, 170, 118,  90,
            1414,1198, 634,1022, 620, 438, 442, 332, 220,
            3636,3184,2000,1568,1070, 566,1140, 884, 360)

y      <- ciclos
n      <- length(y)
X      <- cbind(rep(1,n),diseno)

par(mfrow = c(2,2),mar = c(4, 4, 2, 2))
plot(X[,2], y, pch = 20,xlab = "x1", cex.axis = .9, xlim = c(-1.1,1.1))
plot(X[,3], y, pch = 20,xlab = "x2", cex.axis = .9, xlim = c(-1.1,1.1))
plot(X[,4], y, pch = 20,xlab = "x3", cex.axis = .9, xlim = c(-1.1,1.1))


# Modelo lineal
out <- lm(y ~ -1 + X)
summary(out)

#> 
#> Call:
#> lm(formula = y ~ -1 + X)
#> 
#> Residuals:
#>    Min     1Q Median     3Q    Max 
#> -644.6 -279.1 -150.1  199.6 1268.0 
#> 
#> Coefficients:
#>       Estimate Std. Error t value Pr(>|t|)    
#> X       861.33      93.95   9.168 3.83e-09 ***
#> XVar1   660.00     115.06   5.736 7.66e-06 ***
#> XVar2  -535.89     115.06  -4.657 0.000109 ***
#> XVar3  -310.78     115.06  -2.701 0.012751 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 488.2 on 23 degrees of freedom
#> Multiple R-squared:  0.8639, Adjusted R-squared:  0.8402 
#> F-statistic: 36.49 on 4 and 23 DF,  p-value: 1.201e-09

par(mfrow=c(3,2),mar=c(5, 5, 2, 2))

plot(out, which=1:6, caption=c("Residuales vs Ajustados", "QQ Normal", "Escala-Localizacion","D de Cook", "Residuales vs Influencia", 
                               "D de Cook vs Influencia"))

# Función transformación Box-Cox
yt  <- function(lam){  
  if(lam!=0){
    return((-1+y^lam)/lam)
  }else{
    return(log(y))
  }
}

# Función verosimilitud Perfil lambda
logperfil <- function(lam){
  ylam <- yt(lam) 
  bet  <- solve( t(X) %*% X , t(X) %*% ylam )
  sig2 <- mean((ylam - X %*% bet)^2)
  return(cte + (lam - 1) * sum(log(y)) - n * log(sig2)/2) 
}
 

cte  <- -n * (1 + log(2*pi))/2
m    <- 201
lamb <- seq(-1, 1, length=m)
perf <- rep(0, m)

for(i in 1:m){
  perf[i] <- logperfil(lamb[i])
}


# sobre el EMV de lambda e intervalos
aa     <- optimize(logperfil,interval=c(-1,1),maximum=TRUE)
hatlambda <- aa$maximum                        # [1] -0.05927981
hatlambda

#> [1] -0.05927981

linf   <- aa$objective-qchisq(.95,1)/2
rango  <- lamb[(abs(perf-linf)< .15)]

# Grafica de la log verosimilitud perfil de lambda
par(mfrow=c(1,1))
plot(lamb,perf,xlab=expression(lambda),type="l",col="blue",lwd=2,
 ylab=expression(paste("Perfil( ",lambda," )")),mgp=c(1.5,.5,0))
segments(hatlambda, min(perf), hatlambda, aa$objective, col = gray(.7))
segments(rango[1], linf, rango[2], linf, col = "red", lwd = 1.5)
segments(rango[1], min(perf), rango[1], linf, col = "red", lwd = 1.5)
segments(rango[2], min(perf), rango[2], linf, col = "red", lwd = 1.5)

Ahora transformamos los datos con el valor de lambda que maximiza la verosimilitud perfil

ytran  <- (-1 + y^hatlambda)/hatlambda

par(mfrow=c(2,2),mar = c(4, 4, 2, 2))
plot(X[,2], ytran, pch = 20, xlab = "x1", cex.axis = .9, xlim = c(-1.1,1.1))
plot(X[,3], ytran, pch = 20, xlab = "x2", cex.axis = .9, xlim = c(-1.1,1.1))
plot(X[,4], ytran, pch = 20, xlab = "x3", cex.axis = .9, xlim = c(-1.1,1.1))

Verificamos el ajuste del modelo usando los valores transformados

out2 <- lm(ytran ~ -1+X)
summary(out2)

#> 
#> Call:
#> lm(formula = ytran ~ -1 + X)
#> 
#> Residuals:
#>       Min        1Q    Median        3Q       Max 
#> -0.282038 -0.085499  0.008526  0.087029  0.183945 
#> 
#> Coefficients:
#>       Estimate Std. Error t value Pr(>|t|)    
#> X      5.26369    0.02412 218.202  < 2e-16 ***
#> XVar1  0.57160    0.02954  19.347 9.96e-16 ***
#> XVar2 -0.43318    0.02954 -14.662 3.68e-13 ***
#> XVar3 -0.26939    0.02954  -9.118 4.24e-09 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.1253 on 23 degrees of freedom
#> Multiple R-squared:  0.9995, Adjusted R-squared:  0.9994 
#> F-statistic: 1.207e+04 on 4 and 23 DF,  p-value: < 2.2e-16

par(mfrow = c(3,2), mar = c(5, 5, 2, 2))
plot(out2, which=1:6, caption=c("Residuales vs Ajustados", "QQ Normal", "Escala-Localizacion","D de Cook", "Residuales vs Influencia", 
                                "D de Cook vs Influencia"))

Alternativamente

library(MASS)
boxcox(out, plotit = T)

boxcox(out, plotit = T,lambda = seq(-0.5,0.5,by = 0.05))

4.4.2 Mínimo Cuadrados Generalizados

En el marco del modelo clásico, los supuestos de homocedasticidad, \(E\left(\epsilon_i^2\right)=\sigma^2\ \ (i=1,2, \ldots n)\), y ausencia de autocorrelación, \(E\left(\epsilon_i \epsilon_j\right)=0\) \(\ \forall i \neq j \ (i, j=1,2, \ldots n)\), implican que la matriz de varianzas y covarianzas de \(\boldsymbol{\epsilon}\) es

\[ \begin{align} \operatorname{VarCov}(\boldsymbol{\epsilon}) & =\mathbb{E}\left[(\boldsymbol{\epsilon}-\mathbb{E}(\boldsymbol{\epsilon}))\left(\boldsymbol{\epsilon}-\mathbb{E}\left(\boldsymbol{\epsilon}\right)\right)^{\top}\right]=\mathbb{E}\left[\left(\begin{array}{c} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_n \end{array}\right)\left(\begin{array}{llll} \epsilon_1 & \epsilon_2 & \ldots & \epsilon_n \end{array}\right)\right] \\ & =\left(\begin{array}{cccc} \mathbb{E}\left(\epsilon_1^2\right) & \mathbb{E}\left(\epsilon_1 \epsilon_2\right) & \ldots & \mathbb{E}\left(\epsilon_1 \epsilon_n\right) \\ \mathbb{E}\left(\epsilon_2 \epsilon_1\right) & \mathbb{E}\left(\epsilon_2^2\right) & \ldots & \mathbb{E}\left(\epsilon_2 \epsilon_n\right) \\ \vdots & \vdots & \ddots & \vdots \\ \mathbb{E}\left(\epsilon_n \epsilon_1\right) & \mathbb{E}\left(\epsilon_n \epsilon_2\right) & \ldots & \mathbb{E}\left(\epsilon_n^2\right) \end{array}\right)=\left(\begin{array}{cccc} \sigma^2 & 0 & \ldots & 0 \\ 0 & \sigma^2 & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & \sigma^2 \end{array}\right)=\sigma^2 I_n \end{align} \]

El supuesto de homocedasticidad implica que todos los elementos de la diagonal principal de \(\operatorname{VarCov}(\boldsymbol{\epsilon})\), las varianzas, son iguales a un escalar \(\sigma^2\), mientras que el de autocorrelación implica que los elementos situados fuera de la diagonal principal de \(\operatorname{VarCov}(\boldsymbol{\epsilon})\), las covarianzas, son iguales a cero. Si relajamos estos dos supuestos, entonces la matriz de varianzas y covarianzas deja de ser escalar

\[ \operatorname{VarCov}(\boldsymbol{u})=\left(\begin{array}{cccc} \mathbb{E}\left(\epsilon_1^2\right) & \mathbb{E}\left(\epsilon_1 \epsilon_2\right) & \ldots & \mathbb{E}\left(\epsilon_1 \epsilon_n\right) \\ \mathbb{E}\left(\epsilon_2 \epsilon_1\right) & \mathbb{E}\left(\epsilon_2^2\right) & \ldots & \mathbb{E}\left(\epsilon_1 \epsilon_n\right) \\ \vdots & \vdots & \ddots & \vdots \\ \mathbb{E}\left(\epsilon_n \epsilon_1\right) & \mathbb{E}\left(\epsilon_n \epsilon_2\right) & \ldots & \mathbb{E}\left(\epsilon_n^2\right) \end{array}\right)=\left(\begin{array}{cccc} \sigma_1^2 & \sigma_{12} & \ldots & \sigma_{1 n} \\ \sigma_{21} & \sigma_2^2 & \ldots & \sigma_{2 n} \\ \vdots & \vdots & \ddots & \vdots \\ \sigma_{n 1} & \sigma_{n 2} & \ldots & \sigma_n^2 \end{array}\right)=\Sigma \]

Definición 4.2 (El modelo lineal general con perturbaciones no esféricas) El modelo lineal general con perturbaciones no esféricas es \[ Y_i=\beta_1+\beta_2 X_{2 i}+\cdots+\beta_k X_{k i}+\epsilon_i, \quad i=1, \ldots, n \] en donde \(\mathbb{E}\left(\epsilon_i\right)=0\), \(\mathbb{E}\left(\epsilon_i^2\right)=\sigma_i^2\)  y  \(\mathbb{E}\left(\epsilon_i \epsilon_j\right)=\sigma_{i j}\)   \(\forall i \neq j\), que podemos escribir en notación matricial como \[ \boldsymbol{Y}=X \boldsymbol{\beta}+\boldsymbol{u},\quad \text{ con } \quad \mathbb{E}(\boldsymbol{\epsilon})=\boldsymbol{0}, \quad \mathbb{E}\left(\boldsymbol{\epsilon \epsilon}^{\top}\right) = \Sigma. \]

Conviene escribir \(\Sigma=\sigma^2 \Omega\), para obtener el modelo lineal general con perturbaciones esféricas como un caso especial del modelo con perturbaciones no esféricas, \(\Omega = I_n\).

A continuación vamos a demostrar que en el modelo lineal general con perturbaciones no esféricas, hay un estimador alternativo y superior al estimador de mínimos cuadrados ordinarios:   el estimador de mínimos cuadrados generalizados, que es equivalente al estimador de máxima verosimilitud cuando el vector de errores sigue una distribución normal multivariada con matriz de Var-Cov \(\Sigma\).

Estimador de mínimos cuadrados ordinarios (EMO)

Proposición 4.2 En el modelo lineal general con heterocedasticidad y/o autocorrelación, el estimador de mínimos cuadrados ordinarios es \[ \widehat{\boldsymbol{\beta}}_{\text{MCO}} = \left(X^{\top} X\right)^{-1} X^{\top} \boldsymbol{Y}. \]

Esto es trivial demostrar ya que el método de mínimos cuadrados ordinarios no tiene en cuenta la matriz de varianzas y covarianzas de los errores al minimizar la suma de cuadrados de los residuos: \[ Q=(\boldsymbol{Y}-X \widehat{\boldsymbol{\beta}})^{\top}(\boldsymbol{Y}-X \widehat{\boldsymbol{\beta}}). \]

Proposición 4.3 En el modelo lineal general con heterocedasticidad y/o autocorrelación, el estimador \(\text{MCO}\) es insesgado.

Esta proposición es lógica porque que la propiedad de insesgadez se basa en los supuestos de regresores no estocásticos y \(\mathbb{E}(\boldsymbol{\epsilon})=\boldsymbol{0}\), pero no tiene en cuenta la matriz de varianzas y covarianzas de los errores. Así que \[ \mathbb{E}\left(\widehat{\boldsymbol{\beta}}_{\text{MCO}} \right) =\left[\left(X^{\top} X\right)^{-1} X^{\top} (X\boldsymbol{\beta})\right]=\boldsymbol{\beta}. \]

Proposición 4.4 En el modelo lineal general con heterocedasticidad y/o autocorrelación, la matriz de varianzas y covarianzas de \(\widehat{\boldsymbol{\beta}}_{\text{MCO}}\) es \[ \operatorname{Var}\left(\widehat{\boldsymbol{\beta}}_{M C O}\right)=\sigma^2\left(X^{\top} X\right)^{-1}\left(X^{\top} \Omega X\right)\left(X^{\top} X\right)^{-1} \]

Demostración

Por definición \[ \operatorname{VarCov}\left(\widehat{\boldsymbol{\beta}}_{\text{MCO}}\right)=\mathbb{E}\left[\left(\widehat{\boldsymbol{\beta}}_{\text{MCO}}-\mathbb{E}\left(\widehat{\boldsymbol{\beta}}_{\text{MCO}}\right)\right)\left(\widehat{\boldsymbol{\beta}}_{\text{MCO}}-\mathbb{E}\left(\widehat{\boldsymbol{\beta}}_{\text{MCO}}\right)\right)^{\top}\right]. \]

Como el estimador es insesgado, \[ \begin{align} \widehat{\boldsymbol{\beta}}_{\text{MCO}}-\mathbb{E}\left(\widehat{\boldsymbol{\beta}}_{\text{MCO}}\right)=\widehat{\boldsymbol{\beta}}_{\text{MCO}}-\boldsymbol{\beta} &= \left(X^{\top} X\right)^{-1} X^{\top}\boldsymbol{Y}-\boldsymbol{\beta}\\ &=\left(X^{\top} X\right)^{-1} X^{\top}(\boldsymbol{Y}-X\boldsymbol{\beta})\\ &=\left(X^{\top} X\right)^{-1} X^{\top} \boldsymbol{\epsilon}. \end{align} \] Entonces \[ \begin{aligned} \operatorname{VarCov}\left(\widehat{\boldsymbol{\beta}}_{\text{MCO}}\right) & =\mathbb{E}\left[\left(X^{\top} X\right)^{-1} X^{\top} \boldsymbol{\epsilon} \boldsymbol{\epsilon}^{\top} X\left(X^{\top} X\right)^{-1}\right]=\left(X^{\top} X\right)^{-1} X^{\top} \mathbb{E}\left[\boldsymbol{\epsilon} \boldsymbol{\epsilon}^{\top}\right] X\left(X^{\top} X\right)^{-1} \\ & =\left(X^{\top} X\right)^{-1} X^{\top}\left[\sigma^2\ \Omega\right] X\left(X^{\top} X\right)^{-1} =\sigma^2\left(X^{\top} X\right)^{-1}\left(X^{\top} \Omega X\right)\left(X^{\top} X\right)^{-1} \end{aligned} \]

Proposición 4.5 (Distribución del estimador) Si suponemos que \(\boldsymbol{\epsilon} \sim N\left(\mathbf{0}, \sigma\ \Omega\right)\), el estimador MCO tiene una distribución normal \[ \widehat{\boldsymbol{\beta}}_{\text{MCO}} \sim N\left(\boldsymbol{\beta},\ \ \sigma^2 \ \left(X^{\top} X\right)^{-1}\left(X^{\top} \Omega X\right)\left(X^{\top} X\right)^{-1}\right). \]

Proposición 4.6 El estimador \(\hat{\sigma}^2=\boldsymbol{e}^{\top} \boldsymbol{e} /(n-p)\) es un estimador sesgado.

Demostración

\[ \begin{align} \mathbb{E}\left(\boldsymbol{e}^{\top} \mathbf{e}\right)&=\mathbb{E}\left(\boldsymbol{Y}^{\top} (I_n-P) \boldsymbol{Y}\right)\\ &=\operatorname{tr}\left[(I_n-P)\Sigma\right]+(X\boldsymbol{\beta})^{\top}(I_n-P) (X\boldsymbol{\beta}),\quad \text{donde }\ \Sigma=\sigma^2\Omega\\ &=\sigma^2\operatorname{tr}\left[(I_n-P)\Omega\right]+0\\ &\neq \sigma^2(n-p), \qquad \text{ en general.} \end{align} \]

El estimador de mínimos cuadrados generalizados

Ahora planteamos la siguiente pregunta: ¿es posible transformar un modelo lineal general con perturbaciones no esféricas en un modelo lineal general con perturbaciones esféricas?

Si la respuesta es afirmativa, entonces el modelo transformado cumplirá las hipótesis básicas y todos los resultados establecidos en los temas anteriores serán de aplicación directa. El estimador de mínimos cuadrados ordinarios (MCO) en el modelo transformado se denomina estimador de mínimos cuadrados generalizados (MCG).

Para encontrar un modelo transformado con las hipótesis básicas, premultiplicamos el modelo lineal general con perturbaciones no esféricas por una matriz \(T\) no estocástica

\[ T\boldsymbol{Y}=TX \boldsymbol{\beta}+T\boldsymbol{\epsilon} \]

Este modelo transformado puede escribirse como \[ \boldsymbol{Y}_*=X_* \boldsymbol{\beta}+\boldsymbol{\epsilon}_* \]

El término de error en el modelo transformado \(\boldsymbol{\epsilon}_*\) cumple las siguientes propiedades:

1. \(\mathbb{E}\left(\boldsymbol{\epsilon}_*\right)=\mathbb{E}(T\boldsymbol{\epsilon})=T \mathbb{E}(\boldsymbol{\epsilon})=\boldsymbol{0}\)
2. \(\mathbb{E}\left(\boldsymbol{\epsilon}_* \boldsymbol{\epsilon}_*^{\top}\right)=\mathbb{E}\left(T\boldsymbol{\epsilon \epsilon}^{\top} T^{\top}\right)=T \mathbb{E}\left(\boldsymbol{\epsilon \epsilon}^{\top}\right) T^{\top} = \sigma^2\ T \Omega T^{\top}\).

Si la matriz \(T\) es tal que \(\sigma^2\, T\, \Omega\, T^{\top} = \sigma^2 I_n\), entonces el modelo transformado:

1. contiene los parámetros de interés \(\boldsymbol{\beta}\) y \(\sigma\)
2. cumple las hipótesis básicas.

De aquí, el estimador de mínimos cuadrados ordinarios en el modelo transformado proporciona el estimador lineal, insesgado y eficiente de \(\boldsymbol{\beta}\) y el estimador insesgado de \(\sigma^2\).

Proposición 4.7 (Existencia) Existe una matriz \(T\) tal que \(T\,\Omega\, T^{\top} = I_n\)

Demostración

Usando el teorema espectral, podemos escribir \[ \Omega \, C = C \, \Lambda \] en donde \(\Lambda=\operatorname{diag}\left(\lambda_1, \ldots, \lambda_n\right)\) es la matriz diagonal de eigenvalores y \(C\) es la matriz ortogonal con los eigenvectores. De aquí, podemos escribir \[ \Omega = C \Lambda C^{\top} \]

Como \(\Lambda^{1/2}=\operatorname{diag}\left(\sqrt{\lambda_1}, \ldots, \sqrt{\lambda_n}\right)\), tenemos que

\[ \Omega = C \Lambda^{1 / 2} \Lambda^{1 / 2} C^{\top} \]

Premultiplicando \(\Omega\) por \(\Lambda^{-1 / 2} C^{\top}\) y postmultiplicando por \(C \Lambda^{-1 / 2}\), \[ \Lambda^{-1 / 2} C^{\top} \Omega C \Lambda^{-1 / 2} = \Lambda^{-1 / 2} C^{\top} C \Lambda^{1 / 2} \Lambda^{1 / 2} C^{\top} C \Lambda^{-1 / 2} = I_n \] De aquí, vemos que la matriz buscada es \[ T = \Lambda^{-1 / 2} C^{\top}. \]

De la demostración anterior, se derivan las dos siguientes relaciones que son de interés:

1. \(\Omega^{-1}= T^{\top} T\) \[ \begin{align} \because \ \Omega^{-1}&=C\Lambda^{-1}C^{\top}\\ &=C\Lambda^{-1/2}\Lambda^{-1/2}C^{\top}\\ &=(\Lambda^{-1/2}C^{\top})^{\top}\Lambda^{-1/2}C^{\top} \end{align} \]
2. \(\Omega = T^{-1} \left(T^{\top}\right)^{-1}\) \[ \begin{align} \because \ T^{-1} \left(T^{\top}\right)^{-1}&=C^{\top}\Lambda^{1/2}\left(C\Lambda^{-1/2}\right)^{-1}\\ &=C^{\top}\Lambda^{1/2}(\Lambda^{1/2})C^{\top},\quad \text{porque }\ C^{-1}=C^{\top}. \end{align} \]

Proposición 4.8 (estimador de Mínimos Cuadrados Generalizados, MCG) El estimador lineal, insesgado y óptimo de \(\boldsymbol{\beta}\) es \[ \widehat{\boldsymbol{\beta}}_{\text{MCG}}=\left(X^{\top} \boldsymbol{\Omega}^{-1} X\right)^{-1} X^{\top} \Omega^{-1} \mathbf{Y} \] que se denomina estimador de Mínimos Cuadrados Generalizados o estimador de Aitken.

Demostración

Como el modelo transformado cumple los supuestos del modelo clásico, el estimador el estimador de mínimos cuadrados ordinarios \[ \widehat{\boldsymbol{\beta}}=\left(X_*^{\top} X_*\right)^{-1} X_*^{\top} \mathbf{Y}_* \] será el estimador lineal, insesgado y óptimo, que podemos expresarse en términos de los datos originales \[ \widehat{\boldsymbol{\beta}}=\left(X^{\top} T^{\top} T X\right)^{-1} X^{\top} T^{\top} T \, \mathbf{Y}=\left(X^{\top} \Omega^{-1} X\right)^{-1} X^{\top} \Omega^{-1} \mathbf{Y}. \]

Proposición 4.9 La matriz de varianzas y covarianzas del estimador de MCG es

\[ \operatorname{Var}\left(\widehat{\boldsymbol{\beta}}_{\text{MCG}}\right)=\sigma^2\left(X^{\top} \Omega^{-1} X\right)^{-1} \]

Demostración

La matriz de varianzas y covarianzas del estimador de MCO de \(\boldsymbol{\beta}\) en el modelo transformado es \[ \operatorname{VarCov}\left(\widehat{\boldsymbol{\beta}}_{\text{MCG}}\right)=\sigma^2\left(X_*^{\top} X_*\right)^{-1}=\sigma^2\left(X^{\top} T^{\top} T X\right)^{-1}=\sigma^2\left(X^{\top} \Omega^{-1} X\right)^{-1}. \]

El Teorema de Gauss-Markov aplicado en el modelo transformado, garantiza que el estimador MCO es el estimador con varianza m¶‡nima entre los estimadores lineales e insesgados.

Proposición 4.10 El estimador insesgado de \(\sigma^2\) es \[ \widehat{\sigma}_{\text{MCG}}^2=\frac{\left(\boldsymbol{Y}-X \widehat{\boldsymbol{\beta}}_{\text{MCG}}\right)^{\top} \Omega^{-1}\left(\boldsymbol{Y}-X \widehat{\boldsymbol{\beta}}_{\text{MCG}}\right)}{n-p} \]

Demostración

El estimador insesgado de \(\sigma^2\) en el modelo transformado es \[ \begin{align} \widehat{\sigma}^2 &=\frac{\left(\boldsymbol{Y}_*-X_* \widehat{\boldsymbol{\beta}}_{\text{MCG}}\right)^{\top}\left(\boldsymbol{Y}_*-X_* \widehat{\boldsymbol{\beta}}_{\text{MCG}}\right)}{n-p} =\frac{\left(\boldsymbol{Y}-X \widehat{\boldsymbol{\beta}}_{\text{MCG}}\right)^{\top} T^{\top}T\left(\boldsymbol{Y}-X \widehat{\boldsymbol{\beta}}_{\text{MCG}}\right)}{n-p}. \end{align} \]

Estimadores máximo verosímiles

Como \(\boldsymbol{Y}\sim N_n(X\boldsymbol{\beta}, \sigma^2\Omega)\), entonces \[ \ell\left(\boldsymbol{\beta}, \sigma^2\right)=-\frac{n}{2} \ln (2 \pi)-\frac{n}{2} \ln \left(\sigma^2\right)-\frac{1}{2} \ln (\Omega)-\frac{1}{2 \sigma^2}(\boldsymbol{Y}-X \boldsymbol{\beta})^{\top} \Omega^{-1}(\boldsymbol{Y}-X \boldsymbol{\beta}), \] ya que \(\left|\,\sigma^2 \Omega\,\right|=\left(\sigma^2\right)^n|\Omega|\).

Proposición 4.11 (EMVs) Los estimadores de máximo verosímiles de \(\boldsymbol{\beta}\) y \(\sigma_u^2\) son \[ \begin{aligned} \widehat{\boldsymbol{\beta}} & =\left(X^{\top} \Omega^{-1} X\right)^{-1} X^{\top} \Omega^{-1} \boldsymbol{Y} \\[5pt] \widehat{\sigma}^2 & =\frac{(\boldsymbol{Y}-X \widehat{\boldsymbol{\beta}})^{\top} \Omega^{-1}(\boldsymbol{Y}-X \widehat{\boldsymbol{\beta}})}{n} \end{aligned} \]

Bajo el supuesto normalidad, es estimador de mínimos cuadrados generalizados coincide con el estimador de máxima verosimilitud.

La demostración es trivial a partir de obtener las derivadas parciales de \(\ell\left(\boldsymbol{\beta}, \sigma^2\right)\) respecto a \(\boldsymbol{\beta}\) y \(\sigma^2\): \[ \begin{aligned} & \frac{\partial \ell\left(\boldsymbol{\beta}, \sigma^2\right)}{\partial \boldsymbol{\beta}}=-\frac{1}{\sigma^2}\left(-X^{\top} \Omega^{-1} \boldsymbol{Y}+X^{\top} \Omega^{-\mathbf{1}} X \boldsymbol{\beta}\right) \\ & \frac{\partial \ell\left(\boldsymbol{\beta}, \sigma^2\right)}{\partial \sigma^2}=-\frac{n}{2 \sigma^2}+\frac{1}{2 \sigma^4}(\boldsymbol{Y}-X \boldsymbol{\beta})^{\top}\Omega^{-1}(\boldsymbol{Y}-X \boldsymbol{\beta}) \end{aligned} \]

igualandolas a cero y resoviendo simultáneamente.

El estimador de \(\text{MCG}\) supone que la matriz \(\Omega\) es conocida.

Mínimo Cuadrados Generalizados Factibles

Hasta ahora hemos supuesto que conocemos \(\Sigma=\sigma^2 \Omega\) o al menos \(\Omega\). El estimador de MCG en este caso es lineal, insesgado y de varianza mínima. Pero en la práctica la mayoría de las veces \(\Omega\) o \(\Sigma\) son desconocidas. En este caso el estimador MCG no es directamente calculable. La solución habitual es sustituir \(\Omega\) (o \(\Sigma\) ) por una estimación suya en la expresión del estimador de MCG dando lugar al estimador MCGF: \[ \begin{aligned} \widehat{\boldsymbol{\beta}}_{\text{MCGF}} & =\left(X^{\top} \widehat{\Omega}^{-1} X\right)^{-1} X^{\top} \widehat{\Omega}^{-1} \boldsymbol{Y} \\ & =\left(X^{\top} \widehat{\Sigma}^{-1} X\right)^{-1} X^{\top} \widehat{\Sigma}^{-1} \boldsymbol{Y}. \end{aligned} \]

---

Bibliografía:

• Christensen, R. (2011). Plane answers to complex questions: The Theory of Linear Models. 4th Ed. New York: Springer.
• Seber, G. A., & Lee, A. J. (2003). Linear regression analysis. John Wiley & Sons.
• Radhakrishna R., Toutenburg, H.(1997) Linear Models Least Squares and Alternatives. Springer. (Cap 3.10)

Referencias:

• Box, G.E.P. & Cox, D.R. (1964) Analysis of Transformations. JRSS-B, Vol.26, No.2, 211-252.
• Thome, N. (2019). La inversa generalizada de Moore-Penrose y Aplicaciones. Publicaciones electrónicas Sociedad Matemática Mexicana, 21, 81. Liga a ebook
• Bott, R. y Duffin, R. J. (1953)On the algebra of networks, Trans. Amer. Math. Soc. 74:99-109.
• Drazin, M. P. (1958). “Pseudo-inverses in associative rings and semigroups”. The American Mathematical Monthly. 65 (7): 506–514.
• Rawlings, J. O., Pantula, S. G., & Dickey, D. A. (Eds.). (1998). Applied regression analysis: a research tool. New York, NY: Springer New York.
• Jolliffe, I. T. (1972). Discarding variables in a principal component analysis. I: Artificial data. Journal of the Royal Statistical Society Series C: Applied Statistics, 21(2), 160-173.