Introducción

En esta sección se explica únicamente la parte de la estandarización del código de preprocesamiento para un arreglo de tres vías.

La idea general es:

  1. partir de la matricización en la vía A, denotada por Xa,
  2. centrar los datos en la vía A,
  3. reconstruir las matricizaciones centradas en las vías B y C,
  4. calcular matrices diagonales de normalización para las vías B y C,
  5. obtener una matriz final estandarizada.

Para ilustrar el procedimiento, se usarán matrices pequeñas simuladas.

Datos simulados

Supongamos que se tienen dos ocasiones, cada una representada por una matriz de 3 individuos y 2 variables. Se tendrá por tanto un arreglo de tres vías \(\mathbf{\underline{X}} \in \mathbb{R}^{3 \times 2 \times 2}\). De esta manera, las matrices simuladas son:

\[\begin{equation} M_1 = \begin{pmatrix} 1 & 2 \\ 3 & 4 \\ 5 & 6 \end{pmatrix} \qquad M_2 = \begin{pmatrix} 10 & 20 \\ 30 & 40 \\ 50 & 60 \end{pmatrix} \end{equation}\]

M1 <- matrix(c(
  1, 2,
  3, 4,
  5, 6
), nrow = 3, byrow = TRUE)

M2 <- matrix(c(
  10, 20,
  30, 40,
  50, 60
), nrow = 3, byrow = TRUE)

rownames(M1) <- rownames(M2) <- c("A1", "A2", "A3")
colnames(M1) <- colnames(M2) <- c("V1", "V2")

lista <- list(T1 = M1, T2 = M2)

M1
##    V1 V2
## A1  1  2
## A2  3  4
## A3  5  6
M2
##    V1 V2
## A1 10 20
## A2 30 40
## A3 50 60

Funciones auxiliares

Se necesita reconstruir las matricizaciones a partir de \(Xa\).

Se va a crear una función llamada matricizar_Xa une las matrices por columnas. El resultado si se tiene \(K\) matrices de tamaño \(I\times J\) se obtiene una matriz \(X_a \in \mathbb{R}^{I \times (J\cdot K)}\)

matricizar_Xa <- function(lista_matrices) {
  do.call(cbind, lista_matrices)
}

Se va a crear una funnción llamada matricizar_Xb que aplicará la transpuesta, \(t\), a cada matriz y convierte cada matriz de tamaño \(I\times J\) en \(J\times I\), luego se une todas las matrices transpuestas por columnas. Como resultado se obtiene \(X_b \in \mathbb{R}^{J \times (I\cdot K)}\)

matricizar_Xb <- function(lista_matrices) {
  do.call(cbind, lapply(lista_matrices, t))
}

Se crea la función matriciar_Xc, la cual se convierte en vector columna, vectorizado por columnas \((vec(m))\), luego se uden todos los vectores como filas, de lo cual se obtiene \(X_c \in \mathbb{R}^{K \times (I\cdot J)}\)

matricizar_Xc <- function(lista_matrices) {
  do.call(rbind, lapply(lista_matrices, function(m) c(as.matrix(m))))
}

Se crea la función reconstruir_desde_Xa; pero ¿por qué existe está función?

Cuando se trabaja con datos de tres vías, el arreglo o tensor:

\[\begin{equation} \mathbf{\underline{X}} \in \mathbb{R}^{I \times J \times K} \end{equation}\]

pero para poder operar computacionalmente en R, se transforma en matrices:

\[\begin{equation} X_a \in \mathbb{R}^{I \times JK}, \; X_b \in \mathbb{R}^{J \times IK}, \; X_c \in \mathbb{R}^{K \times IJ} \end{equation}\]

Esto es la matricización.

Como el preprocesamiento se establece sobre \(X_a\), por ejemplo \(X_{ac}=PX_a\), donde \(P\) es una matriz de centrado. El problema clave aquí es que despues del centrado, \(X_{ac}\) ya no se tiene directamente el arreglor orifinal, sino solo la versión matricizada. Pero se necesita También: \(X_{bac}\) y \(X_{cac}\) para calcular las funciones \(W_j\) Y \(W_j\), como se vera más adelante. Entonces, la función hace exactamente esto:

\[\begin{equation} X_a \rightarrow \mathcal{X} \rightarrow X_b, X_c$ \end{equation}\]

Es decir, (matriz) \(\rightarrow\) (arreglo implícito) \(\rightarrow\) (otras matricizaciones)

Matematicamente se puede explicar de la siguiente manera. Sea un arreglo de tres vías:

