El Análisis de Perfil Latente (LPA) intenta identificar grupos de individuos (es decir, perfiles latentes) basados en respuestas a una serie de variables continuas (es decir, indicadores). LPA asume que hay perfiles latentes no observados que generan patrones de respuestas en elementos de indicadores.
Aquí, pasaré por un ejemplo rápido de LPA para identificar grupos de personas en función de sus intereses/pasatiempos. Los datos provienen de la Encuesta de Jóvenes, disponible gratuitamente en Kaggle.com.
Este es un adelanto de lo que buscamos:
Nota terminológica: Las personas usan los términos clusters, perfiles, clases y grupos indistintamente, pero hay diferencias sutiles. Me ceñiré principalmente al perfil para referirme a una agrupación de casos, de acuerdo con la terminología de LPA. Debemos tener en cuenta que el LPA es una rama del Modelado de Mezclas Finitas Gaussianas, que incluye el Análisis de Clases Latentes (LCA). La diferencia entre LPA y LCA es conceptual, no computacional: LPA utiliza indicadores continuos y LCA utiliza indicadores binarios. LPA es un modelo probabilístico, lo que significa que modela la probabilidad de que el caso pertenezca a un perfil. Esto es superior a un enfoque como K-means que utiliza algoritmos de distancia.
Con eso a un lado, carguemos los datos.
library(tidyverse)survey <- read_csv("https://raw.githubusercontent.com/whipson/tidytuesday/master/young_people.csv") %>% select(History:Pets)Los datos es de 32 intereses/aficiones. Cada artículo está clasificado del 1 (no está interesado) al 5 (muy interesado).
La descripción en Kaggle sugiere que puede haber respuestas descuidadas (por ejemplo, participantes que seleccionaron el mismo valor una y otra vez). Podemos usar el paquete descuidado para identificar «cadena que responde». También busquemos valores atípicos multivariantes con Distancia de Mahalanobis (vea mi publicación anterior en Mahalanobis para identificar valores atípicos).
library(careless)library(psych)interests <- survey %>% mutate(string = longstring(.)) %>% mutate(md = outlier(., plot = FALSE))Limitaremos la cadena que responde a un máximo de 10 y usaremos un corte D Mahalanobis de alpha = .001.
cutoff <- (qchisq(p = 1 - .001, df = ncol(interests)))interests_clean <- interests %>% filter(string <= 10, md < cutoff) %>% select(-string, -md)El paquete mclust realiza varios tipos de agrupamiento basado en modelos y reducción de dimensiones. Además, es muy intuitivo de usar. Requiere datos completos (no faltan), por lo que en este ejemplo eliminaremos los casos con NAs. Este no es el enfoque preferido; sería mejor imputar. Pero para fines ilustrativos, esto funciona bien. También voy a estandarizar todos los indicadores para que cuando trazemos los perfiles sea más claro ver las diferencias entre los grupos. Ejecutar este código llevará unos minutos.
library(mclust)interests_clustering <- interests_clean %>% na.omit() %>% mutate_all(list(scale))BIC <- mclustBIC(interests_clustering)Comenzaremos trazando Criterios de Información Bayesianos para todos los modelos con perfiles que van del 1 al 9.
plot(BIC)
No está claro de inmediato qué modelo es el mejor, ya que el eje y es tan grande y muchos de los modelos puntúan juntos. resumen (BIC) muestra los tres modelos principales basados en BIC.
summary(BIC)## Best BIC values:## VVE,3 VEE,3 EVE,3## BIC -75042.7 -75165.1484 -75179.165## BIC diff 0.0 -122.4442 -136.461El BIC más alto proviene de VVE, 3. Esto dice que hay 3 grupos con volumen variable, forma variable, orientación igual y distribución elipsodial (vea la Figura 2 de este artículo para una imagen visual). Sin embargo, VEE, 3 no se queda atrás y en realidad puede ser un modelo más útil teóricamente, ya que restringe la forma de la distribución para que sea igual. Por esta razón, vamos a ir con VEE, 3.
Si queremos ver este modelo más de cerca, lo guardamos como un objeto e inspeccionamos con summary().
mod1 <- Mclust(interests_clustering, modelNames = "VEE", G = 3, x = BIC)summary(mod1)## ---------------------------------------------------- ## Gaussian finite mixture model fitted by EM algorithm ## ---------------------------------------------------- ## ## Mclust VEE (ellipsoidal, equal shape and orientation) model with 3 components: ## ## log-likelihood n df BIC ICL## -35455.83 874 628 -75165.15 -75216.14## ## Clustering table:## 1 2 3 ## 137 527 210La salida se describen las características geométricas de los perfiles y el número de casos clasificados en cada uno de los tres grupos.
El BIC es uno de los mejores índices de ajuste, pero siempre se recomienda buscar más pruebas de que la solución que hemos elegido es la correcta. También podemos comparar los valores del criterio Integrated Completed Likelikood (ICL). Consulte este documento para obtener más detalles. ICL no es muy diferente de BIC, excepto que agrega una penalización en soluciones con mayor entropía o incertidumbre de clasificación.
ICL <- mclustICL(interests_clustering)plot(ICL)
summary(ICL)## Best ICL values:## VVE,3 VEE,3 EVE,3## ICL -75134.69 -75216.13551 -75272.891## ICL diff 0.00 -81.44795 -138.203vemos resultados similares. ICL sugiere que el modelo VEE, 3 encaja bastante bien. Finalmente, realizaremos la Prueba de Razón de Verosimilitud de Arranque (BLRT) que compara el ajuste del modelo entre los modelos de clúster k-1 y k. En otras palabras, busca ver si un aumento en los perfiles aumenta el ajuste. Basado en simulaciones de Nylund, Asparouhov y Muthén (2007), BIC y BLRT son los mejores indicadores de cuántos perfiles hay. Esta línea de código tardará mucho tiempo en ejecutarse, por lo que si solo estás siguiendo, te sugiero que la omitas a menos que quieras salir a tomar un café.
mclustBootstrapLRT(interests_clustering, modelName = "VEE")## Warning in mclustBootstrapLRT(interests_clustering, modelName = "VEE"): some## model(s) could not be fitted!## ------------------------------------------------------------- ## Bootstrap sequential LRT for the number of mixture components ## ------------------------------------------------------------- ## Model = VEE ## Replications = 999 ## LRTS bootstrap p-value## 1 vs 2 197.0384 0.001## 2 vs 3 684.8743 0.001## 3 vs 4 -124.1935 1.000BLRT también sugiere que un 3-perfil solución es ideal.
Visualización de LPA
Ahora que confiamos en nuestra elección de una solución de 3 perfiles, vamos a trazar los resultados. En concreto, queremos ver cómo difieren los perfiles en los indicadores, es decir, los elementos que componen los perfiles. Si la solución es teóricamente significativa, deberíamos ver diferencias que tengan sentido.
Primero, extraeremos los medios para cada perfil (recuerde, elegimos estos para estandarizarlos). Luego, usamos pivot_longer para organizarlo en forma larga. Tenga en cuenta que estoy recortando valores que exceden +1 SD, de lo contrario, nos encontraremos con problemas de trazado.
means <- data.frame(mod1$parameters$mean) %>% rownames_to_column() %>% rename(Interest = rowname) %>% pivot_longer(cols = c(X1, X2, X3), names_to = "Profile", values_to = "Mean") %>% mutate(Mean = round(Mean, 2), Mean = ifelse(Mean > 1, 1, Mean))Aquí está el código para la gráfica. Estoy reordenando los indicadores para que las actividades similares estén juntas.
means %>% ggplot(aes(Interest, Mean, group = Profile, color = Profile)) + geom_point(size = 2.25) + geom_line(size = 1.25) + scale_x_discrete(limits = c("Active sport", "Adrenaline sports", "Passive sport", "Countryside, outdoors", "Gardening", "Cars", "Art exhibitions", "Dancing", "Musical instruments", "Theatre", "Writing", "Reading", "Geography", "History", "Law", "Politics", "Psychology", "Religion", "Foreign languages", "Biology", "Chemistry", "Mathematics", "Medicine", "Physics", "Science and technology", "Internet", "PC", "Celebrities", "Economy Management", "Fun with friends", "Shopping", "Pets")) + labs(x = NULL, y = "Standardized mean interest") + theme_bw(base_size = 14) + theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "top")
Tenemos muchos indicadores (más de los típicos para LPA), pero vemos algunas diferencias interesantes. Claramente, el grupo rojo está interesado en la ciencia y el grupo azul muestra un mayor interés en las artes y las humanidades. El grupo de los Verdes parece desinteresado tanto por la ciencia como por el arte, pero moderadamente interesado en otras cosas.
Podemos hacer que esta gráfica sea más informativa conectando nombres y proporciones de perfil. También voy a guardar esta trama como un objeto para que podamos hacer algo realmente genial con él!
p <- means %>% mutate(Profile = recode(Profile, X1 = "Science: 16%", X2 = "Disinterest: 60%", X3 = "Arts & Humanities: 24%")) %>% ggplot(aes(Interest, Mean, group = Profile, color = Profile)) + geom_point(size = 2.25) + geom_line(size = 1.25) + scale_x_discrete(limits = c("Active sport", "Adrenaline sports", "Passive sport", "Countryside, outdoors", "Gardening", "Cars", "Art exhibitions", "Dancing", "Musical instruments", "Theatre", "Writing", "Reading", "Geography", "History", "Law", "Politics", "Psychology", "Religion", "Foreign languages", "Biology", "Chemistry", "Mathematics", "Medicine", "Physics", "Science and technology", "Internet", "PC", "Celebrities", "Economy Management", "Fun with friends", "Shopping", "Pets")) + labs(x = NULL, y = "Standardized mean interest") + theme_bw(base_size = 14) + theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "top")p
algo realmente bueno que quiero hacer es hacer un interactivo de la parcela. ¿Por qué querría hacer esto? Bueno, uno de los problemas con la gráfica estática es que con tantos indicadores es difícil leer los valores de cada indicador. Una trama interactiva permite al lector concentrarse en indicadores específicos o perfiles de interés. Usaremos plotly para convertir nuestra trama estática en una interactiva.
library(plotly)ggplotly(p, tooltip = c("Interest", "Mean")) %>% layout(legend = list(orientation = "h", y = 1.2))Hay un ejemplo rápido de la LPA. En general, creo que el LPA es una gran herramienta para el Análisis Exploratorio, aunque cuestiono su reproducibilidad. Lo importante es que el estadístico tenga en cuenta tanto los índices de ajuste como la teoría al decidir el número de perfiles.
Referencias & Recursos
Bertoletti, M., Friel, N., & Rastelli, R. (2015). Elegir el número de clústeres en un modelo de mezcla finita utilizando un criterio de Verosimilitud Completa Integrada exacta. https://arxiv.org/pdf/1411.4257.pdf.
Nylund, K. L., Asparouhov, T., & Muthén, B. O. (2007). Deciding on the Number of Classes in Latent Class Analysis and Growth Mixture Modeling: A Monte Carlo Simulation Study. Modelado de ecuaciones estructurales, 14, 535-569.Scrucca, L., Fop, M., Murphy, T. B., & Raftery, A. E. (2016). mclust5: Clustering, Classification and Density Estimation Using Gaussian Finite Mixture Models (en inglés). The R Journal, 8, 289-317.