Ando falto de ideas, no sé si es la pandemia, el teletrabajo ( o la esclavitud en tiempos modernos como me gusta llamarlo) u otra cosa. Total, que me he puesto a bichear un post antiguo de mi amigo Carlos Gil sobre NMF (factorización no negativa de matrices). Cómo siempre el lo cuenta mucho mejor que yo.
Total, que puede que en breve me toque tener algo a lo que quizá se pueda aplicar este tipo de técnicas, a saber, tener clientes y productos.
De hecho voy a usar su mismo ejemplo.
Nota: La librería es mejor si se instala desde BioConductor con BiocManager.
La librería NMF está bastante bien, utiliza paralelización, por debajo está escrita en C, pero tiene el incoveniente de que aún no está implementado un método predict para nuevos datos
library(MASS)
library(NMF) # BiocManager::install("NMF")
a <- as.matrix(caith)
res <- nmf(a, rank = 2)
a
#> fair red medium dark black
#> blue 326 38 241 110 3
#> light 688 116 584 188 4
#> medium 343 84 909 412 26
#> dark 98 48 403 681 85El caso es factorizar esta matriz en dos matrices no negativas. Otras formas de factorizar esta matriz podría ser con análisis de correspondencias simples. Pero vamos a la descomposición por nmf
w <- res@fit@W
h <- res@fit@Hw
#> [,1] [,2]
#> blue 34.41271 386.48565
#> light 32.29480 896.90510
#> medium 438.78664 576.80290
#> dark 720.24703 5.38822h
#> fair red medium dark black
#> [1,] 0.1141435 0.06226416 0.6667132 0.8736829 9.626828e-02
#> [2,] 0.7049219 0.11239402 0.7074373 0.1715769 1.971208e-13Reconstruimos a
w %*% h
#> fair red medium dark black
#> blue 276.37021 45.58136 296.3578 96.37782 3.312852
#> light 635.93433 102.81758 656.0355 182.10365 3.108965
#> medium 456.68567 92.14988 700.5967 482.32648 42.241237
#> dark 86.00979 45.45118 484.0100 630.19204 69.336946Y si comparamos con a
a - (w %*% h)
#> fair red medium dark black
#> blue 49.62979 -7.581356 -55.35776 13.622177 -0.3128523
#> light 52.06567 13.182418 -72.03547 5.896346 0.8910349
#> medium -113.68567 -8.149881 208.40327 -70.326485 -16.2412371
#> dark 11.99021 2.548819 -81.01004 50.807961 15.6630545Bueno, la reconstrucción no es perfecta, pero bueno, no está tan mal.
Bien, tal y como cuenta Carlos en su entrada ahora podemos normalizar las filas de W y de H, de forma que tengamos probabilidades. Dónde entonces H sería funciones de probabilidad sobre las filas de la matriz original y W serán ponderaciones. O como dice él, H es un menú de preferencias (imaginemos que tenemos usuarios en filas y productos en columnas), en este caso hemos hecho una reducción de dimensión para quedarnos en 2 preferencias, (sería el equivalente conceptual al número de componentes en un PCA o en un CA), y W serían las ponderaciones que cada usuario da a cada una de las preferencias (sus coordenadas en un correspondencias siguiendo el símil)
Normalicemos
w_hat <- w / rowSums(w)
w_hat
#> [,1] [,2]
#> blue 0.08176014 0.918239864
#> light 0.03475549 0.965244506
#> medium 0.43205116 0.567948839
#> dark 0.99257448 0.007425522h_hat <- h / rowSums(h)
h_hat
#> fair red medium dark black
#> [1,] 0.06295585 0.03434180 0.3677257 0.4818799 5.309678e-02
#> [2,] 0.41555704 0.06625716 0.4170398 0.1011460 1.162043e-13Así, el primer “menú” está compuesto por los “productos” fair, red, etc, en proporción a como indica la primera fila de h_hat. Y el individuo “blue” prefiere el primer menú en casi un 0.9 de probabilidad vs un alrededor de 0.1 de preferencia del menú 2. En un PCA diríamos que esos son los “loadings”.
Las filas de W a veces se asocian con arquetipos o individuos promedio. Los individuos “blue” tienen esos pesos los dos factores latentes.
En este caso dónde tenemos color de ojos (fila) y color del pelo (columnas), vemos que en las dos distribuciones multinomiales que hay en sendas filas de h_hat (si, eso es lo que son, dos distribuciones multinomiales), la probabilidad de tener el pelo negro es bastante pequeña (tiene que ver con la tabla de contingencia original, hay muy pocos con el pelo negro). Pero vemos que hay un arquetipo, (el de os ojos oscuros) para el cual el peso que da al menú de preferencias dónde la probabilidad de tener el pelo negro sea mayor. Es decir, al final es una mixtura de distribuciones multinomiales.
En realidad lo que hace NMF es descubrir la estructura subyacente de los datos en un espacio de menor dimensión que el original. Bueno, pues con W y H normalizadas podemos construir una matriz diagonal D que simplemente nos genere muestras de individuos y en qué columnas caen.
Podemos utilizar como matriz diagonal la suma de las filas de a, y así obtener
(d <- diag(rowSums(a)))
#> [,1] [,2] [,3] [,4]
#> [1,] 718 0 0 0
#> [2,] 0 1580 0 0
#> [3,] 0 0 1774 0
#> [4,] 0 0 0 1315Y podemos hacer \(A \approx D \cdot W \cdot H\)
d %*% w_hat %*% h_hat
#> fair red medium dark black
#> [1,] 277.67093 45.69909 296.5397 94.97332 3.116981
#> [2,] 637.21749 102.93373 656.2149 180.71812 2.915739
#> [3,] 466.94392 93.07840 702.0314 471.24977 40.696489
#> [4,] 86.22994 45.47111 484.0408 629.95432 69.303794Que se parece bastante a w %*% h , o podríamos usar otra D, en este caso para obtener qué matriz se obtendría para 10 casos de cada fila.
d <- diag(c(10,10, 10, 10))
d %*% w_hat %*% h_hat
#> fair red medium dark black
#> [1,] 3.867283 0.6364776 4.130079 1.322748 0.04341200
#> [2,] 4.033022 0.6514793 4.153259 1.143786 0.01845405
#> [3,] 2.632153 0.5246809 3.957336 2.656425 0.22940524
#> [4,] 0.655741 0.3457879 3.680919 4.790527 0.52702505y bueno, la verdad es que me pregunto si esto se parece o no a un análisis de correspondencias. Veamos
library(FactoMineR)
res_ca <- CA (a, ncp = 2, graph = FALSE)
factoextra::fviz_ca(res_ca)
Lo primero que hay que darse cuenta es que ambas técnicas no son del todo comparables, el correspondencias busca encontrar dimensiones que expliquen la mayor cantidad de inercia (distancia Chi-cuadrado) y es parecido al PCA en el sentido de que la primera dimensión es la que más explica, etc.. De hecho el CA, diagonaliza 2 matrices derivadas de la tabla de contingencia, una la de los perfiles filas y otra la de los perfiles columna. Y las pinta junta de acuerdo a algún teorema baricéntrico que tuve que demostrar en algún examen allá por los lejanos 90’s.
Pero en realidad si nos fijamos en las coordenadas de las filas en el CA
res_ca$row$coord
#> Dim 1 Dim 2
#> blue -0.40029985 0.16541100
#> light -0.44070764 0.08846303
#> medium 0.03361434 -0.24500190
#> dark 0.70273880 0.13391383No es más que ver las filas en un subespacio (el calculado por el CA) del espacio definido por las columnas y de forma análoga pasa con las columnas. Estas coordenadas podrían ser una forma de codificar la variable categórica. Cabe preguntarse si tienen relación con la estructura obtenida por el NMF.
coordenadas_filas <- cbind(res_ca$row$coord, w_hat)
colnames(coordenadas_filas)[3:4] <- paste0("nmf_", 1:2)
coordenadas_filas
#> Dim 1 Dim 2 nmf_1 nmf_2
#> blue -0.40029985 0.16541100 0.08176014 0.918239864
#> light -0.44070764 0.08846303 0.03475549 0.965244506
#> medium 0.03361434 -0.24500190 0.43205116 0.567948839
#> dark 0.70273880 0.13391383 0.99257448 0.007425522y
cor(coordenadas_filas)
#> Dim 1 Dim 2 nmf_1 nmf_2
#> Dim 1 1.00000000 -0.05155554 0.99991352 -0.99991352
#> Dim 2 -0.05155554 1.00000000 -0.04510151 0.04510151
#> nmf_1 0.99991352 -0.04510151 1.00000000 -1.00000000
#> nmf_2 -0.99991352 0.04510151 -1.00000000 1.00000000Resultado coherente, ¿no? . En este ejemplo de juguete una única dimensión del correspondencias explica el 86,5% de la inercia.
Cosas buenas del nmf.
Ah, se me olvidaba. ¿qué pasa si tengo una nueva fila/usario?, la librería NMF no permite predecir, y aunque se podría implementar, buscando un poco se encuentra la forma
library(NNLM)
res_nmf2 <- nnmf(a, k = 2, loss = "mkl")res_nmf2$W
#> [,1] [,2]
#> blue 0.6928050 13.3441587
#> light 0.1273311 30.9726568
#> medium 13.0923170 19.8721597
#> dark 22.4117740 0.1060281
res_nmf2$H
#> fair red medium dark black
#> [1,] 3.740643 2.012339 21.49645 28.091524 3.09335848
#> [2,] 20.516743 3.311351 21.09274 5.764011 0.08766068Nuevas filas
b <- matrix( data = c(20, 30, 40, 0,20, 10, 10, 30,10, 90), nrow=2, byrow = TRUE)
colnames(b) <- colnames(a)
rownames(b) <- c("tipo_n1", "tipo_n2")
b
#> fair red medium dark black
#> tipo_n1 20 30 40 0 20
#> tipo_n2 10 10 30 10 90Y tiene un método predict
predict(res_nmf2, newdata = b, which = "W")
#> $coefficients
#> [,1] [,2]
#> tipo_n1 0.734150 1.32159
#> tipo_n2 2.566985 0.00000
#>
#> $n.iteration
#> [1] 30
#>
#> $error
#> MSE MKL target.error
#> 1305.05201 26.74997 26.74997
#>
#> $options
#> $options$method
#> [1] "scd"
#>
#> $options$loss
#> [1] "mkl"
#>
#> $options$max.iter
#> [1] 10000
#>
#> $options$rel.tol
#> [1] 1e-12
#>
#>
#> $call
#> nnlm(x = t(object$H), y = t(newdata), method = method, loss = loss)
#>
#> attr(,"class")
#> [1] "nnlm"Lo dicho, una técnica muy interesante y útil.