Comprensión del método DIFFIT con datos simulados

Introducción

El propósito de este documento es comprender la lógica del método DIFFIT en el contexto del modelo Tucker3. Más que profundizar en todos los aspectos del algoritmo de estimación, el interés principal es mostrar cómo DIFFIT permite seleccionar una dimensionalidad razonable a partir de una colección de modelos estimados.

En un modelo Tucker3, el investigador debe decidir cuántos componentes conservar en cada uno de los tres modos. Es decir, debe elegir una terna:

\[\begin{equation} (P,Q,R), \end{equation}\]

donde:

• \(P\) es el número de componentes del modo A,
• \(Q\) es el número de componentes del modo B,
• \(R\) es el número de componentes del modo C.

A medida que estos valores aumentan, el ajuste del modelo generalmente mejora. Sin embargo, una mejora en ajuste no siempre justifica una solución más compleja. Por ello, se necesita un criterio que permita identificar el punto en el que seguir aumentando componentes deja de aportar una mejora sustancial. Ese es precisamente el papel de DIFFIT.

Objetivo

El objetivo de este documento es entender la lógica y la utilidad del método DIFFIT mediante un ejemplo con datos simulados.

En particular, se busca:

1. simular un arreglo tridimensional con estructura Tucker3;
2. estimar múltiples modelos con distintas ternas \((P,Q,R)\);
3. construir la secuencia de mejores soluciones por cada nivel de complejidad \(S=P+Q+R\);
4. calcular los incrementos de ajuste \(dif(s)\);
5. identificar los máximos secuenciales;
6. seleccionar la solución final mediante DIFFIT;
7. resumir los resultados en una tabla final.

Idea general del método DIFFIT

DIFFIT se apoya en la siguiente lógica:

1. Se estiman muchos modelos Tucker3 para distintas ternas \((P,Q,R)\).
2. Para cada valor de \[\begin{equation} S=P+Q+R, \end{equation}\] se conserva solo el modelo con mayor ajuste.
3. Luego se calcula cuánto mejora el ajuste al pasar de un valor de \(S\) al siguiente. Esa mejora se denota por: \[\begin{equation} dif(s). \end{equation}\]
4. Se identifican aquellos incrementos que siguen siendo grandes en comparación con los incrementos posteriores.
5. Finalmente, se selecciona la solución donde todavía hay un salto importante en el ajuste, pero a partir de la cual los beneficios adicionales se reducen considerablemente.

En términos intuitivos, DIFFIT busca el punto donde la curva de ajuste comienza a mostrar rendimientos decrecientes.

1. Limpieza del entorno

rm(list = ls())

2. Simulación de datos de tres vías

2.1 Dimensiones del arreglo y modelo verdadero

Se considera un arreglo de dimensión \(I\times J \times K = 30 \times 8 \times 6\), de manera que se tienen \(30\) unidade en la vía A, \(8\) unidades en la vía B y \(6\) unidades en la C. Además, se establece como estructura verdadera \((P, Q, R) = (2, 2, 2)\). Es decir, el arreglo \(\mathbf {\underline{X}}\), será generado a partir de dos componentes en cada vía.

El parámetro ruido_sd = \(0.25\) controla la magnitud del ruido aleatorio. Mientras mayor sea este valor, más difícil será recuperar la estructura verdadera.

Ejemplo numérico simple

Si en lugar de \(30 \times 8\times 6\) se trabajara con un arreglo pequeño \(2\times 3\times 2\), entonces cada observación \(X_{IJK}\) estaría indexada por: \(i=1, 2\); \(j=1,2,3\); \(k=1,2\). Y si el modelo verdadero fuera \((1,2,1)\), significaría que la estructura latente se resume con: 1 componente en la primera vía, 2 en la segunda, 1 en la tercera.

# Dimensiones del arreglo
I <- 30
J <- 8
K <- 6

# Modelo verdadero
P_true <- 2
Q_true <- 2
R_true <- 2

# Nivel de ruido
ruido_sd <- 0.25

