Simulación y estimación del modelo Tucker3 en R

1 Introducción

El análisis de datos de tres vías constituye una extensión natural de los métodos multivariados clásicos cuando la información no se organiza únicamente en una matriz de dos dimensiones, sino en un arreglo tridimensional. Este tipo de estructura aparece con frecuencia en aplicaciones empíricas donde se observa simultáneamente un conjunto de individuos, variables y ocasiones o condiciones.

Dentro de este marco, el modelo Tucker3 representa uno de los enfoques más importantes para estudiar la estructura subyacente de dichos arreglos. Su propósito es resumir la información del arreglo original mediante un conjunto reducido de componentes para cada modo y un núcleo central que describe la interacción entre ellos.

En este documento se presenta, de manera reproducible, un procedimiento completo para:

1. simular datos de tres vías bajo un modelo Tucker3,
2. construir sus matricizaciones,
3. estimar modelos Tucker3 candidatos en R,
4. seleccionar la dimensionalidad mediante el scree plot multiway,
5. interpretar los resultados de ajuste.

2 Planteamiento del modelo Tucker3

Sea un arreglo de datos tridimensional

\[ \mathcal{X} = (x_{ijk}), \qquad i = 1,\dots,I,\; j = 1,\dots,J,\; k = 1,\dots,K, \]

donde:

• \(I\) representa el número de individuos,
• \(J\) representa el número de variables,
• \(K\) representa el número de ocasiones o condiciones.

El modelo Tucker3 expresa cada elemento del arreglo como:

\[ x_{ijk} = \sum_{p=1}^{P} \sum_{q=1}^{Q} \sum_{r=1}^{R} g_{pqr}\, a_{ip}\, b_{jq}\, c_{kr} + e_{ijk}, \]

donde:

• \(A = (a_{ip})\) es la matriz de componentes del modo de individuos,
• \(B = (b_{jq})\) es la matriz de componentes del modo de variables,
• \(C = (c_{kr})\) es la matriz de componentes del modo de ocasiones,
• \(\mathcal{G} = (g_{pqr})\) es el arreglo núcleo,
• \(e_{ijk}\) es el término residual.

En términos conceptuales, el modelo reduce el arreglo original de dimensión \(I \times J \times K\) a una representación más parsimoniosa definida por tres matrices de componentes y un núcleo central.

2.1 Ejemplo numérico ilustrativo

Supóngase un arreglo muy pequeño con:

• \(I=2\) individuos,
• \(J=2\) variables,
• \(K=2\) ocasiones,

y que el modelo verdadero posee dimensión:

• \(P=1\),
• \(Q=1\),
• \(R=1\).

En ese caso, la expresión del modelo se reduce a:

\[ x_{ijk} = g_{111}\, a_{i1}\, b_{j1}\, c_{k1} + e_{ijk}. \]

Si se tuviera, por ejemplo,

\[ g_{111}=2,\qquad a_{11}=1,\; a_{21}=3,\qquad b_{11}=2,\; b_{21}=4,\qquad c_{11}=1,\; c_{21}=5, \]

entonces:

• para \(i=1, j=1, k=1\),

\[ x_{111} = 2(1)(2)(1)=4, \]

• para \(i=2, j=2, k=2\),

\[ x_{222} = 2(3)(4)(5)=120. \]

Este ejemplo muestra que, incluso con una estructura muy reducida, el modelo Tucker3 permite generar un arreglo tridimensional completo a partir de factores latentes.

3 Simulación de datos bajo un modelo Tucker3

En esta sección se simulan datos a partir de un modelo Tucker3 verdadero. La idea es generar un arreglo \(\mathcal{X}\) de dimensión \(I \times J \times K\), con rangos verdaderos \((P,Q,R)\), y posteriormente verificar si el procedimiento de estimación es capaz de recuperar una estructura cercana a la generadora.

3.1 Definición de dimensiones y parámetros

rm(list = ls())

set.seed(123)

tol <- 1e-10
max_iter <- 100

# Dimensiones del arreglo
I <- 24   # individuos
J <- 6    # variables
K <- 4    # ocasiones

# Rangos verdaderos del modelo Tucker3
P_verd <- 2
Q_verd <- 2
R_verd <- 1

En este ejemplo se considera un arreglo con:

• 24 individuos,
• 6 variables,
• 4 ocasiones.

Además, se asume que la estructura verdadera subyacente corresponde a un modelo Tucker3 con dimensionalidad \((2,2,1)\).

3.2 Etiquetas de filas y columnas

ind <- c(
  "AMA","ANC","APU","ARE","AYA","CAJ","CUS","HUV","HUC","ICA","JUN","LAL",
  "LAM","LIM","LOR","MMD","MOQ","PAS","PIU","PUN","SAN","TAC","TUM","UCA"
)

var <- paste0("Var", 1:J)
cond <- paste0("Ocasion", 1:K)

Estas etiquetas son útiles para interpretar los resultados y para mantener una estructura parecida a la empleada en aplicaciones reales.

3.3 Generación de matrices factoriales verdaderas

A0 <- qr.Q(qr(matrix(rnorm(I * P_verd), I, P_verd)))
B0 <- qr.Q(qr(matrix(rnorm(J * Q_verd), J, Q_verd)))
C0 <- qr.Q(qr(matrix(rnorm(K * R_verd), K, R_verd)))

A0

