はじめに
今回の記事では、数学的な説明のほとんどを省き、R(Rstudio)を用いた質的（カテゴリ）変数（以下、「カテゴリ変数」とする）を含むクラスター分析の手法に関して記述していきます。
質的変数の扱い（ダミー変数化）
クラスタリングは数値データに関して実施できるので、カテゴリ変数に関してはダミー変数化し数値に変換します。
| 顧客ID | 年齢 | 性別 | 注文経路 | 年間購入金額 |
|---|---|---|---|---|
| 1 | 25 | 女性 | アプリのみ | 600,000 |
| 2 | 62 | 女性 | アプリのみ | 800,000 |
| 3 | 46 | 女性 | 店舗・アプリ併用 | 1,500,000 |
| 4 | 50 | 男性 | 店舗・アプリ併用 | 800,000 |
| 5 | 34 | 男性 | アプリのみ | 450,000 |
上記のようなデータがあった場合に、まずはカテゴリ変数をfactor型に変更します（変更するまではcharacter型になっているはずです）
# データの構造を確認
str(df.test)
# カテゴリ変数をfactor化
df.test$性別 <- as.factor(df.target$性別)
df.test$注文経路 <- as.factor(df.target$注文経路)
このようにfactor化できたら、実際にダミー変数に変えていきます。
# ダミー変数に変換してくれる優れものパッケージ
install.packages("fastDummies")
library(fastDummies)
# ここでダミー変数化
df.test <- dummy_cols(df.test, remove_first_dummy = F, remove_most_frequent_dummy = F)
# ダミー変数化したら元の列は不要なので削除
df.target <- select(df.test, -性別, -注文経路)
| 顧客ID | 年齢 | 性別_男性 | 性別_女性 | 注文経路_アプリのみ | 注文経路_店舗・アプリ併用 | 年間購入金額 |
|---|---|---|---|---|---|---|
| 1 | 25 | 0 | 1 | 1 | 0 | 600,000 |
| 2 | 62 | 0 | 1 | 1 | 0 | 800,000 |
| 3 | 46 | 0 | 1 | 0 | 1 | 1,500,000 |
| 4 | 50 | 1 | 0 | 0 | 1 | 800,000 |
| 5 | 34 | 1 | 0 | 1 | 0 | 450,000 |
これでダミー変数化が終わったのですが、多重共線性の問題が出てきました。
要は「性別_男性」「性別_女性」の列は完全に相関係数が-1になるような組み合わせなのでそもそも片方の列だけでいいのでは？というようなことです。
相関係数の確認
そこで多重共線性をひとまず視覚的に確認するためにcorrplot関数を使用します。
これを使用することで、以下のように相関係数と有意なものを一気に確認することができます。
# パッケージのインストールとライブラリの読み込み
install.packages("corrplot")
library(corrplot)
c <- cor(test2)
p <- cor.mtest(test2)
# 小数第２位まで表示させるような設定
round(c,2)
corrplot(c, method = "color", addCoef.col = T, p.mat = p$p, sig.level = 0.001)
もし、これで出力した際に（上記の画像のように）日本語が豆腐化していた場合は、par(family= "HiraKakuProN-W3")を実行してから再度出力してください。
ダミー変数化（多重共線性の解消）
多重共線性を解消するために便利なのがcarというパッケージです。VIF統計量（多重共線性がどれくらいあるのかの尺度）を出力してくれます。
# 多重共線性への影響考慮
library(car)
vif_values <- vif(lm(cluster~., data = df.test))
VIF統計量は一般的にに10以下であれば多重共線性がないとされていますが、理想値は2以下です。
VIF統計量が10を超えた変数がある場合にはモデルからその変数を外してもう一度VIF統計量を計算するなど、モデルを再構成する必要が出てきます。
正規化の実施
次に正規化を実施します。年齢は基本的に20から80程度の範囲で収まりますが、年間購入額は0〜1,500,000とかなり幅広くなってしまいます。また年齢は２０スタートにも関わらず金額は0からスタートしています。そこで活躍するのが正規化です。
正規化とは、平均を0、分散を1にスケーリングすることです。
これはscale関数でかなり簡単に実施できます。
# df.testを正規化
df.test <- data.frame(scale(df.test))
これのポイントはdata.frame関数を使ってデータフレームとして格納することを明記することです。
ここまできたらクラスタリングの前処理は終わりなので、今から実際にクラスタリングに入っていきます。今回は大規模データ（40万行程度）を想定しているので、非階層クラスタリングをメインに実施します。
非階層クラスタリング
クラスタ数の決定（エルボー法）
クラスタ数を決定するためにエルボー法をしばしば使用することがあります。
上記画像は横軸がクラスタ数になっており、これがいきなり緩やかになっているポイントが適切なクラスタ数と判断できます。
このエルボー法によるクラスタ数の決定は以下のようにkmeans関数を使用して実施します。
wss <- sapply(1:10, function(k) {
kmeans(test.df.target, k, nstart=10)$tot.withinss
})
# エルボープロットを作成
elbow_plot <- data.frame(
K = 1:10,
WSS = wss
)
# ggplotで書き出し
ggplot(elbow_plot, aes(x=K, y=WSS)) +
geom_line() +
geom_point() +
labs(title="Elbow Method for Optimal k", x="Number of clusters (k)", y="Total within-cluster sum of squares (WSS)")
クラスタ数の決定（エルボー法含め全て）
エルボー法以外にも、シルエット分析やギャップ統計量など複数の指標があるのですが、これらを一気に計算して、多数決的に最適なクラスタ数を提案してくれる神のような関数NbClustを使用します。
こちらは、階層的・非階層的クラスタリングを問わず使用できます。
今回は参考として、凝集型階層的クラスタリングではウォード法を用いて、
非階層的クラスタリングではＫ平均法を用いて計算を実施します。
そして、それぞれから得た結果をfviz_nbclusterで描画します。
# パッケージのインストールとライブラリの読み込み
install.packages("NbClust")
library(NbClust)
# ここで読み込む
myAHCnum = NbClust(test2, method = "ward.D", index = "all")
myNHCnum = NbClust(test2, method = "kmeans", index = "alllong")
fig1 = fviz_nbclust(myAHCnum, method = "silhouette")
fig2 = fviz_nbclust(myNHCnum, method = "gap_stat", nboot = 100)
# 描画する
multiplot(fig1, fig2)
k-means法を使用した非階層クラスタリング
k-means法に関しては巷に情報が溢れているので説明は省略して一気にコードを紹介します。
# 非階層クラスタリング
df.test <- select(df.test, -cluster)
km <- kmeans(df.test, 4, iter.max = 30) # この４の部分がクラスタ数
head(km$cluster)
table(km$cluster)
# クラスタリング結果をデータフレームに追加
df.test$cluster <- as.factor(km$cluster)
# 主成分分析（PCA）を実行
pca <- prcomp(df.test[, -ncol(df.test)], scale. = TRUE)
# PCAの結果とクラスタリング結果を含むデータフレームを作成
pca_df <- data.frame(PC1 = pca$x[, 1], PC2 = pca$x[, 2], Cluster = df.test$cluster)
# 散布図を作成
ggplot(pca_df, aes(x = PC1, y = PC2, color = Cluster)) +
geom_point(size = 2, alpha = 0.7) +
labs(x = "PC1", y = "PC2", title = "非階層クラスタリングの結果") +
theme_minimal() +
theme(legend.position = "bottom")
PCA（主成分分析）は、散布図による可視化を実施するために行っています。
すると以下のような画像が出力されます。
また、PCAで各因子にどのような要素が含まれているのかを可視化する際には以下のようなコードを実施します。
# 主成分の負荷量を取得
loadings <- pca$rotation[, 1:2]
# 変数名を取得
var_names <- rownames(loadings)
# 負荷量の絶対値を計算
abs_loadings <- abs(loadings)
# 負荷量の降順にソート
sorted_loadings <- abs_loadings[order(abs_loadings[, 1], decreasing = TRUE), ]
# 上位10個の変数を選択
top_vars <- head(var_names[order(abs_loadings[, 1], decreasing = TRUE)], 10)
# 負荷量の絶対値を含むデータフレームを作成
loadings_df <- data.frame(Variable = rep(var_names, 2),
PC = rep(c("PC1", "PC2"), each = length(var_names)),
Loading = c(abs_loadings[, 1], abs_loadings[, 2]))
# 上位10個の変数に絞り込む
loadings_df <- loadings_df[loadings_df$Variable %in% top_vars, ]
# 棒グラフを作成
ggplot(loadings_df, aes(x = reorder(Variable, Loading), y = Loading, fill = PC)) +
geom_bar(stat = "identity", position = "dodge") +
coord_flip() +
labs(x = "Variable", y = "Absolute Loading", title = "Variable Contribution to Principal Components") +
theme_minimal() +
theme(legend.position = "bottom") +
theme_gray (base_family = "HiraKakuPro-W3")
おまけ：次元削減方法に関するその他手法
PCA以外にも、UMAP法やt-sne法などがあります。
しかしこれらに関しては、重複行があると使えないといったデメリットもあるので正直あまり使えないかなあという印象です。
（顧客情報がたまたま重複していることなんてデータが多ければ可能性は高いですからね。。。）
<details><summary>UMAP法のサンプルコード</summary>
以下のようなコードで実施できます
# ライブラリの読み込み
library(umap)
# umap法を用いた次元削減
test2 <- test.df.target |> select(-cluster)
umap_result <- umap(test2)
# 可視化
plot(umap_result$layout, col = clusters, pch = 19, main = "UMAP")
</details>
階層クラスタリングに関して
階層クラスタリングは一瞬で終わります。
（ただし計算量は多いので、データ数によっては時間がかかります）
以下のコードを実施します。
# 距離行列の作成
dist_matrix <- dist(df.test)
# 階層的クラスタリングの実行
hc <- hclust(dist_matrix, method = "ward.D2")
# デンドログラムのプロット
plot(hc)
以上です。
さいごに
綺麗にクラスタリングできなかった場合に関しては、勇気を持ってクラスターが綺麗に分割できなかったと判定することが大事です。
あくまで、クラスタリングは１つの方法なので。
Discussion