# Semilla
semilla <- 2026
set.seed(semilla)

2.2 Generación de matrices de componentes verdaderas

Se general las matrices verdaderas de los componentes:

\(\mathbf{A} \in \mathbb{R}^{I\times P}\), \(\mathbf{B} \in \mathbb{R}^{J\times Q}\), y \(\mathbf{C} \in \mathbb{R}^{K\times R}\).

En este caso\(\mathbf{A} \in \mathbb{R}^{30\times 2}\), \(\mathbf{B} \in \mathbb{R}^{8\times 2}\), y \(\mathbf{C} \in \mathbb{R}^{6\times 2}\).

A_true <- qr.Q(qr(matrix(rnorm(I * P_true), I, P_true)))
B_true <- qr.Q(qr(matrix(rnorm(J * Q_true), J, Q_true)))
C_true <- qr.Q(qr(matrix(rnorm(K * R_true), K, R_true)))

La expresión: qr.Q(qr(…)) se usa para ortonormalizar las columnas. Esto implica que, por ejemplo, en \(\mathbf{A}\):

\[\begin{equation} \mathbf{A'A} = \mathbf{I}_p \end{equation}\]

lo mismo ocurre con \(\mathbf{B}\) y \(\mathbf{C}\).

Ejemplo numérico

Si una matriz inicial fuera:

\[\begin{equation} M= \begin{pmatrix} 1 & 2\\ 2 & 1\\ 1 & 1 \end{pmatrix}. \end{equation}\]

la descomposición QR produce una matriz \(Q\) con columnas ortogonales y de norma \(1\). Esa será la base utilizada como matriz de componentes.

2.3 Generación del núcleo verdadero

Aquí se genera el núcleo:

\[\begin{equation} \mathbf{G} \in \mathbb{R}^{p\times Q \times R} \end{equation}\]

Como en este caso \(P=Q=R=2\), entonces:

\[\begin{equation} \mathbf{G} \in \mathbb{R}^{2\times 2 \times 2} \end{equation}\]

El núcleo contiene la relación entre los componentes de las tres vías.

Ejemplo numérico

Si \(P=Q=R=2\), el núcleo tiene 8 elementos:

\[\begin{equation} g_{111}, g_{112}, g_{121}, g_{122}, g_{211}, g_{212}, g_{221}, g_{222} \end{equation}\]

Cada uno de ellos pondera la interacción entre un componente de , uno de y uno de .

G_true <- array(rnorm(P_true * Q_true * R_true, mean = 0, sd = 1),
                dim = c(P_true, Q_true, R_true))

2.4 Reconstrucción de la parte estructural

Este bloque implementa directamente el modelo Tucker3:

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

Cada elemento del arreglo estructural \(\mathbf {\underline{X}}\)_true es una suma ponderada de productos de componentes.

Ejemplo numérico

Supongamos un caso extremadamente simple con: \(P=Q=R=1\)

Entonces el mdoelo se reduce a:

\[\begin{equation} x_{ijk} = a_{i1}b_{j1}c_{k1}g_{111} \end{equation}\]

Si \(a_{i1}= 1, b_{j1}= 3, c_{k1}= 4, g_{111} = 5\)

entonces:

\[\begin{equation} x_{111} = 2\cdot 3 \cdot4 \cdot5 = 120 \end{equation}\]

X_true <- 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_true) {
        for (q in 1:Q_true) {
          for (r in 1:R_true) {
            suma <- suma +
              A_true[i, p] *
              B_true[j, q] *
              C_true[k, r] *
              G_true[p, q, r]
          }
        }
      }
      X_true[i, j, k] <- suma
    }
  }
}

2.5 Agregar ruido

Los datos observados se generan como:

\[\begin{equation} \mathbf {\underline{X}} = \mathbf {\underline{X}}_{true} + \mathbf {\underline{E}} \end{equation}\]

donde es un arreglo de errores aleatorios.

Luego se centra globalmente el arreglo restando su media:

Esto desplaza el arreglo para que su promedio total sea cero.