##              [,1]         [,2]
##  [1,] -0.12187608 -0.125248315
##  [2,] -0.05005236 -0.371046013
##  [3,]  0.33894293  0.147062710
##  [4,]  0.01533213  0.032461511
##  [5,]  0.02811377 -0.257709601
##  [6,]  0.37294287  0.236014978
##  [7,]  0.10022676  0.083433020
##  [8,] -0.27508903 -0.033334888
##  [9,] -0.14935695  0.217770463
## [10,] -0.09690971  0.207754569
## [11,]  0.26617801  0.152067586
## [12,]  0.07824193  0.144636523
## [13,]  0.08714822  0.113469763
## [14,]  0.02406808 -0.016690915
## [15,] -0.12086830 -0.054053376
## [16,]  0.38856633 -0.131104324
## [17,]  0.10825816 -0.168103748
## [18,] -0.42764318  0.004230988
## [19,]  0.15251065 -0.300900791
## [20,] -0.10280904  0.496956046
## [21,] -0.23219950  0.297511913
## [22,] -0.04739890 -0.245398008
## [23,] -0.22310586 -0.063594524
## [24,] -0.15849824 -0.085507220

B0

##             [,1]        [,2]
## [1,] -0.48793669  0.04309609
## [2,]  0.05215467 -0.66649214
## [3,] -0.15847298  0.66838171
## [4,]  0.01785850 -0.25720176
## [5,]  0.02681924 -0.05179705
## [6,] -0.85618094 -0.19586010

C0

##            [,1]
## [1,] -0.3054290
## [2,]  0.4041311
## [3,]  0.2680732
## [4,]  0.8194681

Aquí se construyen las matrices de componentes verdaderas:

• \(A_0\) de dimensión \(I \times P\),
• \(B_0\) de dimensión \(J \times Q\),
• \(C_0\) de dimensión \(K \times R\).

Se utiliza la descomposición QR para obtener columnas ortonormales, lo cual facilita la interpretación y evita problemas numéricos.

3.4 Generación del núcleo verdadero

G_array <- array(
  rnorm(P_verd * Q_verd * R_verd, sd = 2),
  dim = c(P_verd, Q_verd, R_verd)
)

G_array

## , , 1
## 
##            [,1]      [,2]
## [1,] -2.1435825 0.8964196
## [2,]  0.6070573 0.1060085

El arreglo núcleo \(\mathcal{G}\) contiene las interacciones entre los componentes de los tres modos.

En este caso, como \((P,Q,R)=(2,2,1)\), el núcleo tiene dimensión \(2 \times 2 \times 1\).

4 Construcción del arreglo tridimensional

El arreglo observado \(\mathcal{X}\) se genera según el modelo Tucker3:

\[ x_{ijk} = \sum_{p=1}^{P} \sum_{q=1}^{Q} \sum_{r=1}^{R} g_{pqr}\, a_{ip}\, b_{jq}\, c_{kr}. \]

4.1 Cálculo del arreglo sin ruido

X_array <- array(0, dim = c(I, J, K))

for (i in 1:I) {
  for (j in 1:J) {
    for (k in 1:K) {
      suma <- 0
      for (p in 1:P_verd) {
        for (q in 1:Q_verd) {
          for (r in 1:R_verd) {
            suma <- suma + G_array[p, q, r] * A0[i, p] * B0[j, q] * C0[k, r]
          }
        }
      }
      X_array[i, j, k] <- suma
    }
  }
}

Este bloque implementa literalmente la definición matemática del modelo. Para cada celda \(x_{ijk}\), se suma el producto de:

• el elemento del núcleo \(g_{pqr}\),
• el peso del individuo \(a_{ip}\),
• el peso de la variable \(b_{jq}\),
• el peso de la ocasión \(c_{kr}\).

4.2 Adición de ruido aleatorio

sigma_error <- 0.20

E_array <- array(
  rnorm(I * J * K, mean = 0, sd = sigma_error),
  dim = c(I, J, K)
)

X_array <- X_array + E_array

En aplicaciones empíricas, los datos observados rara vez siguen exactamente la estructura teórica. Por ello, se incorpora un término de error aleatorio con desviación estándar igual a 0.20.

5 Organización del arreglo en una lista de matrices

En muchos casos, los datos de tres vías se almacenan como una colección de matrices \(I \times J\), una por cada ocasión. Esta representación resulta útil para la lectura, almacenamiento e inspección de la información.

lista_datos <- vector("list", K)

for (k in 1:K) {
  mat_k <- X_array[, , k]
  rownames(mat_k) <- ind
  colnames(mat_k) <- var
  lista_datos[[k]] <- as.data.frame(mat_k)
}

names(lista_datos) <- cond

lista_datos[[1]]

##            Var1         Var2         Var3        Var4         Var5
## AMA  0.21366944  0.019853062  0.055113783  0.12972106 -0.086827167
## ANC  0.39354646 -0.140842899 -0.116661233 -0.05841075  0.059279629
## APU -0.19738409  0.347306083 -0.265994899 -0.28585978  0.137587510
## ARE -0.46402117 -0.116343963 -0.208970961 -0.30151173 -0.096376454
## AYA  0.16887949  0.440487585  0.013472441 -0.31929155  0.105114048
## CAJ -0.24435738  0.390122401 -0.294609835 -0.07437422  0.084852510
## CUS -0.16337084 -0.024434612 -0.126205782 -0.28370275 -0.040169809
## HUV  0.29327029 -0.265272855  0.027405837  0.11482799  0.004436993
## HUC  0.01191877 -0.171842486  0.413286805  0.40885022 -0.012271791
## ICA -0.19353587  0.032857764 -0.100992475 -0.26432135  0.421929750
## JUN -0.03836743  0.010133987 -0.028072473  0.18016705 -0.140319735
## LAL -0.04081309 -0.050836657 -0.005724087  0.16095863 -0.217192469
## LAM -0.01760841 -0.170093694 -0.216482773  0.07416566  0.009949865
## LIM  0.06759661 -0.003990533 -0.021292413 -0.19978277  0.062915089
## LOR -0.03890825 -0.183807965  0.323150926 -0.03408643  0.083646408
## MMD -0.01151823 -0.250975882 -0.022136412 -0.02483076 -0.078907707
## MOQ -0.09493238 -0.054595722 -0.024098135  0.12064392 -0.208675062
## PAS  0.20839327  0.091210716  0.038160967 -0.10958120  0.239049472
## PIU  0.14204496 -0.085615089 -0.456711777  0.20640818 -0.064097677
## PUN  0.16535861  0.105240198  0.259595600 -0.08086506 -0.178003526
## SAN  0.03823219 -0.370333950 -0.223240974  0.19296837 -0.055606512
## TAC  0.22360409 -0.024302852  0.159681356 -0.21495846 -0.040130925
## TUM  0.26694161  0.054793338  0.445304434 -0.27066976  0.215112139
## UCA  0.15456662  0.024877343 -0.243990929  0.63476442  0.012198434
##              Var6
## AMA  0.1919159783
## ANC -0.1357409521
## APU -0.1046519143
## ARE -0.0673504311
## AYA -0.0378798267
## CAJ -0.3291641372
## CUS -0.2991941300
## HUV  0.5333895445
## HUC  0.2318062530
## ICA -0.1668302173
## JUN -0.2320614025
## LAL -0.2528811380
## LAM  0.4143168379
## LIM  0.2475263529
## LOR -0.0006811923
## MMD -0.1099798008
## MOQ -0.1654989856
## PAS  0.1222330113
## PIU -0.2847076188
## PUN  0.0152349807
## SAN  0.4970055141
## TAC -0.0272902403
## TUM  0.1264487171
## UCA  0.1149682726

