Search

latente profielanalyse (LPA) probeert clusters van individuen (d.w.z. latente profielen) te identificeren op basis van reacties op een reeks continue variabelen (d.w.z. indicatoren). LPA gaat ervan uit dat er onopgemerkte latente profielen zijn die patronen van reacties op indicatoritems genereren.

Hier ga ik door een snel voorbeeld van LPA om groepen mensen te identificeren op basis van hun interesses/hobby ‘ s. De gegevens zijn afkomstig van de Jongerenenquête, die vrij beschikbaar is op Kaggle.com.

Hier is een voorproefje van wat we gaan doen:

terminologie opmerking: mensen gebruiken de termen clusters, profielen, klassen en groepen door elkaar, maar er zijn subtiele verschillen. Ik zal me vooral houden aan het profiel om te verwijzen naar een groep van gevallen, in overeenstemming met LPA terminologie. We moeten opmerken dat LPA is een tak van Gaussiaanse eindige mengsel modellering, die latente Klasse analyse (LCA) omvat. Het verschil tussen LPA en LCA is conceptueel, niet computationeel: LPA gebruikt continue indicatoren en LCA gebruikt binaire indicatoren. LPA is een probabilistisch model, wat betekent dat het model de waarschijnlijkheid van geval behorend tot een profiel. Dit is superieur aan een benadering zoals K-betekent dat gebruik maakt van afstand algoritmen.

met dat terzijde, laten we de gegevens laden.

library(tidyverse)

survey <- read_csv("https://raw.githubusercontent.com/whipson/tidytuesday/master/young_people.csv") %>% select(History:Pets)

de gegevens hebben betrekking op 32 interesses/hobby ‘ s. Elk item is gerangschikt 1 (niet geïnteresseerd) tot 5 (zeer geïnteresseerd).

de beschrijving op Kaggle suggereert dat er onzorgvuldig kan worden gereageerd (bijvoorbeeld deelnemers die dezelfde waarde steeds opnieuw selecteerden). We kunnen het onzorgvuldige pakket gebruiken om te identificeren “string reageren”. Laten we ook kijken naar multivariate uitschieters met Mahalanobis afstand (zie mijn vorige post op Mahalanobis voor het identificeren van uitschieters).

library(careless)library(psych)interests <- survey %>% mutate(string = longstring(.)) %>% mutate(md = outlier(., plot = FALSE))

we cap tekenreeks die reageert tot een maximum van 10 en gebruiken een Mahalanobis d cutoff van alpha = .001.

cutoff <- (qchisq(p = 1 - .001, df = ncol(interests)))interests_clean <- interests %>% filter(string <= 10, md < cutoff) %>% select(-string, -md)

het pakket mclust voert verschillende soorten modelgebaseerde clustering en dimensiereductie uit. Plus, het is echt intuïtief te gebruiken. Het vereist volledige gegevens (geen ontbrekende), dus voor dit voorbeeld zullen we gevallen verwijderen met NAs. Dit is niet de voorkeur aanpak; we zouden beter af toerekenen. Maar voor illustratieve doeleinden werkt dit prima. Ik ga ook alle indicatoren standaardiseren, zodat wanneer we de profielen plotten, het duidelijker is om de verschillen tussen clusters te zien. Het uitvoeren van deze code duurt een paar minuten.

library(mclust)interests_clustering <- interests_clean %>% na.omit() %>% mutate_all(list(scale))BIC <- mclustBIC(interests_clustering)

we beginnen met het plotten van Bayesiaanse informatiecriteria voor alle modellen met profielen variërend van 1 tot 9.

plot(BIC)

Het is niet direct duidelijk welk model het beste is omdat de y-as zo groot is en veel modellen dicht bij elkaar scoren. samenvatting (BIC) toont de top drie modellen gebaseerd op 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.461

de hoogste BIC komt uit VVE, 3. Dit zegt dat er 3 clusters zijn met variabel volume, variabele vorm, gelijke oriëntatie en ellipsodiale verdeling (zie Figuur 2 van dit papier voor een visuele). Echter, VEE, 3 is niet ver achter en eigenlijk kan een meer theoretisch nuttig model, omdat het beperkt de vorm van de verdeling gelijk te zijn. Daarom gaan we voor VEE, 3.

als we dit model nader willen bekijken, slaan we het op als een object en inspecteren het met 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 210

De output beschrijft de geometrische kenmerken van de profielen en het aantal gevallen dat in elk van de drie clusters is ingedeeld.

BIC is een van de best passende indices, maar het wordt altijd aanbevolen om meer bewijs te zoeken dat de oplossing die we hebben gekozen de juiste is. We kunnen ook de waarden van het geïntegreerde voltooide Likelikood (ICL) – criterium vergelijken. Zie dit document voor meer details. ICL is niet veel anders dan BIC, behalve dat het voegt een straf op oplossingen met een grotere entropie of classificatie onzekerheid.

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

we zien vergelijkbare resultaten. ICL suggereert dat model VEE, 3 vrij goed past. Tot slot zullen we de Bootstrap Likelihood Ratio Test (BLRT) uitvoeren die model fit vergelijkt tussen K-1 en K cluster modellen. Met andere woorden, het ziet ernaar uit om te zien of een toename van profielen toeneemt passen. Gebaseerd op simulaties door Nylund, Asparouhov, and Muthén (2007) BIC en BLRT zijn de beste indicatoren voor hoeveel profielen er zijn. Deze regel van de code zal een lange tijd duren om te draaien, dus als je gewoon volgen langs Ik stel voor overslaan, tenzij je wilt stap uit voor een koffiepauze.

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

BLRT suggereert ook dat een oplossing met 3 profielen ideaal is.

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))

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")

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

library(plotly)ggplotly(p, tooltip = c("Interest", "Mean")) %>% layout(legend = list(orientation = "h", y = 1.2))