Die latente Profilanalyse (LPA) versucht, Cluster von Individuen (dh latente Profile) basierend auf Antworten auf eine Reihe kontinuierlicher Variablen (dh Indikatoren) zu identifizieren. LPA geht davon aus, dass es unbeobachtete latente Profile gibt, die Antwortmuster auf Indikatorelemente erzeugen.
Hier werde ich ein kurzes Beispiel von LPA durchgehen, um Gruppen von Menschen basierend auf ihren Interessen / Hobbys zu identifizieren. Die Daten stammen aus der Young People Survey, frei verfügbar auf Kaggle.com.
Hier ist ein kleiner Vorgeschmack auf das, was wir wollen:
Terminologie Hinweis: Die Begriffe Cluster, Profile, Klassen und Gruppen werden synonym verwendet, aber es gibt feine Unterschiede. Ich werde mich hauptsächlich an das Profil halten, um mich auf eine Gruppierung von Fällen zu beziehen, in Übereinstimmung mit der LPA-Terminologie. Wir sollten beachten, dass LPA ein Zweig der Gaußschen Finite-Mix-Modellierung ist, der die latente Klassenanalyse (LCA) umfasst. Der Unterschied zwischen LPA und LCA ist konzeptionell, nicht rechnerisch: LPA verwendet kontinuierliche Indikatoren und LCA verwendet binäre Indikatoren. LPA ist ein probabilistisches Modell, dh es modelliert die Wahrscheinlichkeit, dass ein Fall zu einem Profil gehört. Dies ist einem Ansatz wie K-means überlegen, der Entfernungsalgorithmen verwendet.
Abgesehen davon laden wir die Daten ein.
library(tidyverse)survey <- read_csv("https://raw.githubusercontent.com/whipson/tidytuesday/master/young_people.csv") %>% select(History:Pets)Die Daten beziehen sich auf 32 Interessen /Hobbys. Jeder Artikel wird von 1 (nicht interessiert) bis 5 (sehr interessiert) eingestuft.
Die Beschreibung auf Kaggle deutet darauf hin, dass es unvorsichtige Antworten geben kann (z. B. Teilnehmer, die immer wieder denselben Wert ausgewählt haben). Wir können das sorglose Paket verwenden, um „string responding“ zu identifizieren. Suchen wir auch nach multivariaten Ausreißern mit Mahalanobis-Entfernung (siehe meinen vorherigen Beitrag zu Mahalanobis zur Identifizierung von Ausreißern).
library(careless)library(psych)interests <- survey %>% mutate(string = longstring(.)) %>% mutate(md = outlier(., plot = FALSE))Wir begrenzen die Zeichenfolge auf maximal 10 und verwenden einen Mahalanobis D-Cutoff von alpha = .001.
cutoff <- (qchisq(p = 1 - .001, df = ncol(interests)))interests_clean <- interests %>% filter(string <= 10, md < cutoff) %>% select(-string, -md)Das Paket mclust führt verschiedene Arten von modellbasiertem Clustering und Dimensionsreduktion durch. Plus, es ist wirklich intuitiv zu bedienen. Es erfordert vollständige Daten (keine fehlenden), daher entfernen wir in diesem Beispiel Fälle mit NAs. Dies ist nicht der bevorzugte Ansatz; Wir wären besser dran, zu unterstellen. Aber zur Veranschaulichung funktioniert das gut. Ich werde auch alle Indikatoren standardisieren, damit beim Zeichnen der Profile die Unterschiede zwischen den Clustern klarer erkennbar sind. Das Ausführen dieses Codes dauert einige Minuten.
library(mclust)interests_clustering <- interests_clean %>% na.omit() %>% mutate_all(list(scale))BIC <- mclustBIC(interests_clustering)Wir beginnen mit der Darstellung Bayes’scher Informationskriterien für alle Modelle mit Profilen von 1 bis 9.
plot(BIC)
Es ist nicht sofort klar, welches Modell das Beste ist, da die y-Achse so groß ist und viele der Modelle nahe beieinander liegen. zusammenfassung(BIC) zeigt die drei besten Modelle basierend auf 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.461Der höchste BIC kommt von VVE, 3. Dies besagt, dass es 3 Cluster mit variablem Volumen, variabler Form, gleicher Ausrichtung und ellipsodialer Verteilung gibt (siehe Abbildung 2 aus diesem Papier für eine visuelle Darstellung). VEE, 3 ist jedoch nicht weit dahinter und kann tatsächlich ein theoretisch nützlicheres Modell sein, da es die Form der Verteilung auf gleich beschränkt. Aus diesem Grund gehen wir mit VEE, 3.
Wenn wir dieses Modell genauer betrachten wollen, speichern wir es als Objekt und untersuchen es mit 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 210Die Ausgabe beschreibt die geometrischen Eigenschaften der Profile und die Anzahl der Fälle, die in jedem der drei Cluster klassifiziert sind.
BIC ist einer der am besten geeigneten Indizes, aber es wird immer empfohlen, nach weiteren Beweisen zu suchen, dass die von uns gewählte Lösung die richtige ist. Wir können auch Werte des Kriteriums Integrated Completed Likelikood (ICL) vergleichen. Weitere Informationen finden Sie in diesem Dokument. ICL unterscheidet sich nicht wesentlich von BIC, außer dass es Lösungen mit größerer Entropie- oder Klassifikationsunsicherheit bestraft.
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.203Wir sehen ähnliche Ergebnisse. ICL schlägt vor, dass das Modell VEE, 3 recht gut passt. Schließlich führen wir den Bootstrap-Likelihood-Ratio-Test (BLRT) durch, der die Modellanpassung zwischen k-1- und k-Cluster-Modellen vergleicht. Mit anderen Worten, es wird geprüft, ob eine Erhöhung der Profile die Passform erhöht. Basierend auf Simulationen von Nylund, Asparouhov und Muthén (2007) sind BIC und BLRT die besten Indikatoren dafür, wie viele Profile es gibt. Wenn Sie also nur folgen, schlage ich vor, sie zu überspringen, es sei denn, Sie möchten eine Kaffeepause einlegen.
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 schlägt auch vor, dass eine 3-Profil-Lösung ideal ist.
LPA visualisieren
Nachdem wir uns nun für eine 3-Profil-Lösung entschieden haben, zeichnen wir die Ergebnisse auf. Insbesondere möchten wir sehen, wie sich die Profile bei den Indikatoren unterscheiden, dh bei den Elementen, aus denen die Profile bestehen. Wenn die Lösung theoretisch sinnvoll ist, sollten wir Unterschiede sehen, die Sinn machen.
Zuerst extrahieren wir die Mittelwerte für jedes Profil (denken Sie daran, wir haben diese als standardisiert ausgewählt). Dann verwenden wir pivot_longer , um es in eine lange Form zu bringen. Beachten Sie, dass ich Werte über +1 SD beschneide, da sonst Probleme beim Zeichnen auftreten.
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))Hier ist der Code für den Plot. Ich ordne die Indikatoren so an, dass ähnliche Aktivitäten nahe beieinander liegen.
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")
Wir haben viele Indikatoren (mehr als typisch für LPA), aber wir sehen einige interessante Unterschiede. Offensichtlich ist die rote Gruppe an Wissenschaft interessiert und die blaue Gruppe zeigt ein größeres Interesse an Kunst und Geisteswissenschaften. Die grüne Fraktion scheint sowohl an Wissenschaft als auch an Kunst desinteressiert, aber mäßig an anderen Dingen interessiert zu sein.
Wir können dieses Diagramm informativer gestalten, indem wir Profilnamen und Proportionen einfügen. Ich werde diesen Plot auch als Objekt speichern, damit wir etwas wirklich Cooles damit machen können!
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
Das wirklich Coole, das ich machen möchte, ist eine interaktive Handlung. Warum sollte ich das tun wollen? Nun, eines der Probleme mit dem statischen Diagramm ist, dass es bei so vielen Indikatoren schwierig ist, die Werte für jeden Indikator zu lesen. Ein interaktives Diagramm ermöglicht es dem Leser, bestimmte Indikatoren oder Profile von Interesse einzugrenzen. Wir verwenden plotly , um unseren statischen Plot in einen interaktiven zu verwandeln.
library(plotly)ggplotly(p, tooltip = c("Interest", "Mean")) %>% layout(legend = list(orientation = "h", y = 1.2))Es gibt ein kurzes Beispiel für LPA. Insgesamt denke ich, dass LPA ein großartiges Werkzeug für die explorative Analyse ist, obwohl ich seine Reproduzierbarkeit in Frage stelle. Wichtig ist, dass der Statistiker bei der Entscheidung über die Anzahl der Profile sowohl Anpassungsindizes als auch Theorie berücksichtigt.
Referenzen & Ressourcen
Bertoletti, M., Friel, N., & Rastelli, R. (2015). Auswahl der Anzahl der Cluster in einem Finite-Mix-Modell unter Verwendung eines exakten integrierten Wahrscheinlichkeitskriteriums. https://arxiv.org/pdf/1411.4257.pdf.
Nylund, K. L., Asparouhov, T., & Muthén, B. O. (2007). Entscheidung über die Anzahl der Klassen in Latent Class Analysis und Growth Mixture Modeling: Eine Monte-Carlo-Simulationsstudie. Strukturgleichungsmodellierung, 14, 535-569.
Scrucca, L., Fop, M., Murphy, T. B., & Raftery, A. E. (2016). mclust5: Clustering, Klassifizierung und Dichteschätzung mit Gaußschen Finite-Mix-Modellen. Das R Journal, 8, 289-317.