Search

潜在プロファイル分析（LPA）での潜在プロファイル分析の簡単な例は、一連の連続変数（すなわち、指標）に対する応答に基づいて、個人のクラ LPAは、指標項目の応答パターンを生成する観測されていない潜在プロファイルがあることを前提としています。ここでは、興味/趣味に基づいて人々のグループを識別するために、LPAの簡単な例を見ていきます。

ここでは、LPAの簡単な例を見ていきます。

ここでは、 データは、自由に利用可能な若者の調査から来ていますKaggle.com…..

ここでは、私たちが何のために行っているかをスニークピークです：

用語注：人々は用語clusters、profiles、classes、groupsを同じ意味で使用しますが、微妙な違い 私は主にlpaの用語に合わせて、ケースのグループ化を参照するためにプロファイルに固執します。 LPAは、潜在クラス分析（LCA）を含むガウス有限混合モデリングの枝であることに注意する必要があります。 LPAとLCAの違いは概念的であり、計算的ではありません：LPAは連続指標を使用し、LCAはバイナリ指標を使用します。 LPAは確率モデルであり、プロファイルに属するケースの確率をモデル化することを意味します。 これは、距離アルゴリズムを使用するK-meansのようなアプローチよりも優れています。それはさておき、データをロードしてみましょう。

library(tidyverse)

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

データは32の興味/趣味にあります。 各項目は1（興味がない）から5（非常に興味がある）にランク付けされています。

Kaggleの説明は、不注意な応答がある可能性があることを示唆しています（例えば、同じ値を何度も選択した参加者）。 Carelessパッケージを使用して、「文字列応答」を識別できます。 また、Mahalanobis距離を持つ多変量の外れ値を探してみましょう（外れ値を識別するためのMahalanobisに関する私の以前の投稿を参照してください）。

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

最大10に応答する文字列をキャップし、Alpha=のMahalanobis Dカットオフを使用します。001.

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

パッケージmclustは、さまざまなタイプのモデルベースのクラスタリングと次元削減を実行します。 さらに、使用するのは本当に直感的です。 完全なデータ（欠落はありません）が必要なので、この例ではNAsのケースを削除します。 これは好ましいアプローチではありません。 しかし、例示的な目的のために、これは正常に動作します。 私はまた、すべての指標を標準化するつもりですので、プロファイルをプロットすると、クラスター間の違いを見ることがより明確になります。 このコードを実行するには数分かかります。

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

まず、プロファイルが1から9のすべてのモデルのベイズ情報基準をプロットします。p>

plot(BIC)

y軸が非常に大きく、モデルの多くが近くにスコアがあるため、どのモデルが最適かはすぐにはわかりません。 要約（BIC）は、BICに基づいて上位3つのモデルを示しています。最も高いBICはVVE、3から来ています。 これは、可変体積、可変形状、等しい向き、および楕円体分布を持つ3つのクラスターがあることを示しています（視覚的には、この論文の図2を参照）。 しかし、VEE,3はそれほど遅れておらず、実際には分布の形状が等しくなるように制約されるため、理論的にはより有用なモデルかもしれません。 このため、私たちはVEE、3で行くでしょう。このモデルをもっと詳しく見たい場合は、それをオブジェクトとして保存し、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

出力は、プロファイルの幾何学的特性と三つのクラスタのそれぞれに分類されたケースの数を記述します。BICは最良の適合指数の1つですが、選択した解が正しいものであるというより多くの証拠を探すことを常にお勧めします。

また，統合完了Ｌｉｋｅｌｉｋｏｏｄ（ＩＣＬ）基準の値を比較することもできる。 詳細については、この論文を参照してください。 ICLは、より大きなエントロピーまたは分類の不確実性を持つ解にペナルティを追加することを除いて、BICとあまり違いはありません。p>

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

同様の結果が表示されます。 ICLは、モデルVEE、3が非常によく合うことを示唆しています。 最後に、K-1モデルとkクラスターモデルの間のモデル適合を比較するブートストラップ尤度比検定(BLRT)を実行します。 言い換えれば、プロファイルの増加が適合するかどうかを確認するように見えます。 Nylund、Asparouhov、Muthén(2007)によるシミュレーションに基づいて、bicとBLRTは、いくつのプロファイルがあるかを示す最良の指標です。 このコード行は実行に長い時間がかかるため、コーヒーブレイクのためにステップアウトしたい場合を除き、それをスキップすることをお勧めします。

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はまた、3プロファイルのソリューションが理想的であることを示唆しています。

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

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