【R言語】政府統計を用いた階層的クラスタリングの実践例
クラスター分析の復習の記事では、食品ごとのたんぱく質摂取量のパターンにもとづいて欧州25ヵ国をグルーピングしており、地理的に近い地域の国でクラスタを形成しているならば、日々の食事習慣みも似通ってくるのではないかと推測していました。なので、ここではクラスタ分析の使い方を覚えることを主眼に置きつつ、都道府県の主要な食品の消費データを利用して、地理的に近い県で似通った食品の消費傾向を示すクラスタが形成されるかを検証してみます。
目次
- 使用するデータ
- 使用するパッケージ
- データの読込み
- 階層的クラスタ分析
- エルボー曲線
- Calinski-Harabasz (CH) インデックス
- 主成分分析による変数の次元削減
- クラスタの安定度
- おわりに
使用するデータ
使用するパッケージ
library(readr) library(dplyr) library(magrittr) library(reshape2) library(ggplot2) library(fpc)
データの読込み
> data <- read_csv("C:/Users/admin/data/shohi_jittai.csv", col_names = TRUE,locale=locale(encoding="SJIS")) > glimpse(data) Observations: 47 Variables: 6 $ Pref <chr> "北海道", "青森県", "岩手県", "宮城県", "秋田県", "山形県", "福島県", "茨城県", "栃木県... $ Cereal <dbl> 6407, 5505, 6008, 6379, 6378, 6567, 6515, 6387, 6325, 6826, ... $ Fish <dbl> 6526, 6666, 6657, 6484, 6948, 6742, 6063, 6026, 6135, 5974, ... $ Meat <dbl> 5990, 5576, 5470, 5799, 5896, 7237, 5664, 5777, 5679, 5362, ... $ Vegi <dbl> 7816, 7857, 8670, 8831, 8491, 10121, 8406, 8193, 8511, 8416,... $ Fruits <dbl> 2703, 2492, 3020, 3054, 2845, 2901, 2761, 2908, 3149, 3039, ...
Prefは都道府県、Cerealは穀類、Fishは魚介類、Meatは肉類、Vegiは野菜・海藻、Fruitsは果物の支出額（円）を表している。
階層的クラスタ分析
はじめにデンドログラムを作成し、視覚的に地理的に近い県で似通った食品の消費傾向を示すクラスタが形成されるかを確認してみたいと思います。
> pmatrix <- scale(data[,-1]) # 変数の正規化 > d <- dist(pmatrix, method="euclidean") #ユークリッドの距離を使用 > pfit <- hclust(d, method="ward.D2") #ward法による階層的クラスタリング # デンドログラムの表示 > plot(pfit, labels=data$Pref, hang=-1)
上図をみると、細かくグルーピングすれば地理的に近い県でクラスタが形成されているようにも見えそうですが、それだとグループが細かくなりすぎてクラスタ分析する意味がなくなると思います。かといって大きく4つ、5つ程度のクラスタに分けてみてもこの図からは何ともいえなさそうです。
なので、参考までに、エルボー曲線やCalinski-Harabasz (CH) インデックスなどを利用してクラスタ数がいくつくらいがベターなのかを見てみたいと思います。
エルボー曲線
sqr_edist <- function(x, y) { sum((x-y)^2) } wss.cluster <- function(clustermat) { c0 <- apply(clustermat, 2, FUN=mean) sum(apply(clustermat, 1, FUN=function(row){sqr_edist(row,c0)})) } wss.total <- function(dmatrix, labels) { wsstot <- 0 k <- length(unique(labels)) for(i in 1:k){ wsstot <- wsstot + wss.cluster(subset(dmatrix, labels==i)) } wsstot } elbow <- function(dmatrix, kmax) { npts <- dim(dmatrix)[1] wss <- numeric(kmax) wss[1] <- (npts-1)*sum(apply(dmatrix, 2, var)) for(k in 2:kmax) { d <- dist(dmatrix, method="euclidean") pfit <- hclust(d, method="ward.D2") labels <- cutree(pfit, k=k) wss[k] <- wss.total(dmatrix, labels) } list(wss = wss) } elbow <- elbow(pmatrix,10) elbowframe <- data.frame(k=1:10,wss=elbow$wss) elbowframe <- melt(elbowframe, id.vars=c("k"), variable.name="measure", value.name="score") ggplot(elbowframe, aes(x=k, y=score, color=measure)) + geom_point(aes(shape=measure)) + geom_line(aes(linetype=measure)) + theme_dark() + scale_x_continuous(breaks=1:10, labels=1:10) + labs(x = "クラスタ数", y = "クラスタ内距離二乗和")
上図をみるとクラスタ数が4くらいがエルボーなのかと思ってしまいますが、判断してもいいか自分には分かりません（もっとハッキリと図に表れてくれたらよかったのに…）
なので、次にCalinski-Harabasz (CH) インデックスを図示してみます。
Calinski-Harabasz (CH) インデックス
# total sum of squares totss <- function(dmatrix) { grandmean <- apply(dmatrix, 2, FUN=mean) sum(apply(dmatrix, 1, FUN=function(row){sqr_edist(row, grandmean)})) } # wss Pseudo_F <- function(dmatrix, kmax) { npts <- dim(dmatrix)[1] # number of rows. totss <- totss(dmatrix) wss <- numeric(kmax) crit <- numeric(kmax) wss[1] <- (npts-1)*sum(apply(dmatrix, 2, var)) for(k in 2:kmax) { d <- dist(dmatrix, method="euclidean") pfit <- hclust(d, method="ward.D2") labels <- cutree(pfit, k=k) wss[k] <- wss.total(dmatrix, labels) } bss <- totss - wss # bssの計算 crit.num <- bss/(0:(kmax-1)) # bss/k-1 crit.denom <- wss/(npts - 1:kmax) # wss/n-k list(crit = crit.num/crit.denom, wss = wss, totss = totss, test= crit.num)# crit = (bss/k-1) / (wss/n-k) } # 視覚化 Fcrit <- Pseudo_F(pmatrix, 10) Fframe <- data.frame(k=1:10, ch=Fcrit$crit) Fframe <- melt(Fframe, id.vars=c("k"), variable.name="measure", value.name="score") ggplot(Fframe, aes(x=k, y=score, color=measure)) + geom_point(aes(shape=measure)) + geom_line(aes(linetype=measure)) + theme_dark() + scale_x_continuous(breaks=1:10, labels=1:10) + labs(x = "クラスタ数", y = "Pseudo F")
CH基準ではクラスタ数が4のときが好ましいという結果になっています。
なので、とりあえず、クラスタ数が4の場合でデンドログラムを再度確認してみます。
> plot(pfit, labels=data$Pref, hang=-1) > rect.hclust(pfit, k=4, border="red") #都道府県を4つに区切る
1つのクラスタは地理的に近い地域（沖縄県、鹿児島県、宮崎県）で形成されていますが、その他のクラスタは地理的に近い地域で形成されているとはいえなさそうです（クラスタが4つだとそうなるのも仕方がないとは思えますが…）。
ここで、クラスタ及びその構成県ごとの食品の支出額を見てみます。
> groups <- cutree(pfit, k=4) > print_clusters <- function(labels, k) { > for(i in 1:k) { > print(paste("cluster", i)) > print(data[labels==i,c("Pref","Cereal","Fish","Meat","Vegi")]) > } > } > print_clusters(groups, 4) [1] "cluster 1" # A tibble: 13 × 5 Pref Cereal Fish Meat Vegi <chr> <dbl> <dbl> <dbl> <dbl> 1 北海道 6407 6526 5990 7816 2 青森県 5505 6666 5576 7857 3 岩手県 6008 6657 5470 8670 4 宮城県 6379 6484 5799 8831 5 秋田県 6378 6948 5896 8491 6 福島県 6515 6063 5664 8406 7 茨城県 6387 6026 5777 8193 8 栃木県 6325 6135 5679 8511 9 群馬県 6826 5974 5362 8416 10 山梨県 6319 6106 6354 8303 11 長野県 6486 6539 5911 8691 12 鳥取県 6228 6296 6494 7886 13 高知県 5900 6336 6095 7702 [1] "cluster 2" # A tibble: 17 × 5 Pref Cereal Fish Meat Vegi <chr> <dbl> <dbl> <dbl> <dbl> 1 山形県 6567 6742 7237 10121 2 埼玉県 6477 6017 6550 9210 3 千葉県 6373 6402 6547 9223 4 東京都 6163 6510 7091 9744 5 神奈川県 6475 6668 7330 9828 6 新潟県 8308 6600 5817 9869 7 富山県 8054 7155 6843 9548 8 石川県 7092 7099 7107 8790 9 福井県 7981 6630 7722 9235 10 静岡県 6904 6393 6544 8930 11 三重県 6996 6781 7616 8408 12 滋賀県 7925 6343 8096 8546 13 京都府 7056 6529 8499 9316 14 大阪府 6877 6269 7959 8926 15 兵庫県 7210 6227 8246 8840 16 奈良県 7041 6796 8226 8806 17 島根県 7391 6410 6908 9089 [1] "cluster 3" # A tibble: 14 × 5 Pref Cereal Fish Meat Vegi <chr> <dbl> <dbl> <dbl> <dbl> 1 岐阜県 6482 5694 6613 8227 2 愛知県 6811 5673 7155 8535 3 和歌山県 7143 6350 8216 7567 4 岡山県 6666 6234 7109 7627 5 広島県 6994 5916 7434 7780 6 山口県 6823 6434 7342 7963 7 徳島県 7144 5901 7552 7793 8 香川県 7041 5626 6907 7684 9 愛媛県 6951 5752 7193 7542 10 福岡県 6199 5549 7348 7970 11 佐賀県 6240 5606 7166 7648 12 長崎県 6203 5734 6931 7472 13 熊本県 5995 5364 7103 7671 14 大分県 6007 5647 7159 7525 [1] "cluster 4" # A tibble: 3 × 5 Pref Cereal Fish Meat Vegi <chr> <dbl> <dbl> <dbl> <dbl> 1 宮崎県 5128 4978 7232 7259 2 鹿児島県 5320 4933 6859 7105 3 沖縄県 5339 3791 5644 7471
ぱっとみて分かるのはクラスタ2の野菜の支出額が高いことくらいでしょうか。
このままだと傾向がわかりづらいので、図(平面)で県をプロットできるように食品の支出額を示す5つの変数を2つの変数へ次元を削減してみます。
主成分分析による変数の次元削減
まず、食品の支出額を示す5つの変数を2つの変数に次元を削減してみます。
princ <- prcomp(pmatrix) summary(princ) Importance of components: PC1 PC2 PC3 PC4 PC5 Standard deviation 1.6768 1.0909 0.64787 0.60864 0.45635 Proportion of Variance 0.5623 0.2380 0.08395 0.07409 0.04165 Cumulative Proportion 0.5623 0.8003 0.88426 0.95835 1.00000
第2主成分までで全体の約8割を説明できていることが分かります。
ここでの目的は、各県を平面上にプロットして視覚的に各クラスタの傾向を把握することなので、今回は縦軸に第2主成分、横軸に第1主成分とする図を以下のように作成する。
project <- predict(princ, newdata=pmatrix)[,1:2] project.plus <- cbind(as.data.frame(project),cluster=as.factor(groups),Prefecture=data$Pref) ggplot(project.plus, aes(x=PC1, y=PC2)) + geom_point(aes(shape=cluster, colour = cluster)) + theme_dark() + geom_text(aes(label=Prefecture),hjust=0, vjust=1)
上図は縦軸に第2主成分(PC2)、横軸に第1主成分(PC1)をとり、クラスタ1から4に属する各都道府県をプロットしています。
クラスタ1に属する各県をみると、PC1はゼロ近辺に、PC2はゼロより大きな値の領域に分布していることが分かります。クラスタ2はPC1がゼロより大きな値に、PC2は負の値から正の値に満遍なく分布していることが分かります。クラスタ3はPC1とPC2がゼロより小さな値の領域に多くの県が分布しており、クラスタ4はPC1とPC2がゼロより小さな値の領域にすべての県が分布していることが分かります。
ここで、PC1とPC2が何を表しているかを上図と上表から考察してみます。まずPC1が大きな値をとるほど野菜への支出額が高いことが分かります。次にPC2が小さな値をとるほど肉への支出額が高いことが分かります。
以上より各クラスタの特徴は次のように表現できるかもしれません。クラスタ1は肉の消費支出が低い県で構成されており、クラスタ2は野菜の消費支出が高い県で構成されていることが分かります。次にクラスタ3と4は肉の消費支出が高く野菜の消費支出が低い傾向にあり、特にクラスタ4は食品の消費支出額の水準が低い（ただし、沖縄県以外の肉の支出額は除く）傾向にあることが読み取れます。
ところで、穀類や魚は各クラスタでハッキリとした傾向が見れませんでした。これは、ごはんや魚は日本の食文化としてどの地方でも習慣的に食されているためにクラスタ間で傾向がでなかったのではと思います。
クラスタの安定度
> kbest.p<-4 > cboot.hclust <- clusterboot(pmatrix,clustermethod=hclustCBI,method="ward.D2", k=kbest.p) > bootMean.data <- data.frame(cluster = 1:kbest.p, bootMeans = cboot.hclust$bootmean) > round(bootMean.data, digit=2) cluster bootMeans 1 0.82 2 0.78 3 0.71 4 0.62
クラスタ4の安定度が低いようですが、bootMeansの値が0.6を下回るような極端に低い値を有するクラスタが存在しなかったため、とりあえず極端に意味のないクラスタは形成されてはいないといえそうです。
おわりに
ここでは都道府県の主要な食品の消費データを利用して、地理的に近い県で似通った食品の消費傾向を示すかをクラスタ分析を利用して検証してみました。結果は、地理的に近い県で似通った食品の消費傾向を示すことを肯定する証拠は得られませんでしたが、47都道府県を4つのグループに分類することがでました。1つ目のグループは野菜の消費が多いグループであり、2つ目が肉の消費が少ないグループ、3つ目が肉の消費が多く、野菜の消費が少ないグループ、4つ目が3つ目と同様の傾向を示しているが、食品への消費水準が他の県より低いグループであることが示されました。
このようなグループがどのような理由で構成されているかは今回の分析からだと何もいえませんが、分析をしてみた感じたことは、47都道府県を4つのグループに分類するのはやはり無理があるのかなと思いました。さらに、各グループの特徴の説明も苦しいことろがあり、その分野に関する専門的な情報が何もない状態ではやはり、分析結果の解釈が難しいなと思いました。
また、地理的に近い地域で似通った食品の消費傾向を示すかを捉えるのに、都道府県のデータを使用するのは本当に適切なのかという疑問や食品の消費に影響を及ぼす他の要因を全く考慮せずに分析をしている点などが気になりました。とはいえ、これはクラスタ分析の練習であり、分析結果をだすことが目的ではなく、クラスタ分析の使い方を覚えることが目的なので、そいうった意味では目的を達成できたのではと思います。
