Latent Profile Analysis (LPA) cerca di identificare gruppi di individui (cioè profili latenti) in base alle risposte a una serie di variabili continue (cioè indicatori). LPA presuppone che ci siano profili latenti non osservati che generano modelli di risposte sugli elementi dell’indicatore.
Qui, passerò attraverso un rapido esempio di LPA per identificare gruppi di persone in base ai loro interessi/hobby. I dati provengono dall’Indagine Giovani, disponibile gratuitamente su Kaggle.com.
Ecco una sbirciatina a quello che stiamo andando per:
Terminologia nota: Le persone usano i termini cluster, profili, classi e gruppi in modo intercambiabile, ma ci sono sottili differenze. Mi limiterò principalmente al profilo per fare riferimento a un raggruppamento di casi, in linea con la terminologia LPA. Dovremmo notare che LPA è un ramo della modellazione gaussiana della miscela finita, che include l’analisi di classe latente (LCA). La differenza tra LPA e LCA è concettuale, non computazionale: LPA utilizza indicatori continui e LCA utilizza indicatori binari. LPA è un modello probabilistico, il che significa che modella la probabilità di caso appartenente a un profilo. Questo è superiore a un approccio come K-means che utilizza algoritmi di distanza.
A parte questo, carichiamo i dati.
library(tidyverse)survey <- read_csv("https://raw.githubusercontent.com/whipson/tidytuesday/master/young_people.csv") %>% select(History:Pets)I dati sono su 32 interessi / hobby. Ogni articolo è classificato da 1 (non interessato) a 5 (molto interessato).
La descrizione su Kaggle suggerisce che potrebbe esserci una risposta incurante (ad esempio, i partecipanti che hanno selezionato lo stesso valore più e più volte). Possiamo usare il pacchetto careless per identificare “string responding”. Cerchiamo anche valori anomali multivariati con Mahalanobis Distanza (vedere il mio precedente post su Mahalanobis per identificare valori anomali).
library(careless)library(psych)interests <- survey %>% mutate(string = longstring(.)) %>% mutate(md = outlier(., plot = FALSE))Cap string risponde ad un massimo di 10 e utilizzare un Mahalanobis D cutoff di alpha = .001.
cutoff <- (qchisq(p = 1 - .001, df = ncol(interests)))interests_clean <- interests %>% filter(string <= 10, md < cutoff) %>% select(-string, -md)Il pacchetto mclust esegue vari tipi di clustering basato su modelli e riduzione delle dimensioni. Inoltre, è davvero intuitivo da usare. Richiede dati completi (non mancanti), quindi per questo esempio rimuoveremo i casi con NAs. Questo non è l’approccio preferito; sarebbe meglio imputare. Ma a scopo illustrativo, questo funziona bene. Ho anche intenzione di standardizzare tutti gli indicatori in modo che quando tracciamo i profili è più chiaro per vedere le differenze tra i cluster. L’esecuzione di questo codice richiederà alcuni minuti.
library(mclust)interests_clustering <- interests_clean %>% na.omit() %>% mutate_all(list(scale))BIC <- mclustBIC(interests_clustering)Inizieremo tracciando Criteri informativi bayesiani per tutti i modelli con profili che vanno da 1 a 9.
plot(BIC)
Non è immediatamente chiaro quale modello sia il migliore poiché l’asse y è così grande e molti dei modelli sono vicini. riepilogo (BIC) mostra i primi tre modelli basati su 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.461Il BIC più alto proviene da VVE, 3. Questo dice che ci sono 3 cluster con volume variabile, forma variabile, uguale orientamento e distribuzione ellissodiale (vedere la Figura 2 di questo documento per un visual). Tuttavia, VEE, 3 non è molto indietro e in realtà potrebbe essere un modello più teoricamente utile poiché limita la forma della distribuzione ad essere uguale. Per questo motivo, andremo con VEE, 3.
Se vogliamo guardare questo modello più da vicino, lo salviamo come oggetto e lo ispezioniamo 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 210L’output descrive le caratteristiche geometriche dei profili e il numero di casi classificati in ciascuno dei tre cluster.
BIC è uno dei migliori indici di adattamento, ma è sempre consigliabile cercare ulteriori prove che la soluzione che abbiamo scelto sia quella corretta. Possiamo anche confrontare i valori del criterio Integrated Completed Likelikood (ICL). Vedere questo documento per maggiori dettagli. ICL non è molto diverso da BIC, tranne che aggiunge una penalità su soluzioni con maggiore entropia o incertezza di classificazione.
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.203vediamo risultati simili. ICL suggerisce che il modello VEE, 3 si adatta abbastanza bene. Infine, eseguiremo il Bootstrap Likelihood Ratio Test (BLRT) che confronta l’adattamento del modello tra i modelli k-1 e k cluster. In altre parole, sembra vedere se un aumento dei profili aumenta in forma. Sulla base di simulazioni di Nylund, Asparouhov e Muthén (2007) BIC e BLRT sono i migliori indicatori per quanti profili ci sono. Questa riga di codice richiederà molto tempo per essere eseguita, quindi se stai solo seguendo ti suggerisco di saltarla a meno che tu non voglia uscire per una pausa caffè.
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 suggerisce anche che una soluzione a 3 profili è ideale.
Visualizzazione LPA
Ora che siamo fiduciosi nella nostra scelta di una soluzione a 3 profili, tracciamo i risultati. Nello specifico, vogliamo vedere come i profili differiscono sugli indicatori, cioè gli elementi che componevano i profili. Se la soluzione è teoricamente significativa, dovremmo vedere differenze che hanno senso.
Per prima cosa, estrarremo i mezzi per ogni profilo (ricorda, abbiamo scelto questi per essere standardizzati). Quindi, usiamo pivot_longer per disporlo in forma lunga. Si noti che sto tagliando valori superiori a +1 SD, altrimenti ci imbattiamo in problemi di stampa.
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))Ecco il codice per la trama. Sto riordinando gli indicatori in modo che attività simili siano vicine tra loro.
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")
Abbiamo molti indicatori (più che tipici per LPA), ma vediamo alcune differenze interessanti. Chiaramente il gruppo rosso è interessato alla scienza e il gruppo blu mostra un maggiore interesse per le arti e le discipline umanistiche. Il gruppo verde sembra disinteressato sia alla scienza che all’arte, ma moderatamente interessato ad altre cose.
Possiamo rendere questa trama più informativa inserendo i nomi e le proporzioni dei profili. Ho anche intenzione di salvare questa trama come un oggetto in modo che possiamo fare qualcosa di veramente cool con esso!
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
Il qualcosa di veramente bello che voglio fare è creare una trama interattiva. Perche ‘ dovrei volerlo fare? Bene, uno dei problemi con la trama statica è che con così tanti indicatori è difficile leggere i valori per ciascun indicatore. Una trama interattiva consente al lettore di restringere su specifici indicatori o profili di interesse. Useremo plotly per trasformare la nostra trama statica in una interattiva.
library(plotly)ggplotly(p, tooltip = c("Interest", "Mean")) %>% layout(legend = list(orientation = "h", y = 1.2))C’è un rapido esempio di LPA. Nel complesso, penso che LPA sia un ottimo strumento per l’analisi esplorativa, anche se ne dubito della riproducibilità. Ciò che è importante è che lo statistico consideri sia gli indici di adattamento che la teoria al momento di decidere il numero di profili.
Riferimenti& Risorse
Bertoletti, M., Friel, N.,& Rastelli, R. (2015). Scelta del numero di cluster in un modello di miscela finita utilizzando un criterio esatto di verosimiglianza integrata completata. https://arxiv.org/pdf/1411.4257.pdf.
Nylund, KL, Asparouhov, T.,& Muthén, B. O. (2007). Decidere il numero di classi in Latent Class Analysis e Growth Mixture Modeling: uno studio di simulazione Monte Carlo. Modellazione delle equazioni strutturali, 14, 535-569.
Scrucca, L., Fop, M., Murphy, T. B.,& Raftery, A. E. (2016). mclust5: Clustering, classificazione e stima della densità utilizzando modelli gaussiani a miscela finita. Il giornale R, 8, 289-317.