\[\begin{equation} X \leftarrow X - \bar{X}. \end{equation}\]

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

X <- X_true + E

# Centrado global
X <- X - mean(X)

cat("Modelo verdadero simulado: (",
    P_true, ", ", Q_true, ", ", R_true, ")\n", sep = "")

## Modelo verdadero simulado: (2, 2, 2)

3. Generación de combinaciones válidas

No todas las ternas \((P,Q,R)\) son útiles. Se usan las restricciones:

\[\begin{equation} PQ \geq R, PR \geq Q, QR \geq P \end{equation}\]

Estas evitan modelos redundantes o combinaciones cuya complejidad excede lo que las otras dimensiones pueden representar.

Ejemplo numérico

Si \(P=1\), \(Q=1\), \(R=4\), entonces:

\[\begin{equation} PQ = 1 \cdot 1 = 1 < 4 = R \end{equation}\]

Esa terna no se considera válida.

Pmax <- 4
Qmax <- 4
Rmax <- 4

modelos <- NULL

for (P in 1:Pmax) {
  for (Q in 1:Qmax) {
    for (R in 1:Rmax) {
      if (P * Q >= R && P * R >= Q && Q * R >= P) {
        modelos <- rbind(modelos, c(P, Q, R))
      }
    }
  }
}

modelos <- as.data.frame(modelos)
colnames(modelos) <- c("P", "Q", "R")

modelos

##    P Q R
## 1  1 1 1
## 2  1 2 2
## 3  1 3 3
## 4  1 4 4
## 5  2 1 2
## 6  2 2 1
## 7  2 2 2
## 8  2 2 3
## 9  2 2 4
## 10 2 3 2
## 11 2 3 3
## 12 2 3 4
## 13 2 4 2
## 14 2 4 3
## 15 2 4 4
## 16 3 1 3
## 17 3 2 2
## 18 3 2 3
## 19 3 2 4
## 20 3 3 1
## 21 3 3 2
## 22 3 3 3
## 23 3 3 4
## 24 3 4 2
## 25 3 4 3
## 26 3 4 4
## 27 4 1 4
## 28 4 2 2
## 29 4 2 3
## 30 4 2 4
## 31 4 3 2
## 32 4 3 3
## 33 4 3 4
## 34 4 4 1
## 35 4 4 2
## 36 4 4 3
## 37 4 4 4

4. Estimación de todos los modelos Tucker3

En esta sección se estima cada modelo válido. Aunque el objetivo principal no es estudiar el algoritmo ALS, sí se necesita obtener el ajuste de cada terna para poder aplicar DIFFIT.

Se construyen tres versiones matricizadas de \(X\):

• \(X_A \in \mathbb{R}^{I \times JK}\)
• \(X_B \in \mathbb{R}^{J \times IK}\)
• \(X_C \in \mathbb{R}^{K \times IJ}\)

En cada iteración ALS se actualizan \(A\), \(B\) y \(C\) a partir de descomposiciones SVD. Por ejemplo, para \(A\):

\[\begin{equation} M_A = X_A (C \otimes B), \end{equation}\]

y luego se toman las primeras \(P\) columnas ortonormales de la matriz izquierda de la SVD.

Se calcula mediante:

\[\begin{equation} G = A^\top X_A (C \otimes B). \end{equation}\]

Se usa la suma de cuadrados del error:

\[\begin{equation} SSE = \sum_{i,j,k} (x_{ijk} - \hat{x}_{ijk})^2. \end{equation}\]

El ajuste porcentual se calcula como:

\[\begin{equation} Fit = 100 \left( 1 - \frac{SSE}{\sum x_{ijk}^2} \right). \end{equation}\]

Si la suma de cuadrados total es 200 y el error residual es 50, entonces:

\[\begin{equation} Fit = 100 \left( 1 - \frac{50}{200} \right) = 75\%. \end{equation}\]

resultados <- data.frame()

cat("\nEstimando modelos...\n")

## 
## Estimando modelos...

