Rによる階層的クラスタ分析の実践例

クラスター分析の復習の記事では、食品ごとのたんぱく質摂取量のパターンにもとづいて欧州25ヵ国をグルーピングしており、地理的に近い地域の国でクラスタを形成しているならば、日々の食事習慣みも似通ってくるのではないかと推測していました。なので、ここではクラスタ分析の使い方を覚えることを主眼に置きつつ、都道府県の主要な食品の消費データを利用して、地理的に近い県で似通った食品の消費傾向を示すクラスタが形成されるかを検証してみます。

目次

使用するデータ

  • 平成26年全国消費実態調査 都道府県別 家計収支に関する結果
  • 表題:[二人以上の世帯]地域編第42表 地域別1世帯当たり1か月間の収入と支出
  • 利用するデータは穀類、魚介類、肉類、野菜・海藻、果物の支出額(円)

使用するパッケージ

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)

f:id:joure:20170602235853p:plain

上図をみると、細かくグルーピングすれば地理的に近い県でクラスタが形成されているようにも見えそうですが、それだとグループが細かくなりすぎてクラスタ分析する意味がなくなると思います。かといって大きく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 = "クラスタ内距離二乗和")

f:id:joure:20170603235450p:plain

上図をみるとクラスタ数が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") 

f:id:joure:20170603235728p:plain

CH基準ではクラスタ数が4のときが好ましいという結果になっています。
なので、とりあえず、クラスタ数が4の場合でデンドログラムを再度確認してみます。

> plot(pfit, labels=data$Pref, hang=-1)
> rect.hclust(pfit, k=4, border="red") #都道府県を4つに区切る

f:id:joure:20170604000500p:plain

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)

f:id:joure:20170604003123p:plain

上図は縦軸に第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つのグループに分類するのはやはり無理があるのかなと思いました。さらに、各グループの特徴の説明も苦しいことろがあり、その分野に関する専門的な情報が何もない状態ではやはり、分析結果の解釈が難しいなと思いました。

また、地理的に近い地域で似通った食品の消費傾向を示すかを捉えるのに、都道府県のデータを使用するのは本当に適切なのかという疑問や食品の消費に影響を及ぼす他の要因を全く考慮せずに分析をしている点などが気になりました。とはいえ、これはクラスタ分析の練習であり、分析結果をだすことが目的ではなく、クラスタ分析の使い方を覚えることが目的なので、そいうった意味では目的を達成できたのではと思います。