\[\begin{equation} \mathcal{X} \in \mathbb{R}^{I \times J \times K} \end{equation}\]

y su matricización en la vía A:

\[\begin{equation} X_a \in \mathbb{R}^{I \times JK} \end{equation}\]

La función tiene como objetivo recuperar las otras matricizaciones (\(X_b\) y \(X_c\)) a partir de \(X_a\).

La matriz \(X_a\) puede escribirse como una concatenación de bloques:

\[\begin{equation} X_a = \left[ \underbrace{X^{(1)}}_{\text{ocasión 1}} \; \underbrace{X^{(2)}}_{\text{ocasión 2}} \; \cdots \; \underbrace{X^{(K)}}_{\text{ocasión K}} \right] \end{equation}\]

donde cada bloque:

\[\begin{equation} X^{(k)} \in \mathbb{R}^{I \times J}, \quad k = 1, \ldots, K \end{equation}\]

representa la matriz correspondiente a la \(k\)-ésima ocasión.

Cada bloque \(X^{(k)}\) corresponde a una ``rebanada’’ del arreglo original:

\[\begin{equation} \mathcal{X}_{::k} = X^{(k)} \end{equation}\]

De esta forma, la colección de matrices \(\{X^{(k)}\}_{k=1}^K\) permite reconstruir implícitamente el arreglo tridimensional \(\mathcal{X}\).

La matricización en la vía B se obtiene tomando la transpuesta de cada bloque y concatenándolas por columnas:

\[\begin{equation} X_b = \left[ (X^{(1)})^T \; (X^{(2)})^T \; \cdots \; (X^{(K)})^T \right] \end{equation}\]

donde:

\[\begin{equation} (X^{(k)})^T \in \mathbb{R}^{J \times I} \end{equation}\]

Por tanto:

\[\begin{equation} X_b \in \mathbb{R}^{J \times KI} \end{equation}\]

La matricización en la vía C se obtiene vectorizando cada bloque y apilándolos por filas:

\[\begin{equation} X_c = \begin{pmatrix} \operatorname{vec}(X^{(1)})^T \\ \operatorname{vec}(X^{(2)})^T \\ \vdots \\ \operatorname{vec}(X^{(K)})^T \end{pmatrix} \end{equation}\]

donde:

\[\begin{equation} \operatorname{vec}(X^{(k)}) \in \mathbb{R}^{IJ} \end{equation}\]

Por tanto:

\[\begin{equation} X_c \in \mathbb{R}^{K \times IJ} \end{equation}\]

La función implementa el siguiente mapeo:

\[\begin{equation} X_a \longrightarrow \{X^{(1)}, \ldots, X^{(K)}\} \longrightarrow (X_b, X_c) \end{equation}\]

permitiendo recuperar las diferentes matricizaciones a partir de una sola representación matricial del arreglo.

reconstruir_desde_Xa <- function(Xa, J, K) {
  bloques <- lapply(seq_len(K), function(k) {
    Xa[, ((k - 1) * J + 1):(k * J), drop = FALSE]
  })

  Xb <- do.call(cbind, lapply(bloques, t))
  Xc <- do.call(rbind, lapply(bloques, c))

  list(Xb = Xb, Xc = Xc)
}

Matricizaciones iniciales

Primero se construyen las tres matricizaciones: \(X_a\),para la primera vía; \(X_b\) para la segunda vía; y \(X_c\) para la tercera vía. De manera que:

\[\begin{equation} X_a = \begin{pmatrix} 1 & 2 & 10 & 20 \\ 3 & 4 & 30 & 40 \\ 5 & 6 & 50 & 60 \end{pmatrix} \qquad X_b = \begin{pmatrix} 1 & 3 & 5 & 10 & 30 & 50 \\ 2 & 4 & 6 & 20 & 40 & 60 \\ \end{pmatrix} \qquad X_c = \begin{pmatrix} 1 & 3 & 5 & 2 & 4 & 6 \\ 10 & 30 & 50 & 20 & 40 & 60 \\ \end{pmatrix} \end{equation}\]

Xa <- matricizar_Xa(lista)
Xb <- matricizar_Xb(lista)
Xc <- matricizar_Xc(lista)

Xa
##    V1 V2 V1 V2
## A1  1  2 10 20
## A2  3  4 30 40
## A3  5  6 50 60
Xb
##    A1 A2 A3 A1 A2 A3
## V1  1  3  5 10 30 50
## V2  2  4  6 20 40 60
Xc
##    [,1] [,2] [,3] [,4] [,5] [,6]
## T1    1    3    5    2    4    6
## T2   10   30   50   20   40   60