for (m in 1:nrow(modelos)) {
  
  P <- modelos$P[m]
  Q <- modelos$Q[m]
  R <- modelos$R[m]
  
  cat("Modelo: (", P, ", ", Q, ", ", R, ")\n", sep = "")
  
  XA <- matrix(0, nrow = I, ncol = J * K)
  col_ini <- 1
  for (k in 1:K) {
    XA[, col_ini:(col_ini + J - 1)] <- X[, , k]
    col_ini <- col_ini + J
  }
  
  XB <- matrix(0, nrow = J, ncol = I * K)
  col_ini <- 1
  for (k in 1:K) {
    XB[, col_ini:(col_ini + I - 1)] <- t(X[, , k])
    col_ini <- col_ini + I
  }
  
  XC <- matrix(0, nrow = K, ncol = I * J)
  for (k in 1:K) {
    XC[k, ] <- as.vector(X[, , k])
  }
  
  mejor_sse <- Inf
  mejor_A <- NULL
  mejor_B <- NULL
  mejor_C <- NULL
  mejor_G <- NULL
  mejor_Xhat <- NULL
  mejor_iter <- NA
  
  nstart <- 3
  tol <- 1e-7
  max_iter <- 5
  
  for (inicio in 1:nstart) {
    
    A <- qr.Q(qr(matrix(rnorm(I * P), I, P)))
    B <- qr.Q(qr(matrix(rnorm(J * Q), J, Q)))
    C <- qr.Q(qr(matrix(rnorm(K * R), K, R)))
    
    sse_old <- Inf
    
    for (iter in 1:max_iter) {
      
      M_A <- XA %*% kronecker(C, B)
      sv_A <- svd(M_A)
      A <- sv_A$u[, 1:P, drop = FALSE]
      
      M_B <- XB %*% kronecker(C, A)
      sv_B <- svd(M_B)
      B <- sv_B$u[, 1:Q, drop = FALSE]
      
      M_C <- XC %*% kronecker(B, A)
      sv_C <- svd(M_C)
      C <- sv_C$u[, 1:R, drop = FALSE]
      
      Gmat <- t(A) %*% XA %*% kronecker(C, B)
      
      G <- array(0, dim = c(P, Q, R))
      col_ini <- 1
      for (r in 1:R) {
        G[, , r] <- Gmat[, col_ini:(col_ini + Q - 1)]
        col_ini <- col_ini + Q
      }
      
      Xhat <- 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) {
              for (q in 1:Q) {
                for (r in 1:R) {
                  suma <- suma +
                    A[i, p] *
                    B[j, q] *
                    C[k, r] *
                    G[p, q, r]
                }
              }
            }
            Xhat[i, j, k] <- suma
          }
        }
      }
      
      sse <- sum((X - Xhat)^2)
      
      if (abs(sse_old - sse) < tol) {
        break
      }
      
      sse_old <- sse
    }
    
    if (sse < mejor_sse) {
      mejor_sse <- sse
      mejor_A <- A
      mejor_B <- B
      mejor_C <- C
      mejor_G <- G
      mejor_Xhat <- Xhat
      mejor_iter <- iter
    }
  }
  
  ss_total <- sum(X^2)
  fit <- 100 * (1 - mejor_sse / ss_total)
  
  resultados <- rbind(
    resultados,
    data.frame(
      P = P,
      Q = Q,
      R = R,
      S = P + Q + R,
      Fit = fit,
      SSE = mejor_sse,
      Iter = mejor_iter
    )
  )
}

