Søg

Latent profilanalyse (LPA) forsøger at identificere klynger af individer (dvs.latente profiler) baseret på svar på en række kontinuerlige variabler (dvs. indikatorer). LPA antager, at der er ikke-observerede latente profiler, der genererer mønstre af svar på indikatorelementer.

Her vil jeg gennemgå et hurtigt eksempel på LPA for at identificere grupper af mennesker baseret på deres interesser/hobbyer. Dataene kommer fra Ungdomsundersøgelsen, tilgængelig frit på Kaggle.com.

Her er et smugkig på, hvad vi går efter:

terminologi Bemærk: folk bruger udtrykkene klynger, profiler, klasser og grupper om hverandre, men der er subtile forskelle. Jeg vil for det meste holde mig til profilen for at henvise til en gruppering af sager, i overensstemmelse med LPA-terminologi. Vi skal bemærke, at LPA er en gren af Gaussisk endelig Blandingsmodellering, som inkluderer Latent klasseanalyse (LCA). Forskellen mellem LPA og LCA er konceptuel, ikke beregningsmæssig: LPA bruger kontinuerlige indikatorer, og LCA bruger binære indikatorer. LPA er en probabilistisk model, hvilket betyder, at den modellerer sandsynligheden for, at sagen tilhører en profil. Dette er bedre end en tilgang som K-midler, der bruger afstandsalgoritmer.

med det til side, lad os indlæse dataene.

library(tidyverse)

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

dataene er på 32 interesser/hobbyer. Hvert element er rangeret 1 (ikke interesseret) til 5 (meget interesseret).

beskrivelsen på Kaggle antyder, at der kan være skødesløs reaktion (f.eks. deltagere, der valgte den samme værdi igen og igen). Vi kan bruge den skødesløse pakke til at identificere”string svarende”. Lad os også kigge efter multivariate outliers med Mahalanobis Afstand (se mit tidligere indlæg på Mahalanobis for at identificere outliers).

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

Vi vil cap streng reagerer på et maksimum på 10 og bruge en Mahalanobis D cutoff af alpha = .001.

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

pakken mclust udfører forskellige typer modelbaseret klyngedannelse og dimensionreduktion. Plus, det er virkelig intuitivt at bruge. Det kræver komplette data (ingen mangler), så i dette eksempel fjerner vi sager med NAs. Dette er ikke den foretrukne tilgang; vi ville være bedre stillet tilregne. Men til illustrative formål fungerer dette fint. Jeg vil også standardisere alle indikatorerne, så når vi plotter profilerne, er det tydeligere at se forskellene mellem klynger. Det tager et par minutter at køre denne kode.

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

Vi starter med at plotte Bayesian informationskriterier for alle modeller med profiler fra 1 til 9.

plot(BIC)

det er ikke umiddelbart klart, hvilken model der er bedst, da y-aksen er så stor, og mange af modellerne scorer tæt sammen. sammendrag (BIC) viser de tre bedste modeller baseret på 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

den højeste BIC kommer fra VVE, 3. Dette siger, at der er 3 klynger med variabelt volumen, variabel form, lige orientering og ellipsodial fordeling (se figur 2 fra dette papir for en visuel). Imidlertid, VEE, 3 er ikke langt bagefter og kan faktisk være en mere teoretisk nyttig model, da den begrænser formen på fordelingen til at være ens. Af denne grund går vi med VEE, 3.

hvis vi vil se nærmere på denne model, gemmer vi den som et objekt og inspicerer det med resume().

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

udgangen beskriver profilernes geometriske egenskaber og antallet af sager klassificeret i hver af de tre klynger.

BIC er et af de bedste pasningsindekser, men det anbefales altid at kigge efter flere beviser for, at den løsning, vi har valgt, er den rigtige. Vi kan også sammenligne værdier af det integrerede afsluttede Likelikood (ICL) kriterium. Se dette papir for flere detaljer. ICL er ikke meget forskellig fra BIC, bortset fra at det tilføjer en straf på løsninger med større entropi eller klassificeringsusikkerhed.

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

Vi ser lignende resultater. ICL antyder, at model VEE, 3 passer ganske godt. Endelig udfører vi Bootstrap Sandsynlighed Ratio Test (BLRT), som sammenligner modelpasning mellem k-1 og k cluster modeller. Med andre ord, det ser ud til at se, om en stigning i profiler øger pasformen. Baseret på simuleringer af Nylund, Asparouhovog Muth Larn (2007) BIC og BLRT er de bedste indikatorer for, hvor mange profiler der er. Denne kodelinje vil tage lang tid at køre, så hvis du bare følger med, foreslår jeg at springe den over, medmindre du vil træde ud for en kaffepause.

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 foreslår også, at en 3-profil løsning er ideel.

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