Conversión a matrices

Estas líneas convierten los objetos a matrices numéricas, para poder operar con álgebra matricial.

Xa_mat <- as.matrix(Xa)
Xb_mat <- as.matrix(Xb)
Xc_mat <- as.matrix(Xc)

Dimensiones de cada vía

na <- nrow(Xa_mat)
nb <- nrow(Xb_mat)
nc <- nrow(Xc_mat)

na #tamaño de la vía A
## [1] 3
nb #tamaño de la vía B
## [1] 2
nc #tamaño de la vía C
## [1] 2

Construcción de la matriz centradora

Esta matriz corresponde a:

\[\begin{equation} P = I - \frac{1}{n_a}\mathbf 11' \end{equation}\]

como \(na = 3\), se obtiene:

\[\begin{equation} \begin{pmatrix} 2/3 & -1/3 & -1/3 \\ -1/3 & 2/3 & -1/3 \\ -1/3 & -1/3 & 2/3 \\ \end{pmatrix} \end{equation}\]

Su función es centrar las columnas de Xa respecto a la vía A, es decir, restar a cada columna su media.

Port <- diag(na) - (1 / na) * matrix(1, nrow = na, ncol = na)
Port
##            [,1]       [,2]       [,3]
## [1,]  0.6666667 -0.3333333 -0.3333333
## [2,] -0.3333333  0.6666667 -0.3333333
## [3,] -0.3333333 -0.3333333  0.6666667

Centrado en la vía A

la matriz original era:

\[\begin{equation} X_a = \begin{pmatrix} 1 & 2 & 10 & 20 \\ 3 & 4 & 30 & 40 \\ 5 & 6 & 50 & 60 \\ \end{pmatrix} \end{equation}\]

Las medias por columna son: \((3, 4, 30, 40)\)

Entonces la matriz centrada queda:

\[\begin{equation} X_{ac} = \begin{pmatrix} -2 & -2 & -20 & -20 \\ 0 & 0 & 0 & 0 \\ 2 & 2 & 20 & 20 \\ \end{pmatrix} \end{equation}\]

Es decir, a cada valor se le restó la media de su columna.

  Xac <- Port %*% Xa_mat
  Xac
##                 V1            V2            V1            V2
## [1,] -2.000000e+00 -2.000000e+00 -2.000000e+01 -2.000000e+01
## [2,]  2.220446e-16  4.440892e-16  7.105427e-15  7.105427e-15
## [3,]  2.000000e+00  2.000000e+00  2.000000e+01  2.000000e+01

Reconstrucción de las matricizaciones centradas

Aquí se toma Xac y se reconstruyen sus versiones matricizadas en las vías B y C.

La matriz Xac se divide en dos bloques de 2 columnas:

Primer bloque:

\[\begin{equation} \begin{pmatrix} -2 & -2 \\ 0 & 0 \\ 2 & 2 \\ \end{pmatrix} \end{equation}\]

Segundo bloque:

\[\begin{equation} \begin{pmatrix} -20 & -20 \\ 0 & 0 \\ 20 & 20 \\ \end{pmatrix} \end{equation}\]

A partir de ellos se obtienen:

\[\begin{equation} X_{bac} \begin{pmatrix} -2 & 0 & 2 & -20 & 0 & 20\\ -2 & 0 & 2 & -20 & 0 & 20 \\ \end{pmatrix} \end{equation}\]

y

\[\begin{equation} X_{cac} \begin{pmatrix} -2 & 0 & 2 & -2 & 0 & 2\\ -20 & 0 & 20 & -20 & 0 & 20 \\ \end{pmatrix} \end{equation}\]

recons_centrado <- reconstruir_desde_Xa(Xac, J = nb, K = nc)
Xbac <- recons_centrado$Xb
Xcac <- recons_centrado$Xc

Xbac
##    [,1]         [,2] [,3] [,4]         [,5] [,6]
## V1   -2 2.220446e-16    2  -20 7.105427e-15   20
## V2   -2 4.440892e-16    2  -20 7.105427e-15   20
Xcac
##      [,1]         [,2] [,3] [,4]         [,5] [,6]
## [1,]   -2 2.220446e-16    2   -2 4.440892e-16    2
## [2,]  -20 7.105427e-15   20  -20 7.105427e-15   20

Cálculo de la matriz de normalización de la vía B

Paso 1: producto \(X_{bac}X'_{bac}\)

La diagonal de esa matriz contiene las sumas de cuadrados de cada fila de Xbac.

Para la primera fila:

\[\begin{equation} (−2)^2+0^2+2^2+(−20)^2+0^2+20^2=808 \end{equation}\]

La segunda fila da el mismo valor. Entonces:

\[\begin{equation} diag(X_{bac}X'_{bac})= (808, 808) \end{equation}\]

Xbac %*% t(Xbac)
##     V1  V2
## V1 808 808
## V2 808 808
Paso 2: dividir por \(na\cdot nc\)

Como:

\[\begin{equation} na\cdot nc = 3 \cdot 2 = 6 \end{equation}\]

se obtiene:

\[\begin{equation} (808/6,808/6)=(134.6667,134.6667) \end{equation}\]

Paso 3: raíz cuadrada e inverso

\[\begin{equation} \frac{1}{\sqrt{134.6667}}\approx 0.08617 \end{equation}\]

por tanto:

\[\begin{equation} W_{j} \begin{pmatrix} 0.08617 & 0 \\ 0 & 0.08617 \\ \end{pmatrix} \end{equation}\]

Esta matriz diagonal sirve para escalar el modo B.

Wj <- diag(1 / sqrt(diag(Xbac %*% t(Xbac)) / (na * nc)))
Wj
##            [,1]       [,2]
## [1,] 0.08617275 0.00000000
## [2,] 0.00000000 0.08617275

Cálculo de la matriz de normalización de la vía C

Paso 1: producto \(X_{cac}X'_{cac}\)

La diagonal contiene las sumas de cuadrados de cada fila de Xcac.

Primera fila:

\[\begin{equation} (−2)^2+0^2+2^2+(−2)^2+0^2+2^2=16 \end{equation}\]

Segunda fila:

\[\begin{equation} (−20)^2+0^2+20^2+(−20)^2+0^2+20^2=1600 \end{equation}\]

Entonces:

\[\begin{equation} diag(X_{cac}X'_{cac})= (16, 1600) \end{equation}\]

Xcac %*% t(Xcac)
##      [,1] [,2]
## [1,]   16  160
## [2,]  160 1600
Paso 2: dividir por \(na\cdot nb\)

Como:

\[\begin{equation} na\cdot nb = 3 \cdot 2 = 6 \end{equation}\]

resulta:

\[\begin{equation} (16/6,1600/6)=(2.6667,266.6667) \end{equation}\]

Paso 3: raíz cuadrada e inverso

\[\begin{equation} \frac{1}{\sqrt{2.6667}}\approx 0.6124 \\ \frac{1}{\sqrt{266.6667}}\approx 0.0612 \end{equation}\]

Así:

\[\begin{equation} W_{k} \begin{pmatrix} 0.6124 & 0 \\ 0 & 0.0612 \\ \end{pmatrix} \end{equation}\]

Esta matriz diagonal sirve para escalar la vía C.

Wk <- diag(1 / sqrt(diag(Xcac %*% t(Xcac)) / (na * nb)))
Wk
##           [,1]       [,2]
## [1,] 0.6123724 0.00000000
## [2,] 0.0000000 0.06123724

Obtención de la matriz estandarizada final

Aquí se utiliza el producto de Kronecker:

\[\begin{equation} W_k \otimes W_j \end{equation}\]

Como ambas matrices son diagonales, el resultado también es diagonal:

\[\begin{equation} W_k \otimes W_j \approx diag(0.0528,0.0528,0.00528,0.00528) \end{equation}\]

luego:

\[\begin{equation} Y_a = X_{ac}(W_k \otimes W_j) \end{equation}\]

y se obtiene aproximadamente:

\[\begin{equation} Y_{a} = \begin{pmatrix} −0.1056 & −0.1056 & −0.1056 & −0.1056 \\ 0 & 0 & 0 & 0 \\ 0.1056 & 0.1056 & 0.1056 & 0.1056 \\ \end{pmatrix} \end{equation}\]

Interpretación final

La estandarización realizada en este procedimiento consta de dos etapas:

  1. Centrado en la vía A: Se resta, a cada columna de Xa, su promedio entre individuos.

  2. Normalización en las vías B y C Se construyen matrices diagonales Wj y Wk, que ajustan la escala de las variables y de las ocasiones.

Finalmente, la matriz:

Ya <- Xac %*% kronecker(Wk, Wj)
Ya
##               [,1]          [,2]          [,3]          [,4]
## [1,] -1.055396e-01 -1.055396e-01 -1.055396e-01 -1.055396e-01
## [2,]  1.171725e-17  2.343451e-17  3.749521e-17  3.749521e-17
## [3,]  1.055396e-01  1.055396e-01  1.055396e-01  1.055396e-01