Cada elemento de lista_datos corresponde a una ocasión y contiene la información de los individuos por variables.

5.1 Ejemplo conceptual de esta organización

Si \(K=3\), entonces el arreglo tridimensional puede verse como tres “rebanadas”:

\[ X^{(1)},\quad X^{(2)},\quad X^{(3)}, \]

donde cada \(X^{(k)}\) es una matriz \(I \times J\).
En otras palabras, el tercer modo actúa como un índice que identifica cada una de las matrices.

6 Matricizaciones del arreglo

La estimación de Tucker3 suele apoyarse en las matricizaciones del arreglo tridimensional. Estas consisten en reorganizar los mismos datos en matrices de dos dimensiones, preservando la información original.

6.1 Matricización por el modo A

La primera matricización, denotada aquí como \(X_a\), apila horizontalmente las matrices \(I \times J\) correspondientes a cada ocasión:

\[ X_a \in \mathbb{R}^{I \times JK}. \]

Xa <- do.call(cbind, lapply(1:K, function(k) {
  as.matrix(lista_datos[[k]])
}))
rownames(Xa) <- ind
colnames(Xa) <- as.vector(sapply(cond, function(cc) paste(var, cc, sep = "_")))

dim(Xa)

## [1] 24 24

Xa[1:5, 1:8]

##     Var1_Ocasion1 Var2_Ocasion1 Var3_Ocasion1 Var4_Ocasion1 Var5_Ocasion1
## AMA     0.2136694    0.01985306    0.05511378    0.12972106   -0.08682717
## ANC     0.3935465   -0.14084290   -0.11666123   -0.05841075    0.05927963
## APU    -0.1973841    0.34730608   -0.26599490   -0.28585978    0.13758751
## ARE    -0.4640212   -0.11634396   -0.20897096   -0.30151173   -0.09637645
## AYA     0.1688795    0.44048759    0.01347244   -0.31929155    0.10511405
##     Var6_Ocasion1 Var1_Ocasion2 Var2_Ocasion2
## AMA    0.19191598    0.20783783   0.000322177
## ANC   -0.13574095   -0.08141971   0.103990117
## APU   -0.10465191   -0.06727334  -0.034608404
## ARE   -0.06735043    0.33803366  -0.161213511
## AYA   -0.03787983   -0.04553665  -0.161721640

6.2 Matricización por el modo B

La segunda matricización, \(X_b\), reorganiza la información de manera que las filas correspondan a las variables:

\[ X_b \in \mathbb{R}^{J \times IK}. \]

Xb <- do.call(cbind, lapply(1:K, function(k) {
  t(as.matrix(lista_datos[[k]]))
}))
rownames(Xb) <- var
colnames(Xb) <- as.vector(sapply(cond, function(cc) paste(ind, cc, sep = "_")))

dim(Xb)

## [1]  6 96

Xb[1:6, 1:8]

##      AMA_Ocasion1 ANC_Ocasion1 APU_Ocasion1 ARE_Ocasion1 AYA_Ocasion1
## Var1   0.21366944   0.39354646   -0.1973841  -0.46402117   0.16887949
## Var2   0.01985306  -0.14084290    0.3473061  -0.11634396   0.44048759
## Var3   0.05511378  -0.11666123   -0.2659949  -0.20897096   0.01347244
## Var4   0.12972106  -0.05841075   -0.2858598  -0.30151173  -0.31929155
## Var5  -0.08682717   0.05927963    0.1375875  -0.09637645   0.10511405
## Var6   0.19191598  -0.13574095   -0.1046519  -0.06735043  -0.03787983
##      CAJ_Ocasion1 CUS_Ocasion1 HUV_Ocasion1
## Var1  -0.24435738  -0.16337084  0.293270289
## Var2   0.39012240  -0.02443461 -0.265272855
## Var3  -0.29460983  -0.12620578  0.027405837
## Var4  -0.07437422  -0.28370275  0.114827994
## Var5   0.08485251  -0.04016981  0.004436993
## Var6  -0.32916414  -0.29919413  0.533389545

6.3 Matricización por el modo C

La tercera matricización, \(X_c\), organiza la información de forma que las filas correspondan a las ocasiones:

\[ X_c \in \mathbb{R}^{K \times IJ}. \]

Xc <- matrix(NA, nrow = K, ncol = I * J)