## Modelo: (1, 1, 1)
## Modelo: (1, 2, 2)
## Modelo: (1, 3, 3)
## Modelo: (1, 4, 4)
## Modelo: (2, 1, 2)
## Modelo: (2, 2, 1)
## Modelo: (2, 2, 2)
## Modelo: (2, 2, 3)
## Modelo: (2, 2, 4)
## Modelo: (2, 3, 2)
## Modelo: (2, 3, 3)
## Modelo: (2, 3, 4)
## Modelo: (2, 4, 2)
## Modelo: (2, 4, 3)
## Modelo: (2, 4, 4)
## Modelo: (3, 1, 3)
## Modelo: (3, 2, 2)
## Modelo: (3, 2, 3)
## Modelo: (3, 2, 4)
## Modelo: (3, 3, 1)
## Modelo: (3, 3, 2)
## Modelo: (3, 3, 3)
## Modelo: (3, 3, 4)
## Modelo: (3, 4, 2)
## Modelo: (3, 4, 3)
## Modelo: (3, 4, 4)
## Modelo: (4, 1, 4)
## Modelo: (4, 2, 2)
## Modelo: (4, 2, 3)
## Modelo: (4, 2, 4)
## Modelo: (4, 3, 2)
## Modelo: (4, 3, 3)
## Modelo: (4, 3, 4)
## Modelo: (4, 4, 1)
## Modelo: (4, 4, 2)
## Modelo: (4, 4, 3)
## Modelo: (4, 4, 4)

resultados <- resultados[order(resultados$S, -resultados$Fit), ]
rownames(resultados) <- NULL

resultados

##    P Q R  S      Fit      SSE Iter
## 1  1 1 1  3 11.39860 96.37621    5
## 2  2 2 1  5 19.03203 88.07294    5
## 3  2 1 2  5 15.99182 91.37993    5
## 4  1 2 2  5 12.68992 94.97157    5
## 5  2 2 2  6 22.41652 84.39147    5
## 6  3 2 2  7 24.44519 82.18477    5
## 7  2 3 2  7 23.05179 83.70044    5
## 8  2 2 3  7 22.76124 84.01649    5
## 9  3 3 1  7 21.87077 84.98510    5
## 10 3 1 3  7 18.49484 88.65728    5
## 11 1 3 3  7 13.48567 94.10600    5
## 12 4 2 2  8 25.77844 80.73453    5
## 13 3 3 2  8 25.48422 81.05457    5
## 14 3 2 3  8 25.10779 81.46403    5
## 15 2 3 3  8 23.71989 82.97372    5
## 16 2 4 2  8 23.13279 83.61234    5
## 17 2 2 4  8 22.79789 83.97663    5
## 18 3 3 3  9 27.81601 78.51817    5
## 19 4 2 3  9 27.39862 78.97218    5
## 20 4 3 2  9 27.30676 79.07211    5
## 21 3 4 2  9 27.03230 79.37065    5
## 22 3 2 4  9 26.12848 80.35377    5
## 23 4 4 1  9 24.48828 82.13790    5
## 24 2 4 3  9 24.25846 82.38789    5
## 25 2 3 4  9 24.05744 82.60655    5
## 26 4 1 4  9 20.52694 86.44685    5
## 27 1 4 4  9 13.80791 93.75548    5
## 28 4 3 3 10 30.91305 75.14936    5
## 29 4 4 2 10 29.71202 76.45578    5
## 30 3 4 3 10 28.91667 77.32092    5
## 31 3 3 4 10 28.62903 77.63381    5
## 32 4 2 4 10 28.22889 78.06906    5
## 33 2 4 4 10 24.79477 81.80452    5
## 34 4 4 3 11 32.67800 73.22953    5
## 35 4 3 4 11 31.59967 74.40249    5
## 36 3 4 4 11 30.07567 76.06022    5
## 37 4 4 4 12 34.33618 71.42585    5

5. Primer paso de DIFFIT: mejor modelo dentro de cada S

DIFFIT no compara directamente todas las ternas. Primero agrupa por:

\[\begin{equation} S = P + Q + R \end{equation}\]

Luego, dentro de cada valor de \(S\), conserva solamente la solución con mayor ajuste.

Ejemplo numérico

Supongamos que para \(S=5\) aparecen tres ternas:\((1, 2, 2)\) con ajuste \(52\), \((2, 1, 2)\) con ajuste \(49\), y \((2, 2, 1)\) con ajuste \(68\).

Entonces, DIFFIT conserva para \(S=5\) únicamente la terna \((2,2,1)\), porque tiene el mayor ajuste.

valores_S <- sort(unique(resultados$S))

