マハラノビスの距離を応用したクラスタリングの話
つい先日、日本マーケティング・リサーチ協会（JMRA）のリサーチ・イノベーション委員会より「マハラノビス研究会報告」というレポートがひっそりと出ています。
www.jmra-net.or.jp
これはタイトルの通り、マハラノビス距離について、マーケティングリサーチでの応用を念頭に、朝野熙彦先生が中心となってまとめられたものです。
業界団体の活動として非常に興味深い試みだと思うんですがイマイチ人の目に触れてない気がしてもったいない。
特に６.４ クラスター分析と汎距離の節は応用しやすそうなので記事を書いてみます。
マハラノビス距離とは？
この図がわかりやすいます。
左の図の円はcenterからのユークリッド距離、我々が普段使っている直線距離を示しています。
これに対して右の図の楕円はマハラノビス距離を示しています。
右の図ではxとyに相関関係があって、その相関関係を考慮した距離がマハラノビスの距離です。
従来、マハラノビスの距離は異常値検知などに活用されてきており、MT法などの応用が知られています。
クラスター分析でのマハラビノスの距離の応用
さて、「マハラノビス研究会報告」にはマーケティング・リサーチにおける従来のクラスター分析の問題点として以下のように書かれています。
マーケティングリサーチの実務では、消費者をセグメントする際に、あらかじめ原データを因子得点や主成分得点に変換してからクラスター分析することが行われてきた。因子分析の主因子解および主成分は直交しているので、ユークリッド距離で消費者間の距離を測ってよいはずだ、というのが従来の通説だったようである。
しかし空間が複数の集団から成り、しかも集団によって分散共分散行列が異なるという一般的な事態に、はたしてこの理屈は通用するのだろうか。なぜなら、データを一括して直交化したとしても、個々の集団内の分布まで無相関になるわけではない。したがって直交化後のユークリッド距離と集団ごとに測った汎距離は一致しない。集団が複数の場合の距離の測定はこれまで見過ごされてきたのではないだろうか。本節ではクラスターごとに汎距離を測る方法を検討する。
うむ、よく分からん。
これ、どういうことを言っているかというと、下図のように異なった相関関係を持ついくつかのグループが
混在する場合、これまでよく使われてきた直交回転の因子分析からのクラスター分析では歯が立たないよね、という話です。
これはちょっと前に界隈で話題になった「k-meansのクラスタ数決めるのにエルボー法使うな！」という論文の中でも指摘されています。
これ、非常に面白い論文なので読むとよいです。
そもそも今どきエルボー法なんて使ってる人がいるのか？？という疑問もありますが……
arxiv.org
で、こういったケースではマハラノビスの距離を応用すると上手いことクラスタリングできるんじゃない？というのが「マハラノビス研究会報告」で言われていること、と私は理解しました。
ということで、さっそくRで試してみましょう。
準備
まず、必要な各種パッケージの呼び出しと、後述するマハラノビスの距離を応用したk-means法の実装を読み込みます。
library(conflicted) #関数名が衝突する場合にエラーを返す library(tidyverse) #データ処理全般 library(psych) # 因子分析 library(cluster) #pam library(mclust) # モデルベースのクラスタリングのための混合ガウスモデル(GMM) library(withr) #乱数種固定 library(patchwork) # 複数のグラフを組み合わせる library(mvtnorm) # mkmeans関数の内部で使用 library(beepr) #音を鳴らす # 改良k-meansの実装。以下の記事から。 # https://qiita.com/sbtseiji/items/3ec0f26bad6a72a698a7 source("mkmeans.r")
次に3つのデモデータを生成します。
シンプソンのパラドクス、分布が交差するケース、バラツキが異なるケースの3パターンです。
n <- 100 demo <- MASS::mvrnorm( n = n, mu = rep(0, 2), Sigma = matrix(c(3, 2, 3, 2), 2, 2) ) |> with_seed(123, code = _) |> as.data.frame() # デモデータその１シンプソンのパラドクス demo1 <- demo |> bind_rows( demo |> mutate(V1 = V1 + 2, V2 = V2 - 2), demo |> mutate(V1 = V1 - 2, V2 = V2 + 2), ) |> mutate( V1 = V1 |> scale(), V2 = V2 |> scale(), cluster = c(rep(1, n), rep(2, n), rep(3, n)) ) # デモデータその２分布が交差するケース demo2 <- demo |> bind_rows( demo |> mutate(V1 = V1 * -1.5 + 2, V2 = V2 * 1 + 5), demo |> mutate(V1 = V1, V2 = V2 + 4) ) |> mutate( V1 = V1 |> scale(), V2 = V2 |> scale(), cluster = c(rep(1, n), rep(2, n), rep(3, n)) ) # デモデータその3バラツキが異なるケース demo3 <- rbind( rnorm(n * 4, 0, 1.5) |> with_seed(1234, code = _) |> matrix(ncol = 2) |> as.data.frame(), rnorm(n * 1, 0, 0.5) |> with_seed(1234, code = _) |> matrix(ncol = 2) |> as.data.frame() |> mutate(V1 = V1 + 5, V2 = V2 + 1), rnorm(n * 1, 0, 0.5) |> with_seed(1234, code = _) |> matrix(ncol = 2) |> as.data.frame() |> mutate(V1 = V1 + 5, V2 = V2 - 1) ) |> mutate( V1 = V1 |> scale(), V2 = V2 |> scale(), cluster = c(rep(1, n * 2), rep(2, n / 2), rep(3, n / 2)) )
どんなデータなのか可視化してみましょう。
plot_data <- function(data, cl, title) { data |> ggplot( aes( x = V1, y = V2, color = as.factor(cl), fill = as.factor(cl), shape = as.factor(cl) ) ) + geom_point() + stat_ellipse(type = "norm", geom = "polygon", alpha = 0.3) + theme(legend.position = "none") + labs(title = title) } p1 <- plot_data(demo1, demo1$cluster, "demo1, シンプソンのパラドクス") p2 <- plot_data(demo2, demo2$cluster, "demo2, 分布が交差するケース") p3 <- plot_data(demo3, demo3$cluster, "demo3, バラツキが異なるケース") p1 + p2 + p3 + plot_layout(ncol = 2) + plot_annotation(title = "3つのデモデータ")
検証
これらのデモデータに対していくつかのクラスタリング手法を適用して、結果を比べてみます。
試すクラスタリング手法は以下の5つです。
このうち改良k-means法は豊田先生が開発されたマハラノビスの距離を応用したクラスタリング手法で、こちらに論文があります。
これをRで実装されている方がいらっしゃるので今回はこちらを使わせていただきます。
qiita.com
ではやってみましょう。
# 各手法をまとめて実行できるように関数化 clustering <- function(data) { data_subset <- data |> dplyr::select(!cluster) res <- data |> mutate( cl_km = data_subset |> kmeans(3) |> with_seed(1234, code = _) |> pluck("cluster"), cl_fakm = data_subset |> fa(nfactors = 2, rotate = "varimax") |> pluck("scores") |> kmeans(3) |> with_seed(1234, code = _) |> pluck("cluster"), cl_pam = data_subset |> pam(3) |> with_seed(1234, code = _) |> pluck("clustering"), cl_mkm = data_subset |> mkmeans(3, iter.max = 100) |> with_seed(1234, code = _) |> pluck("cluster"), cl_mclust = data_subset |> Mclust() |> with_seed(1234, code = _) |> pluck("classification") ) beep(2) #終わったら音を鳴らす return(res) } # まとめて描画する関数 plot_cluster <- function(res, subt) { p1 <- plot_data(res, res$cluster, "正解") p2 <- plot_data(res, res$cl_km, "k-means法") p3 <- plot_data(res, res$cl_fakm, "因子分析 + k-means法") p4 <- plot_data(res, res$cl_pam, "pam") p5 <- plot_data(res, res$cl_mkm, "改良k-means法") p6 <- plot_data(res, res$cl_mclust, "モデルベースクラスタリング") p1 + p2 + p3 + p4 + p5 + p6 + plot_layout(ncol = 3) + plot_annotation( title = "クラスタリング手法比較", subtitle = subt) } # デモデータを各クラスタリング手法にかける res1 <- clustering(demo1) res2 <- clustering(demo2) res3 <- clustering(demo3) # 結果を可視化 plot_cluster(res1, "Demo1, シンプソンのパラドクス") plot_cluster(res2, "Demo2, 分布が交差するケース") plot_cluster(res3, "Demo3, バラツキが異なるケース")
このデータはいわゆるシンプソンのパラドクスで、各グループが同じ相関を持つ場合です。
因子分析+k-means法、pam、モデルベースクラスタリングはうまくグループ化してくれていますが、k-means法と改良k-means法はダメなようです。
こちらは相関が異なるグループが混在し、なおかつ分布が交差しているケースです。
こちらは改良k-means法とモデルベースクラスタリングがいい感じです。
最後は分散が異なるグループが混在するケース。
モデルベースクラスタリング以外は全然です。
まとめ
こうなると「モデルベースクラスタリング一択じゃん！」と言いたくなりますが、そうは問屋が卸しません。
「マハラノビス研究会報告」では以下のように指摘されています。
クラスター分析には正規混合法(GMM)のようなモデルベースのクラスター分析もあるが、現実のデータにモデルベースの理論が適合する保証はない
ということで振出しに戻りますｗ
今回は便宜上、正解のあるデータをもとにそれぞれの手法を検証しましたが、クラスタリングは判別のための手法ではなく「いい感じにグループ分けしてくれる」手法でしかありません。
特にマーケティングの領域に置いては、外部に判断基準を置いてどのクラスタ分けを採用するのが適切か判断するべきでしょう。
そのためにも複数のクラスタリング手法を使えるようにしておくことは大切です。
ということで、皆さんのお役に立てば幸いです。
