Ejemplillo con NMF

estadística
correspondencias
factorización
nmf
2020
Author

José Luis Cañadas Reche

Published

October 21, 2020

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

Mostrar / ocultar código
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    85

El 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

Mostrar / ocultar código
w <-  res@fit@W
h <-  res@fit@H
Mostrar / ocultar código
w 
#>             [,1]      [,2]
#> blue    34.41271 386.48565
#> light   32.29480 896.90510
#> medium 438.78664 576.80290
#> dark   720.24703   5.38822
Mostrar / ocultar código
h
#>           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-13

Reconstruimos a

Mostrar / ocultar código
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.336946

Y si comparamos con a

Mostrar / ocultar código
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.6630545

Bueno, 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

Mostrar / ocultar código
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.007425522
Mostrar / ocultar código
h_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-13

Así, 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

Mostrar / ocultar código
(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 1315

Y podemos hacer \(A \approx D \cdot W \cdot H\)

Mostrar / ocultar código
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.303794

Que 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.

Mostrar / ocultar código
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.52702505

y bueno, la verdad es que me pregunto si esto se parece o no a un análisis de correspondencias. Veamos

Mostrar / ocultar código
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

Mostrar / ocultar código
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.13391383

No 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.

Mostrar / ocultar código
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.007425522

y

Mostrar / ocultar código
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.00000000

Resultado 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

Mostrar / ocultar código
library(NNLM)

res_nmf2 <-  nnmf(a, k = 2, loss = "mkl")
Mostrar / ocultar código
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.08766068

Nuevas filas

Mostrar / ocultar código

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    90

Y tiene un método predict

Mostrar / ocultar código
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.