for (k in 1:K) {
  Xc[k, ] <- as.vector(as.matrix(lista_datos[[k]]))
}

rownames(Xc) <- cond
colnames(Xc) <- as.vector(sapply(ind, function(ii) paste(ii, var, sep = "_")))

dim(Xc)

## [1]   4 144

Xc[1:4, 1:10]

##             AMA_Var1     AMA_Var2    AMA_Var3   AMA_Var4    AMA_Var5
## Ocasion1  0.21366944  0.393546463 -0.19738409 -0.4640212  0.16887949
## Ocasion2  0.20783783 -0.081419708 -0.06727334  0.3380337 -0.04553665
## Ocasion3  0.05664332  0.007848808 -0.40613180  0.5162115 -0.01273822
## Ocasion4 -0.19876516 -0.154549321  0.47145191  0.1560811 -0.21525745
##              AMA_Var6    ANC_Var1   ANC_Var2    ANC_Var3    ANC_Var4
## Ocasion1 -0.244357380 -0.16337084  0.2932703  0.01191877 -0.19353587
## Ocasion2 -0.008966587 -0.21315798 -0.3735876 -0.20592513  0.05663515
## Ocasion3  0.220217436  0.07737071  0.1275606  0.10308212 -0.08637670
## Ocasion4  0.256024948 -0.11005083 -0.6506734 -0.15476372 -0.15162309

6.4 Ejemplo numérico sencillo de matricización

Supóngase que se tiene un arreglo \(2 \times 2 \times 2\) dado por:

\[ X^{(1)} = \begin{pmatrix} 1 & 2\\ 3 & 4 \end{pmatrix}, \qquad X^{(2)} = \begin{pmatrix} 5 & 6\\ 7 & 8 \end{pmatrix}. \]

Entonces:

• la matricización por el modo A sería

\[ X_a = \begin{pmatrix} 1 & 2 & 5 & 6\\ 3 & 4 & 7 & 8 \end{pmatrix}, \]

• la matricización por el modo B sería

\[ X_b = \begin{pmatrix} 1 & 3 & 5 & 7\\ 2 & 4 & 6 & 8 \end{pmatrix}, \]

• la matricización por el modo C sería

\[ X_c = \begin{pmatrix} 1 & 2 & 3 & 4\\ 5 & 6 & 7 & 8 \end{pmatrix}. \]

Este ejemplo muestra que las matricizaciones no cambian la información, sino únicamente su forma de organización.

7 Estructura de entrada para la estimación

El código de estimación trabaja con una lista llamada X_pre, que contiene las tres matricizaciones:

X_pre <- list(Xa, Xb, Xc)

A partir de esta lista se extraen nombres y matrices para la etapa de estimación.

8 Estimación del modelo Tucker3

8.1 Preparación de objetos de trabajo

lista_est <- X_pre

varcond <- colnames(lista_est[[1]])
indcond <- colnames(lista_est[[2]])
indvar  <- colnames(lista_est[[3]])
cond    <- rownames(lista_est[[3]])
var     <- rownames(lista_est[[2]])
ind     <- rownames(lista_est[[1]])

X1 <- as.matrix(lista_est[[1]])
X2 <- as.matrix(lista_est[[2]])
X3 <- as.matrix(lista_est[[3]])

Aquí:

• X1 = Xa,
• X2 = Xb,
• X3 = Xc.

Estas son las matrices que se utilizarán para la descomposición.

8.2 Valores singulares iniciales

sv_X1 <- svd(X1)
sv_X2 <- svd(X2)
sv_X3 <- svd(X3)

autovalores_iniciales <- list(
  X1_valores_singulares = sv_X1$d,
  X2_valores_singulares = sv_X2$d,
  X3_valores_singulares = sv_X3$d,
  X1_autovalores = sv_X1$d^2,
  X2_autovalores = sv_X2$d^2,
  X3_autovalores = sv_X3$d^2
)

autovalores_iniciales

## $X1_valores_singulares
##  [1] 2.97533643 1.83269029 1.68586201 1.45207219 1.37798403 1.29375781
##  [7] 1.26629037 1.17877410 0.99212608 0.90827490 0.77166523 0.71826535
## [13] 0.69717353 0.65055897 0.57424164 0.54325747 0.45864338 0.41837126
## [19] 0.34430969 0.31612417 0.21724291 0.15846054 0.14386989 0.04359331
## 
## $X2_valores_singulares
## [1] 3.248758 2.281935 2.124688 1.835627 1.673293 1.560726
## 
## $X3_valores_singulares
## [1] 3.536287 2.495977 2.371524 2.126478
## 
## $X1_autovalores
##  [1] 8.852626875 3.358753689 2.842130715 2.108513637 1.898839988 1.673809275
##  [7] 1.603491297 1.389508371 0.984314149 0.824963294 0.595467224 0.515905111
## [13] 0.486050928 0.423226971 0.329753464 0.295128679 0.210353749 0.175034510
## [19] 0.118549163 0.099934490 0.047194482 0.025109744 0.020698544 0.001900377
## 
## $X2_autovalores
## [1] 10.554431  5.207229  4.514298  3.369525  2.799910  2.435866
## 
## $X3_autovalores
## [1] 12.505322  6.229903  5.624126  4.521908

La descomposición en valores singulares de cada matricización proporciona una primera visión de la estructura de los datos.

Debe notarse que:

• sv_X1$d, sv_X2$d, sv_X3$d son valores singulares,
• mientras que sus cuadrados corresponden a los autovalores asociados de \(X'X\) o \(XX'\), según el caso.

9 Construcción del conjunto de modelos candidatos

El procedimiento de estimación evalúa múltiples combinaciones posibles de \((P,Q,R)\) sujetas a restricciones de factibilidad algebraica.