mejores_por_S <- data.frame()

for (s in valores_S) {
  sub <- resultados[resultados$S == s, ]
  fila_mejor <- sub[which.max(sub$Fit), ]
  mejores_por_S <- rbind(mejores_por_S, fila_mejor)
}

rownames(mejores_por_S) <- NULL

mejores_por_S

##   P Q R  S      Fit      SSE Iter
## 1 1 1 1  3 11.39860 96.37621    5
## 2 2 2 1  5 19.03203 88.07294    5
## 3 2 2 2  6 22.41652 84.39147    5
## 4 3 2 2  7 24.44519 82.18477    5
## 5 4 2 2  8 25.77844 80.73453    5
## 6 3 3 3  9 27.81601 78.51817    5
## 7 4 3 3 10 30.91305 75.14936    5
## 8 4 4 3 11 32.67800 73.22953    5
## 9 4 4 4 12 34.33618 71.42585    5

6. Segundo paso: cálculo de dif(s)

Para cada valor de \(S\), se mide cuánto mejora el ajuste respecto a la mejor solución del valor anterior:

\[\begin{equation} dif(s) = Fit(s)- Fit(s-1) \end{equation}\]

Para la primera solución se toma:

\[\begin{equation} dif(s_1) = Fit(s_1) \end{equation}\]

pues no hay una solución anterior con la cual compararla.

Si los mejores ajustes son:

library(knitr)
library(kableExtra)

## Warning: package 'kableExtra' was built under R version 4.5.3

tabla <- data.frame(
  S = c(3, 5, 6, 7),
  Fit = c(40, 68, 74, 77)
)

kable(tabla) %>%
  kable_styling(full_width = FALSE, position = "center")

S	Fit
3	40
5	68
6	74
7	77

Entonces:

\[\begin{equation} dif(3) = 40, \qquad dif(5) = 68 - 40 = 28, \qquad dif(6) = 74 - 68 = 6, \qquad dif(7) = 77 - 74 = 3. \end{equation}\]

La lógica aquí es que el paso de \(S = 3\) a \(S = 5\) produjo una mejora muy grande, mientras que después los incrementos son menores.

mejores_por_S$dif <- NA_real_
mejores_por_S$dif[1] <- mejores_por_S$Fit[1]

if (nrow(mejores_por_S) >= 2) {
  for (i in 2:nrow(mejores_por_S)) {
    mejores_por_S$dif[i] <- mejores_por_S$Fit[i] - mejores_por_S$Fit[i - 1]
  }
}

mejores_por_S

##   P Q R  S      Fit      SSE Iter       dif
## 1 1 1 1  3 11.39860 96.37621    5 11.398598
## 2 2 2 1  5 19.03203 88.07294    5  7.633437
## 3 2 2 2  6 22.41652 84.39147    5  3.384482
## 4 3 2 2  7 24.44519 82.18477    5  2.028676
## 5 4 2 2  8 25.77844 80.73453    5  1.333250
## 6 3 3 3  9 27.81601 78.51817    5  2.037567
## 7 4 3 3 10 30.91305 75.14936    5  3.097039
## 8 4 4 3 11 32.67800 73.22953    5  1.764956
## 9 4 4 4 12 34.33618 71.42585    5  1.658174

7. Tercer paso: máximos secuenciales

Un valor \(dif(s)\) es considerado máximo secuencial si es mayor que todos los incrementos que le siguen.

