Search

Analiza profilu utajonego (LPA) próbuje zidentyfikować klastry osób (tj. profile utajone) na podstawie odpowiedzi na szereg zmiennych ciągłych (tj. wskaźniki). LPA zakłada, że istnieją nieobserwowane profile utajone, które generują wzorce odpowiedzi na elementach wskaźnika.

tutaj przejdę przez szybki przykład LPA, aby zidentyfikować grupy ludzi na podstawie ich zainteresowań / hobby. Dane pochodzą z badania młodych ludzi, dostępnego bezpłatnie na Kaggle.com.

oto rzut oka na to, co zamierzamy:

Uwaga dotycząca terminologii: ludzie używają terminów klastry, profile, klasy i grupy zamiennie, ale istnieją subtelne różnice. Będę głównie trzymać się profilu, aby odnieść się do grupowania przypadków, zgodnie z terminologią LPA. Należy zauważyć, że LPA jest gałęzią Gaussa modelowania skończonej mieszaniny, która obejmuje analizę klas utajonych (LCA). Różnica między LPA a LCA jest pojęciowa, a nie obliczeniowa: LPA wykorzystuje wskaźniki ciągłe, a LCA wykorzystuje wskaźniki binarne. LPA jest modelem probabilistycznym, co oznacza, że modeluje prawdopodobieństwo przypadku należącego do profilu. Jest to lepsze od podejścia takiego jak k-means, które wykorzystuje algorytmy odległości.

pomijając to, załadujmy dane.

library(tidyverse)

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

dane dotyczą 32 zainteresowań/hobby. Każdy przedmiot jest w rankingu od 1 (nie zainteresowany) do 5 (bardzo zainteresowany).

opis na Kaggle sugeruje, że mogą być nieostrożne odpowiedzi (np. uczestnicy, którzy wielokrotnie wybierali tę samą wartość). Możemy użyć nieostrożnego pakietu do identyfikacji „string responding”. Poszukajmy również wielowymiarowych wartości odstających z odległością Mahalanobisa (patrz mój poprzedni post na Mahalanobis w celu identyfikacji wartości odstających).

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

ograniczymy ciąg odpowiadający maksymalnie 10 i użyjemy Mahalanobis d odcięcia alpha = .001.

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

pakiet wykonuje różne rodzaje klastrowania opartego na modelu i redukcji wymiarów. Dodatkowo jest bardzo intuicyjny w obsłudze. Wymaga kompletnych danych (nie brakuje), więc w tym przykładzie usuniemy przypadki za pomocą NAs. Nie jest to preferowane podejście; lepiej byłoby przypisać. Ale dla celów ilustracyjnych, to działa dobrze. Ustandaryzuję też wszystkie wskaźniki, więc kiedy wykreślimy profile, wyraźniej będzie zobaczyć różnice między klastrami. Uruchomienie tego kodu zajmie kilka minut.

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

zaczniemy od wykreślenia bayesowskich kryteriów informacyjnych dla wszystkich modeli z profilami od 1 do 9.

plot(BIC)

nie jest od razu jasne, który model jest najlepszy, ponieważ oś y jest tak duża, a wiele modeli osiąga wyniki blisko siebie. podsumowanie (BIC) pokazuje trzy najlepsze modele oparte na 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

najwyższy BIC pochodzi z VVE, 3. To mówi, że istnieją 3 klastry o zmiennej objętości, zmiennym kształcie, równej orientacji i rozkładzie elipsoidalnym (zobacz rysunek 2 z tego papieru dla wizualizacji). Jednak VEE, 3 nie jest daleko w tyle i faktycznie może być bardziej teoretycznie użytecznym modelem, ponieważ ogranicza kształt rozkładu do równego. Z tego powodu wybierzemy VEE, 3.

jeśli chcemy przyjrzeć się temu modelowi dokładniej, zapisujemy go jako obiekt i sprawdzamy za pomocą metody 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

wynik opisuje charakterystykę geometryczną profili i liczbę przypadków sklasyfikowanych w każdym z trzech klastrów.

BIC jest jednym z najlepiej dopasowanych indeksów, ale zawsze zaleca się szukać więcej dowodów na to, że wybrane przez nas rozwiązanie jest prawidłowe. Możemy również porównać wartości kryterium Integrated Completed Likelikood (ICL). Zobacz ten artykuł, aby uzyskać więcej informacji. ICL niewiele różni się od BIC, poza tym, że dodaje karę dla rozwiązań o większej entropii lub niepewności klasyfikacji.

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

widzimy podobne wyniki. ICL sugeruje, że model VEE, 3 pasuje całkiem dobrze. Na koniec przeprowadzimy Test współczynnika prawdopodobieństwa Bootstrap (BLRT), który porównuje dopasowanie modelu między modelami K-1 i K cluster. Innymi słowy, wygląda na to, aby sprawdzić, czy wzrost profili zwiększa dopasowanie. Na podstawie symulacji Nylund, Asparouhov, and Muthén (2007) BIC i BLRT są najlepszymi wskaźnikami liczby profili. Ta linia kodu zajmie dużo czasu, więc jeśli po prostu podążasz, sugeruję pominięcie jej, chyba że chcesz wyjść na przerwę na kawę.

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 sugeruje również, że rozwiązanie 3-profilowe jest idealne.

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