P <- min(6, nrow(X1) - 1)
Q <- min(5, nrow(X2) - 1)
R <- min(4, nrow(X3) - 1)

c1 <- 1:R
c2 <- 1:Q
c3 <- 1:P

a1 <- rep(R, Q)
a2 <- rep(Q * R, P)

s1 <- rep(c2, a1)
w1 <- rep(c1, P * Q)
w2 <- rep(s1, P)
w3 <- rep(c3, a2)

Z <- cbind(w3, w2, w1)

pq <- w3 * w2
qr <- w2 * w1
pr <- w3 * w1

L <- pq >= w1 & qr >= w3 & pr >= w2
Z <- Z[L, , drop = FALSE]

Ncomp <- rowSums(Z)
Z <- cbind(Z, Ncomp)

maxi <- min((P + 1), ((Q + 1) * (R + 1))) +
  min((Q + 1), ((P + 1) * (R + 1))) +
  min((R + 1), ((P + 1) * (Q + 1))) - 3

if (maxi == (P + Q + R)) {
  Z <- Z[Z[, 4] != max(Z[, 4]), , drop = FALSE]
} else {
  Z <- Z[Z[, 4] < maxi, , drop = FALSE]
}

nZ <- nrow(Z)

if (nZ == 0) stop("No quedaron modelos válidos después del filtrado.")

orden_modelos <- order(Z[, 4])

Z

##       w3 w2 w1 Ncomp
##  [1,]  1  1  1     3
##  [2,]  1  2  2     5
##  [3,]  1  3  3     7
##  [4,]  2  1  2     5
##  [5,]  2  2  1     5
##  [6,]  2  2  2     6
##  [7,]  2  2  3     7
##  [8,]  2  3  2     7
##  [9,]  2  3  3     8
## [10,]  2  4  2     8
## [11,]  2  4  3     9
## [12,]  2  5  3    10
## [13,]  3  1  3     7
## [14,]  3  2  2     7
## [15,]  3  2  3     8
## [16,]  3  3  1     7
## [17,]  3  3  2     8
## [18,]  3  3  3     9
## [19,]  3  4  2     9
## [20,]  3  4  3    10
## [21,]  3  5  2    10
## [22,]  3  5  3    11
## [23,]  4  2  2     8
## [24,]  4  2  3     9
## [25,]  4  3  2     9
## [26,]  4  3  3    10
## [27,]  4  4  1     9
## [28,]  4  4  2    10
## [29,]  4  4  3    11
## [30,]  4  5  2    11
## [31,]  4  5  3    12
## [32,]  5  2  3    10
## [33,]  5  3  2    10
## [34,]  5  3  3    11
## [35,]  5  4  2    11
## [36,]  5  4  3    12
## [37,]  5  5  1    11
## [38,]  5  5  2    12
## [39,]  5  5  3    13
## [40,]  6  2  3    11
## [41,]  6  3  2    11
## [42,]  6  3  3    12
## [43,]  6  4  2    12
## [44,]  6  4  3    13
## [45,]  6  5  2    13

La matriz Z contiene, en cada fila, una combinación candidata de dimensiones:

\[ (P,Q,R,\; P+Q+R). \]

10 Algoritmo iterativo de estimación

La estimación se realiza mediante un procedimiento alternante basado en descomposiciones SVD. Para cada modelo candidato, se actualizan sucesivamente las matrices de componentes hasta que la diferencia en el error cuadrático residual sea suficientemente pequeña.

SCE <- rep(NA_real_, nZ)
A_list  <- vector("list", nZ)
B_list  <- vector("list", nZ)
C_list  <- vector("list", nZ)
G1_list <- vector("list", nZ)
G2_list <- vector("list", nZ)
G3_list <- vector("list", nZ)

D_A_list <- vector("list", nZ)
D_B_list <- vector("list", nZ)
D_C_list <- vector("list", nZ)

tiempo_total <- system.time({
  
  for (i in seq_len(nZ)) {
    p_i <- Z[i, 1]
    q_i <- Z[i, 2]
    r_i <- Z[i, 3]
    
    A_svd <- svd(X1)
    B_svd <- svd(X2)
    C_svd <- svd(X3)
    
    sce_prev <- Inf
    
    for (iter in seq_len(max_iter)) {
      
      A_svd <- svd(
        X1 %*% kronecker(
          C_svd$u[, 1:r_i, drop = FALSE],
          B_svd$u[, 1:q_i, drop = FALSE]
        )
      )
      
      B_svd <- svd(
        X2 %*% kronecker(
          C_svd$u[, 1:r_i, drop = FALSE],
          A_svd$u[, 1:p_i, drop = FALSE]
        )
      )
      
      C_svd <- svd(
        X3 %*% kronecker(
          B_svd$u[, 1:q_i, drop = FALSE],
          A_svd$u[, 1:p_i, drop = FALSE]
        )
      )
      
      G1_curr <- t(A_svd$u[, 1:p_i, drop = FALSE]) %*% X1 %*%
        kronecker(
          C_svd$u[, 1:r_i, drop = FALSE],
          B_svd$u[, 1:q_i, drop = FALSE]
        )
      
      G2_curr <- t(B_svd$u[, 1:q_i, drop = FALSE]) %*% X2 %*%
        kronecker(
          C_svd$u[, 1:r_i, drop = FALSE],
          A_svd$u[, 1:p_i, drop = FALSE]
        )
      
      G3_curr <- t(C_svd$u[, 1:r_i, drop = FALSE]) %*% X3 %*%
        kronecker(
          B_svd$u[, 1:q_i, drop = FALSE],
          A_svd$u[, 1:p_i, drop = FALSE]
        )
      
      X1_est <- A_svd$u[, 1:p_i, drop = FALSE] %*% G1_curr %*%
        kronecker(
          t(C_svd$u[, 1:r_i, drop = FALSE]),
          t(B_svd$u[, 1:q_i, drop = FALSE])
        )
      
      sce_curr <- sum((X1 - X1_est)^2)
      
      if (abs(sce_prev - sce_curr) <= tol) {
        A_list[[i]]  <- A_svd$u[, 1:p_i, drop = FALSE]
        B_list[[i]]  <- B_svd$u[, 1:q_i, drop = FALSE]
        C_list[[i]]  <- C_svd$u[, 1:r_i, drop = FALSE]
        G1_list[[i]] <- G1_curr
        G2_list[[i]] <- G2_curr
        G3_list[[i]] <- G3_curr
        
        D_A_list[[i]] <- A_svd$d[1:p_i]
        D_B_list[[i]] <- B_svd$d[1:q_i]
        D_C_list[[i]] <- C_svd$d[1:r_i]
        
        SCE[i] <- sce_curr
        break
      }
      
      sce_prev <- sce_curr
      
      if (iter == max_iter) {
        A_list[[i]]  <- A_svd$u[, 1:p_i, drop = FALSE]
        B_list[[i]]  <- B_svd$u[, 1:q_i, drop = FALSE]
        C_list[[i]]  <- C_svd$u[, 1:r_i, drop = FALSE]
        G1_list[[i]] <- G1_curr
        G2_list[[i]] <- G2_curr
        G3_list[[i]] <- G3_curr
        
        D_A_list[[i]] <- A_svd$d[1:p_i]
        D_B_list[[i]] <- B_svd$d[1:q_i]
        D_C_list[[i]] <- C_svd$d[1:r_i]
        
        SCE[i] <- sce_curr
      }
    }
  }
})