Formalmente, \(dif(s\) es máximo secuencial si:

\[\begin{equation} dif(s) > dif(t) ~~ para ~ todo ~ t > s \end{equation}\]

Esto permite identificar los saltos importantes que siguen dominando a todos los incrementos posteriores.

Ejemplo numérico

Si:

\[\begin{equation} dif = (40, 28, 6, 3, 1) \end{equation}\]

entonces los máximos secuenciales serían:

• 40, porque es mayor que 28, 6, 3 y 1;
• 28, porque es mayor que 6, 3 y 1;
• 6, porque es mayor que 3 y 1;
• 3, porque es mayor que 1;
• 1, porque es el último.

En cambio, si una secuencia fuera:

\[\begin{equation} dif = (40, 5, 8, 2) \end{equation}\]

entonces 5 no sería máximo secuencial, porque más adelante aparece un valor mayor, 8.

es_max_secuencial <- rep(FALSE, nrow(mejores_por_S))

for (i in 1:nrow(mejores_por_S)) {
  if (i == nrow(mejores_por_S)) {
    es_max_secuencial[i] <- TRUE
  } else {
    es_max_secuencial[i] <- all(mejores_por_S$dif[i] > mejores_por_S$dif[(i + 1):nrow(mejores_por_S)])
  }
}

secuenciales <- mejores_por_S[es_max_secuencial, ]
rownames(secuenciales) <- NULL

secuenciales

##   P Q R  S      Fit      SSE Iter       dif
## 1 1 1 1  3 11.39860 96.37621    5 11.398598
## 2 2 2 1  5 19.03203 88.07294    5  7.633437
## 3 2 2 2  6 22.41652 84.39147    5  3.384482
## 4 4 3 3 10 30.91305 75.14936    5  3.097039
## 5 4 4 3 11 32.67800 73.22953    5  1.764956
## 6 4 4 4 12 34.33618 71.42585    5  1.658174

8. Cuarto paso: cálculo de b

DIFFIT compara los máximos secuenciales mediante:

\[\begin{equation} b_t(m) = \frac{dif_t(m)}{dif_t(m+1)}. \end{equation}\]

Si esta razón es grande, significa que el incremento actual es mucho mayor que el siguiente incremento relevante. Eso sugiere una caída fuerte en la ganancia marginal del ajuste.

Si dos máximos secuenciales consecutivos son:

\[\begin{equation} dif_t(1) = 28, \qquad dif_t(2) = 6, \end{equation}\]

entonces:

\[\begin{equation} b_t(1) = \frac{28}{6} = 4.67. \end{equation}\]

Una razón alta como esta indica que el salto asociado al primer valor fue mucho más importante que el siguiente.

secuenciales$bt <- NA_real_

if (nrow(secuenciales) >= 2) {
  for (i in 1:(nrow(secuenciales) - 1)) {
    if (!is.na(secuenciales$dif[i + 1]) && secuenciales$dif[i + 1] != 0) {
      secuenciales$bt[i] <- secuenciales$dif[i] / secuenciales$dif[i + 1]
    }
  }
}

secuenciales

##   P Q R  S      Fit      SSE Iter       dif       bt
## 1 1 1 1  3 11.39860 96.37621    5 11.398598 1.493246
## 2 2 2 1  5 19.03203 88.07294    5  7.633437 2.255422
## 3 2 2 2  6 22.41652 84.39147    5  3.384482 1.092812
## 4 4 3 3 10 30.91305 75.14936    5  3.097039 1.754740
## 5 4 4 3 11 32.67800 73.22953    5  1.764956 1.064397
## 6 4 4 4 12 34.33618 71.42585    5  1.658174       NA

9. Quinto paso: umbral T

El umbral se define como:

\[\begin{equation} T = \frac{100}{s_{\max} - 3}, \end{equation}\]

donde:

\[\begin{equation} s_{\max} = \min(I, JK) + \min(J, IK) + \min(K, IJ). \end{equation}\]

Este umbral evita seleccionar un incremento demasiado pequeño, incluso si su razón \(b_t\) fuera alta.

No basta con que un salto sea grande comparado con el siguiente: también debe ser sustancial en términos absolutos.

smax <- min(I, J * K) + min(J, I * K) + min(K, I * J)
T_porcentaje <- 100 / (smax - 3)

T_porcentaje

## [1] 2.439024

10. Sexto paso: selección final DIFFIT

DIFFIT selecciona la solución cuya razón \(b_t\) es máxima entre aquellas que además cumplen:

\[\begin{equation} dif_t(m) > T. \end{equation}\]

Es decir, la solución elegida debe:

• corresponder a un incremento relevante;
• y mostrar una caída importante en la ganancia de ajuste a partir de ese punto.

candidatos <- secuenciales[!is.na(secuenciales$bt) & secuenciales$dif > T_porcentaje, ]

if (nrow(candidatos) == 0) {
  candidatos <- secuenciales[!is.na(secuenciales$bt), ]
}

if (nrow(candidatos) == 0) {
  stop("No se pudo seleccionar una solución con DIFFIT.")
}

seleccionado <- candidatos[which.max(candidatos$bt), ]

seleccionado

##   P Q R S      Fit      SSE Iter      dif       bt
## 2 2 2 1 5 19.03203 88.07294    5 7.633437 2.255422

11. Tabla resumen general

tabla_resumen <- data.frame(
  Concepto = c(
    "Modelo verdadero",
    "Modelo seleccionado por DIFFIT",
    "S verdadero",
    "S seleccionado",
    "Fit del modelo seleccionado (%)",
    "SSE del modelo seleccionado",
    "dif(s) seleccionado",
    "bt seleccionado",
    "Umbral T (%)"
  ),
  Valor = c(
    paste0("(", P_true, ", ", Q_true, ", ", R_true, ")"),
    paste0("(", seleccionado$P, ", ", seleccionado$Q, ", ", seleccionado$R, ")"),
    P_true + Q_true + R_true,
    seleccionado$S,
    round(seleccionado$Fit, 4),
    round(seleccionado$SSE, 6),
    round(seleccionado$dif, 6),
    round(seleccionado$bt, 6),
    round(T_porcentaje, 6)
  ),
  stringsAsFactors = FALSE
)

tabla_resumen

##                          Concepto     Valor
## 1                Modelo verdadero (2, 2, 2)
## 2  Modelo seleccionado por DIFFIT (2, 2, 1)
## 3                     S verdadero         6
## 4                  S seleccionado         5
## 5 Fit del modelo seleccionado (%)    19.032
## 6     SSE del modelo seleccionado 88.072938
## 7             dif(s) seleccionado  7.633437
## 8                 bt seleccionado  2.255422
## 9                    Umbral T (%)  2.439024

12. Gráficos

op <- par(no.readonly = TRUE)
par(mfrow = c(1, 2))

plot(
  mejores_por_S$S, mejores_por_S$Fit,
  type = "b", pch = 19,
  xlab = "S = P + Q + R",
  ylab = "Fit (%)",
  main = "Mejor ajuste por S"
)

points(seleccionado$S, seleccionado$Fit, pch = 19, cex = 1.5)
text(
  seleccionado$S, seleccionado$Fit,
  labels = paste0(" (", seleccionado$P, ",", seleccionado$Q, ",", seleccionado$R, ")"),
  pos = 4
)

plot(
  mejores_por_S$S, mejores_por_S$dif,
  type = "b", pch = 19,
  xlab = "S = P + Q + R",
  ylab = "dif(s)",
  main = "DIFFIT con datos simulados"
)

abline(h = T_porcentaje, lty = 2)

points(seleccionado$S, seleccionado$dif, pch = 19, cex = 1.5)
text(
  seleccionado$S, seleccionado$dif,
  labels = paste0(" (", seleccionado$P, ",", seleccionado$Q, ",", seleccionado$R, ")"),
  pos = 4
)

par(op)

13. Tabla resumen de los máximos secuenciales de DIFFIT

Esta es la tabla que resume los elementos centrales del criterio.

tabla_maximos_secuenciales <- secuenciales[, c("P", "Q", "R", "S", "Fit", "dif", "bt")]
tabla_maximos_secuenciales

##   P Q R  S      Fit       dif       bt
## 1 1 1 1  3 11.39860 11.398598 1.493246
## 2 2 2 1  5 19.03203  7.633437 2.255422
## 3 2 2 2  6 22.41652  3.384482 1.092812
## 4 4 3 3 10 30.91305  3.097039 1.754740
## 5 4 4 3 11 32.67800  1.764956 1.064397
## 6 4 4 4 12 34.33618  1.658174       NA