Modelos Mixtos Gaussianos, Algoritmo EM y una aplicación a la detección de cáncer de mama
Daniela Pinto Veizaga y Julio Lezama Amastalli
17/12/2020
Resumen
En este trabajo introducimos los Modelos Gaussianos Mixtos (GMM, por sus siglas en inglés), como una técnica de clasificación que, además de asignar etiquetas de pertenencia a una categoría también asigna probabilidades de pertenencia. Dado que la función de verosimilitud de un modelo GMM es complicada, este tipo de modelos se estiman mediante el Algoritmo EM (Expectation - Maximization). Este algoritmo consta de dos pasos llamados paso E y paso M. En el paso E, tomando como fijas las estimaciones de los parámetros, se calcula el valor esperado condicional de la función de verosimilitud respecto a la distribución de las variables latentes. En el paso M, tomando como dado el valor esperado de las variables latentes, se calculan los parámetros de interés. Además, dado que los GMM tienden a ser multimodales, es importante considerar la inicialización de los mismos, para no quedar atrapados en un máximo local, por ello introducimos tres técnicas de inicialización: 1) basado en distancia Euclidiana, 2) k - means y 3) Aleatorio. Finalmente, aplicamos estas técnicas al conjunto de datos de pacientes con diagnóstico de cáncer maligno o benigno, con el fin de construir un clasificador para este tema. Los resultados son un error en el conjunto de prueba de 4%, con un error de validación cruzada de 5%.
Introducción
Las técnicas para crear clusters tienen como objetivo identificar grupos cuyos elementos tienen características en común dentro de los grupos y características distintas entre los grupos. Un algoritmo muy famoso para identificar clusters es el algoritmo \(k\)-means, que mediante la definición de una distancia logra identificar distintos clusters que son homogéneos al interior de ellos y heterogéneos entre ellos.
En la siguiente Figura se presenta un ejemplo del algoritmo k - means, donde hay definidos dos clusters identificados por los colores azul y rojo. Así como una visualización del algoritmo GMM.
Ejemplo K - means
Destacamos las siguientes dos características de este algoritmo:
Se asigna una y solo una etiqueta a cada observación.
Las fronteras de los clusters están definidas por la distancia que se emplea.
Sin embargo, en casos menos separados, como en el mostrado en el lado izquierdo de la siguiente Figura, muchas veces resulta útil asignar probabilidades de pertenencia a cada cluster, en vez de directamente asignar una etiqueta. Lo anterior lo podemos hacer con Modelos Mixtos Gaussianos (GMM por sus siglas en inglés).
Los GMM son combinaciones lineales convexas de densidades normales, donde cada densidad es asociada con un grupo o cluster. -en la siguiente Figura se ilustran las diferencias entre los algoritmos k - means y Modelos Gaussianos Mixtos.
Diferencias entre K - means GMM
Modelos Gaussianos Mixtos
En esta Sección presentamos formalmente los modelos Gaussianos Mixtos.
Sea X = \(\{X_1, X_2, \ldots, X_n\}\) una muestra i.i.d. proveniente de una distribución de probabilidad mixta de \(K\) componentes, definida por:
\[ f(x_i|\pi, \theta)=\sum_{k=1}^{K}\pi_{k}\phi_k(x|\theta_k). \] Donde
\(K\) es el número de clusters.
\(0<\pi_k\leq 1\) \(\forall k \in\{1,\ldots,K\}\) y \(\sum_{k=1}^{K}\pi_{k}=1\).
\(\theta_k\) es el conjunto de parámetros de la función de la \(k\)- ésima función de densidad.
A \(\pi_k\) lo interpretamos como la probabilidad de que la observación \(X_i\) provenga de la densidad \(\phi_k\).
Es importante notar que, la función de verosimilitud de un GMM es sumamente complicada en términos de optimización.
Ejemplo
Considere el caso con \(K=5\) clusters y distribiciones bivariadas. La especificación es:
\[ f(x,y|\theta)=\pi_1N(\mu_1, \Sigma_1^2)+ \pi_2N(\mu_2, \Sigma_2^2)+\pi_3N(\mu_3, \Sigma_3^2)+ \pi_4N(\mu_4, \Sigma_4^2)+\pi_5N(\mu_5, \Sigma_5^2) \]
Por lo que debemos estimar
- \(\pi_1\), \(\pi_2\), \(\pi_3\), \(\pi_4\) y \(\pi_5\).
- Medias y varianzas \((\mu_i, \Sigma_i^2)\) para \(i=1,\ldots,5\).
Algoritmo EM
Intuición
Paso E, fijamos los parámetros (en este caso, las medias y varianzas) y resolvemos para la distribución posterior de las variables latentes (en este caso \(z_{ik}\), las etiquetas de clase).
Paso M, fijamos la distribución posterior de las variables latentes (en este caso \(z_{ik}\)), y optimizamos para las medias y varianzas, usando el valor esperado de las variables latentes (en este caso las probabilidades de \(z_{ik}\).
Finalmente, repetir hasta satisfacer algún criterio de convergencia.
Algoritmo EM
Formalidad
El algoritmo EM maximiza una cota inferior de la función verosimilitud. De tal forma que cada iteración asegura una mejora en las estimaciones.
EM algorithm
El algoritmo EM es un método iterativo para encontrar las estimaciones de Máxima Verosimilitud y generalmente se emplea cuando los datos están incompletos o tienen valores faltantes.
El método explota el hecho de que, en general, el problema de maximización es más fácil para los datos completos que los datos incompletos. Cada iteración del EM algoritmo consta de dos pasos:
- Esperanza (paso E) y
- Maximización (paso M).
En el caso que estudiamos, los datos no observados son las etiquetas que identifican a que grupo pertenece cada observación.
Por lo que se introduce un conjunto de variables latentes que indican a qué componente pertenece cada dato. Estas variables toman el valor de 0 cuando la etiqueta no corresponde al grupo al que pertenece el dato y toma el valor de 1 cuando sí corresponde.
Por lo que la verosimilitud completa de los datos se define como:
\[ L_c(\theta) = \Pi_{i=1}^{n}\Pi_{k=1}^{K}(\pi_k \phi_k(x_i|\theta_k ))^{z_{ik}} \]
Donde \(z_{ik}=1\) indica que la observación \(x_i\) se origina de la componente \(k\) y en otro caso \(z_{ik}=0\).
De esta manera, la log verosimilitud de datos completos es
\[ \begin{split} l_c(\theta)&=\log L_c(\theta)\\ &=\sum_{i=1}^{n}\sum_{k=1}^{K}z_{ik}[\log \pi_k + \log (\phi_k(x_i|\theta_k)] \end{split} \]
Paso E (Esperanza)
Calcular la esperanza condicional de la log verosimilitud de datos completos, dados los datos observados, y los estimadores del periodo anterior.
Dado que la verosimilitud de datos completos es lineal en las variables latentes y las esperanzas son lineales, entonces la esperanza de la log verosimilitud de datos completos es obtenida reemplazando las variables latentes por sus valores esperados condicionales en los datos observados y estimaciones de los parámetros de la iteración anterior. Estos valores son
\[ \begin{split} \pi_{ik}^{(s)}&=E[z_{ik}|\theta^{s-1}]\\ &= \frac{\pi_{k}^{(s-1)}\phi_{k}(x_i|\theta_{k}^{(s-1)})}{\sum_{k'=1}^{K}\pi_{k'}^{(s-1)}\phi_{k'}(x_i|\theta_{k'}^{(s-1)})} \end{split} \] Para \(i=1,\ldots,n\) y \(k=1,\ldots,K\), donde \(\pi_{k}^{(s-1)}\) corresponde a la probabilidad posterior de que \(x_i\) provenga de la \(k\) - ésima componente del modelo mixto, calculado en la iteración \(s\).
El resultado es la función \(Q\) dada por:
\[ Q(\theta|\theta^{(s-1)})=\sum_{i=1}^n\sum_{k=1}^{K}\pi_{ik}^{(s)}\Big[\log(\pi_k)+\log(\phi_k(x_i|\theta_k))\Big]. \] Es importante notar que el paso E no depende de la forma funcional de las densidades, sino solo de las variables latentes.
Paso M (Maximización)
Las nuevas estimaciones para \(\pi\) y \(\theta\) son obtenidas maximizando la función \(Q\). El problema de optimización puede ser resuelto de manera separada para \(\pi\) y para cada \(\theta_k\).
La estimación para la actualización de \(\pi\) está dada por
\[ \hat{\pi}_{k}^{(s)}=\frac{1}{n}\sum_{i=1}^{n}\pi_{ik}^{(s)}. \] Y la estimación para las medias y varianzas está dada por:
\[ \hat{\mu}_{k}^{(s)}=\frac{\sum_{i=1}^{n}\pi_{ik}^{(s)}x_i}{\sum_{i=1}^{n}\pi_{ik}^{s}} \] y
\[ (\sigma_{k}^2)^{(s)}=\frac{\sum_{i=1}^{n}\pi_{ik}^{(s)}(x_i-\hat{\mu}_{k}^{(s)})^2}{\sum_{i=1}^{n}\pi_{ik}^{s}} \]
Inicialización
El algoritmo necesita un punto de partida y muchas veces el éxito del modelo depende de una buena inicialización. Aquí presentamos tres posibles formas hacer esto:
- Inicialización Euclidiana basada en la distancia
Se elijen \(k\) elementos de la muestra como centros de los \(k\) clusters y el resto de los elementos en las muestras son asignados al cluster más cercano.
- k - means
El algoritmo k - means puede ser empleado para identificar los clusters. Con base en esa información se calculan las medias y varianzas.
- Inicialización aleatoria
Se crea una partición aleatoria con \(K\) grupos. Y se asigna de manera aleatoria, con la misma probabilidad, cada observación en cada grupo. Es forma de inicialización fue propuesta por McLachlan and Peel (2000).
Implementación
Presentamos el ejemplo de Scrucca et al 2016, en el que se utiliza un modelo Gaussiano mixto para clasificar diagnósticos de cáncer como maligno o benigno. Los datos están disponibles en http://archive.ics.uci.edu/ml/datasets/Breast+Cancer+Wisconsin+(Diagnostic)
La base de datos cuenta con 569 pacientes y 30 características de los núcleos celulares. Para cada paciente, el cáncer fue diagnosticado como maligno o benigno. De las 30 características, solo se consideran 3:
- área extrema.
- suavidad extrema, variación local en longitudes de radio.
- textura media, desviación estándar de los valores de la escala de grises
Cuya estadística descriptiva es la siguiente:
## texture.mean area.extreme smoothness.extreme
## Min. : 9.71 Min. : 185.2 Min. :0.07117
## 1st Qu.:16.17 1st Qu.: 515.3 1st Qu.:0.11660
## Median :18.84 Median : 686.5 Median :0.13130
## Mean :19.29 Mean : 880.6 Mean :0.13237
## 3rd Qu.:21.80 3rd Qu.:1084.0 3rd Qu.:0.14600
## Max. :39.28 Max. :4254.0 Max. :0.22260Asignamos 2/3 de las observaciones al conjunto de entrenamiento y el resto al conjunto de prueba.
## Class.train
## B M
## 251 128## Class.test
## B M
## 106 84\(k\) means
## List of 9
## $ cluster : int [1:569] 1 1 1 2 1 2 1 2 2 2 ...
## $ centers : num [1:2, 1:3] 21.69 18.57 1753.02 619.65 0.14 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:2] "1" "2"
## .. ..$ : chr [1:3] "texture.mean" "area.extreme" "smoothness.extreme"
## $ totss : num 1.84e+08
## $ withinss : num [1:2] 35993442 18610175
## $ tot.withinss: num 54603616
## $ betweenss : num 1.3e+08
## $ size : int [1:2] 131 438
## $ iter : int 1
## $ ifault : int 0
## - attr(*, "class")= chr "kmeans"
Estimamos un primer modelo
## ------------------------------------------------
## Gaussian finite mixture model for classification
## ------------------------------------------------
##
## EDDA model summary:
##
## log-likelihood n df BIC
## -2934.09 379 12 -5939.43
##
## Classes n % Model G
## B 251 66.23 VVI 1
## M 128 33.77 VVI 1
##
## Training confusion matrix:
## Predicted
## Class B M
## B 249 2
## M 17 111
## Classification error = 0.0501
## Brier score = 0.0374
##
## Test confusion matrix:
## Predicted
## Class B M
## B 103 3
## M 5 79
## Classification error = 0.0421
## Brier score = 0.0357Observamos lo siguiente:
Utilizamos el Criterio de Información Bayesiana (BIC), para seleccionar el mejor modelo.
El modelo elegido es el VVI: por lo que cada grupo, pacientes con cáncer benigno y maligno, está descrito por una solo distribución normal con volumen y forma variables, pero ambas alineadas con los ejes de coordenadas.
Hay un error en entrenamiento de 5%.
Hay un error en prueba de 4%.
Y obtenemos un error del 5% con una desviación estándar de 0.01.
Ahora implementamos validación cruzada con n=10 particiones.
## error se
## 0.05540897 0.01068636Ahora vamos a mostrar la clasificación gráficamente. Recuerde que se cuenta con 3 dimensiones, área extrema, suavidad extrema y variación local, por lo que presentaremos tres gráficas en el plano para poder visualizar la clasificación.
A continuación, mostramos el modelo en un espacio de dimensión reducida.
## -----------------------------------------------------------------
## Dimension reduction for model-based clustering and classification
## -----------------------------------------------------------------
##
## Mixture model type: EDDA
##
## Classes n Model G
## B 251 VVI 1
## M 128 VVI 1
##
## Estimated basis vectors:
## Dir1 Dir2
## texture.mean -9.5376e-01 -1.0000e+00
## area.extreme 3.0056e-01 2.8068e-03
## smoothness.extreme -3.8247e-05 -4.2721e-05
##
## Dir1 Dir2
## Eigenvalues 0.59101 0.054645
## Cum. % 91.53656 100.000000
Ahora permitimos que cada cluster se pueda modelar con una o más densidades normales.
## ------------------------------------------------
## Gaussian finite mixture model for classification
## ------------------------------------------------
##
## MclustDA model summary:
##
## log-likelihood n df BIC
## -2893.593 379 26 -5941.562
##
## Classes n % Model G
## B 251 66.23 EEI 3
## M 128 33.77 EVI 2
##
## Training confusion matrix:
## Predicted
## Class B M
## B 248 3
## M 8 120
## Classification error = 0.029
## Brier score = 0.0262
##
## Test confusion matrix:
## Predicted
## Class B M
## B 103 3
## M 5 79
## Classification error = 0.0421
## Brier score = 0.0273## error se
## 0.031662269 0.008463335
## -----------------------------------------------------------------
## Dimension reduction for model-based clustering and classification
## -----------------------------------------------------------------
##
## Mixture model type: MclustDA
##
## Classes n Model G
## B 251 EEI 3
## M 128 EVI 2
##
## Estimated basis vectors:
## Dir1 Dir2
## texture.mean 9.8998e-01 -9.9999e-01
## area.extreme 1.4123e-01 3.4221e-03
## smoothness.extreme 4.3501e-05 -4.2716e-05
##
## Dir1 Dir2
## Eigenvalues 0.67062 0.20981
## Cum. % 76.16911 100.00000
## ------------------------------------------------
## Gaussian finite mixture model for classification
## ------------------------------------------------
##
## MclustDA model summary:
##
## log-likelihood n df BIC
## -2910.808 379 26 -5975.992
##
## Classes n % Model G
## B 251 66.23 EEE 2
## M 128 33.77 EEE 2
##
## Training confusion matrix:
## Predicted
## Class B M
## B 247 4
## M 10 118
## Classification error = 0.0369
## Brier score = 0.0297
##
## Test confusion matrix:
## Predicted
## Class B M
## B 104 2
## M 4 80
## Classification error = 0.0316
## Brier score = 0.0248Conclusiones
Los modelos GMM y el algoritmo EM son herramientas poderosas de aprendizaje no supervisado que, permiten implementar clasificadores. En el caso presentado en este trabajo, como era de esperarse, las observaciones que no se diferenciaban bien de un cluster u otro, fueron las que pueden hacer la diferencia en el mejor rendimiento de los modelos GMM, respecto a \(k\) means. También es importante mencionar que el algoritmo EM permite asegurar un avance siempre en cada iteración y que la inicialización del modelo es muy importante para alcanzar máximos globales.
Referencias
[1] Scrucca, Luca & Fop, Michael & Murphy, Thomas & Raftery, Adrian. (2016). mclust 5: Clustering, Classification and Density Estimation Using Gaussian Finite Mixture Models. The R Journal. 8. 205-233. 10.32614/RJ-2016-021.
[2] Tatjana Miljkovic, Bettina Grün, Modeling loss data using mixtures of distributions, Insurance: Mathematics and Economics, Volume 70, 2016, Pages 387-396,
[3] Ramesh Sridharan, Gaussian mixture models and the EM algorithm, https://people.csail.mit.edu/rameshvs/content/gmm-em.pdf
[4] Arnaud Mignan, Mixture modelling from scratch, in R, May 2019.
[5] Breast Cancer Wisconsin (Diagnostic) Data Set, http://archive.ics.uci.edu/ml/datasets/Breast+Cancer+Wisconsin+(Diagnostic)
[6] McLachlan, S., Peel, D., 2000. Finite Mixture Models. John Wiley & Sons, Hobuken, NJ.