tiempo_total <- tiempo_total[["elapsed"]]
tiempo_total

## [1] 0.56

10.1 Interpretación del algoritmo

El esquema anterior puede entenderse como una optimización alternante:

1. se fija inicialmente una aproximación para \(B\) y \(C\),
2. se actualiza \(A\),
3. luego se actualiza \(B\),
4. después se actualiza \(C\),
5. con estas matrices se calcula el núcleo,
6. se reconstruye \(X_1\),
7. se mide el error cuadrático residual.

Este proceso se repite hasta que el cambio en el error sea menor que la tolerancia prefijada.

11 Selección del modelo mediante scree plot multiway

11.1 Construcción de la tabla de resultados

SPM <- cbind(Z, SCE = SCE)
SPM <- SPM[orden_modelos, , drop = FALSE]
SPM <- as.data.frame(SPM)

colnames(SPM) <- c("P", "Q", "R", "No de comp.", "SCE")

SSX <- sum(X1^2)

SCEmin <- tapply(SPM[, "SCE"], SPM[, "No de comp."], min)
ncomp <- as.numeric(names(SCEmin))
ajuste_scree <- (1 - SCEmin / SSX) * 100

ord <- order(ncomp)
ncomp <- ncomp[ord]
SCEmin <- SCEmin[ord]
ajuste_scree <- ajuste_scree[ord]

SPM

##    P Q R No de comp.      SCE
## 1  1 1 1           3 20.72044
## 2  1 2 2           5 20.18617
## 3  2 1 2           5 19.55081
## 4  2 2 1           5 19.20127
## 5  2 2 2           6 18.72687
## 6  1 3 3           7 20.09103
## 7  2 2 3           7 17.96469
## 8  2 3 2           7 17.90240
## 9  3 1 3           7 18.77056
## 10 3 2 2           7 17.64519
## 11 3 3 1           7 18.08197
## 12 2 3 3           8 17.29984
## 13 2 4 2           8 18.00058
## 14 3 2 3           8 16.83524
## 15 3 3 2           8 16.42407
## 16 4 2 2           8 16.64289
## 17 2 4 3           9 17.11437
## 18 3 3 3           9 15.43771
## 19 3 4 2           9 15.92384
## 20 4 2 3           9 15.91861
## 21 4 3 2           9 15.69789
## 22 4 4 1           9 17.26455
## 23 2 5 3          10 17.01004
## 24 3 4 3          10 15.27538
## 25 3 5 2          10 15.88223
## 26 4 3 3          10 13.92960
## 27 4 4 2          10 14.98498
## 28 5 2 3          10 15.12342
## 29 5 3 2          10 15.29794
## 30 3 5 3          11 14.85087
## 31 4 4 3          11 13.34463
## 32 4 5 2          11 14.74612
## 33 5 3 3          11 13.30096
## 34 5 4 2          11 14.13562
## 35 5 5 1          11 16.68590
## 36 6 2 3          11 14.73029
## 37 6 3 2          11 14.52234
## 38 4 5 3          12 13.00840
## 39 5 4 3          12 12.42269
## 40 5 5 2          12 13.83765
## 41 6 3 3          12 12.68152
## 42 6 4 2          12 13.52279
## 43 5 5 3          13 11.84499
## 44 6 4 3          13 11.35981
## 45 6 5 2          13 13.10074

La tabla SPM resume todos los modelos evaluados. Para cada suma total de componentes

\[ S = P + Q + R, \]

se identifica el menor error cuadrático residual.

11.2 Scree plot basado en porcentaje de ajuste

plot(
  ncomp, ajuste_scree, type = "b", pch = 19,
  xlab = "Suma del número de componentes S = P + Q + R",
  ylab = "Porcentaje de ajuste",
  main = "Scree plot multiway"
)

11.3 Scree plot basado en SCE

plot(
  ncomp, SCEmin, type = "b", pch = 19,
  xlab = "Suma del número de componentes S = P + Q + R",
  ylab = "SCE",
  main = "Scree plot multiway basado en SCE"
)

11.4 Interpretación del scree plot

La lógica del scree plot multiway consiste en identificar un punto a partir del cual el incremento en la complejidad del modelo ya no produce mejoras sustanciales en el ajuste. Ese punto suele denominarse codo.

En términos sustantivos, el criterio busca equilibrar:

• parsimonia, y
• capacidad explicativa.

12 Tabla resumen del mejor modelo para cada \(S\)

tabla_mejor_por_S <- do.call(
  rbind,
  lapply(split(SPM, SPM$`No de comp.`), function(df) {
    df_mejor <- df[which.min(df$SCE), , drop = FALSE]
    data.frame(
      S = df_mejor$`No de comp.`,
      P = df_mejor$P,
      Q = df_mejor$Q,
      R = df_mejor$R,
      Varianza_explicada = (1 - df_mejor$SCE / SSX) * 100,
      SCE = df_mejor$SCE
    )
  })
)

tabla_mejor_por_S <- tabla_mejor_por_S[order(tabla_mejor_por_S$S), ]
tabla_mejor_por_S$Varianza_explicada <- round(tabla_mejor_por_S$Varianza_explicada, 4)
tabla_mejor_por_S$SCE <- signif(tabla_mejor_por_S$SCE, 8)

tabla_mejor_por_S

##     S P Q R Varianza_explicada      SCE
## 3   3 1 1 1            28.2564 20.72044
## 5   5 2 2 1            33.5165 19.20127
## 6   6 2 2 2            35.1591 18.72687
## 7   7 3 2 2            38.9044 17.64519
## 8   8 3 3 2            43.1324 16.42407
## 9   9 3 3 3            46.5476 15.43771
## 10 10 4 3 3            51.7694 13.92960
## 11 11 5 3 3            53.9461 13.30096
## 12 12 5 4 3            56.9870 12.42269
## 13 13 6 4 3            60.6672 11.35981

Esta tabla resulta especialmente útil para fines de reporte académico, ya que muestra la mejor combinación \((P,Q,R)\) para cada nivel de complejidad total \(S\).

13 Selección automática del modelo

if (length(ajuste_scree) >= 3) {
  d1 <- diff(ajuste_scree)
  d2 <- diff(d1)
  codo_idx <- which.max(abs(d2)) + 1
  S_opt <- ncomp[codo_idx]
} else {
  S_opt <- ncomp[which.max(ajuste_scree)]
}

fila_seleccionada <- rownames(
  SPM[
    SPM[, "No de comp."] == S_opt &
      SPM[, "SCE"] == min(SPM[SPM[, "No de comp."] == S_opt, "SCE"]),
    , drop = FALSE
  ]
)[1]

j_sel <- as.integer(fila_seleccionada)
jo <- orden_modelos[j_sel]

def <- as.integer(SPM[j_sel, 1:4])

p <- def[1]
q <- def[2]
r <- def[3]

p

## [1] 2

q

## [1] 2

r

## [1] 1

Aquí se selecciona automáticamente el modelo asociado al codo de la curva.

14 Extracción del modelo elegido

A <- A_list[[jo]]
B <- B_list[[jo]]
C <- C_list[[jo]]
G1 <- G1_list[[jo]]
G2 <- G2_list[[jo]]
G3 <- G3_list[[jo]]

autoval_A <- D_A_list[[jo]]
autoval_B <- D_B_list[[jo]]
autoval_C <- D_C_list[[jo]]

res_modelo <- list(
  A = A, B = B, C = C,
  Ga = G1, Gb = G2, Gc = G3,
  Xa = X1,
  varcond = varcond,
  indcond = indcond,
  indvar = indvar,
  ind = ind,
  var = var,
  cond = cond,
  autoval_A = autoval_A,
  autoval_B = autoval_B,
  autoval_C = autoval_C,
  autoval_A_cuadrado = autoval_A^2,
  autoval_B_cuadrado = autoval_B^2,
  autoval_C_cuadrado = autoval_C^2
)

str(res_modelo, max.level = 1)

## List of 19
##  $ A                 : num [1:24, 1:2] -0.09037 -0.00468 0.30858 0.07339 0.00274 ...
##  $ B                 : num [1:6, 1:2] -0.4467 0.1166 -0.4705 0.0805 -0.0855 ...
##  $ C                 : num [1:4, 1] -0.424 0.487 0.212 0.734
##  $ Ga                : num [1:2, 1:2] -2.86 5.85e-08 1.84e-07 1.24
##  $ Gb                : num [1:2, 1:2] -2.86 1.84e-07 5.85e-08 1.24
##  $ Gc                : num [1, 1:4] -2.86 5.85e-08 1.84e-07 1.24
##  $ Xa                : num [1:24, 1:24] 0.214 0.394 -0.197 -0.464 0.169 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ varcond           : chr [1:24] "Var1_Ocasion1" "Var2_Ocasion1" "Var3_Ocasion1" "Var4_Ocasion1" ...
##  $ indcond           : chr [1:96] "AMA_Ocasion1" "ANC_Ocasion1" "APU_Ocasion1" "ARE_Ocasion1" ...
##  $ indvar            : chr [1:144] "AMA_Var1" "AMA_Var2" "AMA_Var3" "AMA_Var4" ...
##  $ ind               : chr [1:24] "AMA" "ANC" "APU" "ARE" ...
##  $ var               : chr [1:6] "Var1" "Var2" "Var3" "Var4" ...
##  $ cond              : chr [1:4] "Ocasion1" "Ocasion2" "Ocasion3" "Ocasion4"
##  $ autoval_A         : num [1:2] 2.86 1.24
##  $ autoval_B         : num [1:2] 2.86 1.24
##  $ autoval_C         : num 3.11
##  $ autoval_A_cuadrado: num [1:2] 8.15 1.53
##  $ autoval_B_cuadrado: num [1:2] 8.15 1.53
##  $ autoval_C_cuadrado: num 9.68

15 Medidas de ajuste del modelo seleccionado

AJtotal <- (sum(res_modelo$Ga^2) / sum(res_modelo$Xa^2)) * 100

contrib.compsA <- (
  diag(as.matrix(res_modelo$Ga) %*% t(as.matrix(res_modelo$Ga))) /
    sum(res_modelo$Xa^2)
) * 100

contrib.compsB <- (
  diag(as.matrix(res_modelo$Gb) %*% t(as.matrix(res_modelo$Gb))) /
    sum(res_modelo$Xa^2)
) * 100

contrib.compsC <- (
  diag(as.matrix(res_modelo$Gc) %*% t(as.matrix(res_modelo$Gc))) /
    sum(res_modelo$Xa^2)
) * 100

contrib.triadacomps <- (as.matrix(res_modelo$Ga)^2) * (100 / sum(res_modelo$Xa^2))

res_ajuste <- list(
  AJtotal = AJtotal,
  contrib.compsA = contrib.compsA,
  contrib.compsB = contrib.compsB,
  contrib.compsC = contrib.compsC,
  contrib.triadacomps = contrib.triadacomps
)

res_ajuste

## $AJtotal
## [1] 33.51651
## 
## $contrib.compsA
## [1] 28.230532  5.285976
## 
## $contrib.compsB
## [1] 28.230532  5.285976
## 
## $contrib.compsC
## [1] 33.51651
## 
## $contrib.triadacomps
##              [,1]         [,2]
## [1,] 2.823053e+01 1.176603e-13
## [2,] 1.184114e-14 5.285976e+00

15.1 Interpretación del ajuste total

La medida

\[ AJ_{total} = \frac{\sum g_{pqr}^2}{\sum x_{ijk}^2}\times 100 \]

representa el porcentaje de la variabilidad total del arreglo original que es explicado por el modelo seleccionado.

Valores altos indican que la estructura Tucker3 retenida captura una gran proporción de la información del arreglo.

16 Resumen final del modelo seleccionado

S_modelo <- p + q + r
SCE_modelo <- SPM[j_sel, "SCE"]
SCEmin_modelo <- SCEmin[as.character(S_modelo)]
ajuste_modelo <- AJtotal

resumen_modelo <- data.frame(
  Indicador = c(
    "Criterio de selección",
    "Fila seleccionada en SPM",
    "Modelo elegido (P, Q, R)",
    "Número total de componentes S",
    "P (modo A)",
    "Q (modo B)",
    "R (modo C)",
    "SCE del modelo elegido",
    "SCE mínimo para ese S",
    "Porcentaje de ajuste (%)",
    "Modelo verdadero simulado"
  ),
  Valor = c(
    "Scree plot multiway",
    j_sel,
    paste0("(", p, ", ", q, ", ", r, ")"),
    S_modelo,
    p,
    q,
    r,
    signif(SCE_modelo, 8),
    signif(SCEmin_modelo, 8),
    round(ajuste_modelo, 4),
    paste0("(", P_verd, ", ", Q_verd, ", ", R_verd, ")")
  ),
  stringsAsFactors = FALSE
)

resumen_modelo

##                        Indicador               Valor
## 1          Criterio de selección Scree plot multiway
## 2       Fila seleccionada en SPM                   4
## 3       Modelo elegido (P, Q, R)           (2, 2, 1)
## 4  Número total de componentes S                   5
## 5                     P (modo A)                   2
## 6                     Q (modo B)                   2
## 7                     R (modo C)                   1
## 8         SCE del modelo elegido           19.201269
## 9          SCE mínimo para ese S           19.201269
## 10      Porcentaje de ajuste (%)             33.5165
## 11     Modelo verdadero simulado           (2, 2, 1)

Esta tabla permite comparar el modelo seleccionado con el modelo verdadero utilizado en la simulación.

17 Discusión de resultados

En un experimento de simulación como el desarrollado aquí, el interés principal radica en evaluar si el procedimiento de estimación y selección es capaz de recuperar una estructura cercana a la verdadera.

Dado que los datos fueron generados a partir de un Tucker3 con dimensión \((2,2,1)\), se espera que el modelo finalmente seleccionado sea igual o próximo a dicha estructura, aunque pequeñas discrepancias pueden surgir debido al ruido incorporado en la simulación.

Desde una perspectiva metodológica, este tipo de ejercicio es útil porque permite:

• validar el algoritmo de estimación,
• estudiar el comportamiento del criterio de selección,
• analizar la estabilidad del ajuste frente a diferentes niveles de ruido y diferentes dimensiones verdaderas.

18 Conclusiones

En este documento se presentó un procedimiento completo y reproducible para la simulación y estimación del modelo Tucker3 en R. El desarrollo incluyó la construcción del arreglo tridimensional, sus matricizaciones, la estimación iterativa de múltiples modelos candidatos y la selección de dimensionalidad mediante el scree plot multiway.

Desde el punto de vista didáctico, los ejemplos numéricos mostraron cómo un arreglo de tres vías puede representarse a partir de componentes latentes y un núcleo central, mientras que desde el punto de vista aplicado, el script desarrollado constituye una base útil para estudios de simulación más amplios.

En términos generales, este enfoque resulta apropiado para investigar la recuperación de estructuras latentes en datos tridimensionales y para evaluar el desempeño de criterios de selección de dimensionalidad en contextos controlados.

19 Anexo: observaciones técnicas

1. En el código se almacenan como autoval_A, autoval_B y autoval_C los elementos d provenientes de la SVD. Estrictamente, estos son valores singulares.
2. Los autovalores asociados corresponden a sus cuadrados:

\[ \lambda = d^2. \]

3. Si se desea profundizar el análisis, puede verse:

   • Reorganización de datos para el análisis de tres vías.
   • Estandarización en datos de tres vías.