はじめに

時系列データをk-meansで分類する。本来であれば、時系列データのクラスタリング手法であるk-Shapeまたは他の手法のほうが妥当なんだろうけども…今回は単純にk-meansをやってみて、時系列の傾向を掴めそうかためしてみる(そんなのデータによるだろうけども…)。

k-Shapeについては、“k-Shapeによる時系列クラスタリングの論文：「k-Shape: Efficient and Accurate Clustering of Time Series」を読んだ”がわかりやすかった。

サンプルデータ

サンプルデータは{gapminder}のデータ。

library(tidyverse)
library(gapminder)

head(gapminder, 10)

## # A tibble: 10 × 6
##    country     continent  year lifeExp      pop gdpPercap
##    <fct>       <fct>     <int>   <dbl>    <int>     <dbl>
##  1 Afghanistan Asia       1952    28.8  8425333      779.
##  2 Afghanistan Asia       1957    30.3  9240934      821.
##  3 Afghanistan Asia       1962    32.0 10267083      853.
##  4 Afghanistan Asia       1967    34.0 11537966      836.
##  5 Afghanistan Asia       1972    36.1 13079460      740.
##  6 Afghanistan Asia       1977    38.4 14880372      786.
##  7 Afghanistan Asia       1982    39.9 12881816      978.
##  8 Afghanistan Asia       1987    40.8 13867957      852.
##  9 Afghanistan Asia       1992    41.7 16317921      649.
## 10 Afghanistan Asia       1997    41.8 22227415      635.

データ加工

このままでは計算できないので、データを加工していきます。{tidyr}のpivot_wider()を使います。pivot_wider()は最近(2019-11-14時点)実装されたので、使い方はドキュメントを参照ください。やっていることは単純で必要な変数を選んで、Long型をWide型に変換しています。

gapminder_wide <- gapminder %>% 
  dplyr::select(country, year, lifeExp) %>% 
  tidyr::pivot_wider(names_from = year, values_from = lifeExp)


head(gapminder_wide, 10)

## # A tibble: 10 × 13
##    country `1952` `1957` `1962` `1967` `1972` `1977` `1982` `1987` `1992` `1997`
##    <fct>    <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
##  1 Afghan…   28.8   30.3   32.0   34.0   36.1   38.4   39.9   40.8   41.7   41.8
##  2 Albania   55.2   59.3   64.8   66.2   67.7   68.9   70.4   72     71.6   73.0
##  3 Algeria   43.1   45.7   48.3   51.4   54.5   58.0   61.4   65.8   67.7   69.2
##  4 Angola    30.0   32.0   34     36.0   37.9   39.5   39.9   39.9   40.6   41.0
##  5 Argent…   62.5   64.4   65.1   65.6   67.1   68.5   69.9   70.8   71.9   73.3
##  6 Austra…   69.1   70.3   70.9   71.1   71.9   73.5   74.7   76.3   77.6   78.8
##  7 Austria   66.8   67.5   69.5   70.1   70.6   72.2   73.2   74.9   76.0   77.5
##  8 Bahrain   50.9   53.8   56.9   59.9   63.3   65.6   69.1   70.8   72.6   73.9
##  9 Bangla…   37.5   39.3   41.2   43.5   45.3   46.9   50.0   52.8   56.0   59.4
## 10 Belgium   68     69.2   70.2   70.9   71.4   72.8   73.9   75.4   76.5   77.5
## # ℹ 2 more variables: `2002` <dbl>, `2007` <dbl>

k-meansを行うときは標準化するほうがいいとか言われますが、今回は同じようなスケールのはずなので、一旦、標準化はなしで。

エルボー法でクラスタ数を考える

tot.withinssはクラスター内平方和のことで、各クラスター内の観測値のばらつき具合みたいなもの。クラスター内平方和が小さなクラスターはコンパクトで、平方和が大きなクラスターは、クラスター内のばらつきが大きくなります。

なので、クラスター数が少ないときは大きくなりやすいので、その状態からクラスター数を増やしていくことで、クラスター内平方和が大きく変化するのであれば、クラスターを増やす意味があるけども、クラスター数を増やしても、クラスター内平方和がそこまで変わらないのであれば、解釈の観点からも少ないほうが良いとなる。その変化曲線が肘みたいなのでエルボー法。

エルボー法で計算した結果を可視化すると、今回は3とか4で良さそう。

# compute Within-Cluster-Sum of Squared Errors(=wss)
n_clust <- 1L:10L
wss <- n_clust %>%
  purrr::map_dbl(
    .x = .,
    .f = function(x) {
      kmeans(x = gapminder_wide %>% select(-country),
             centers = x,
             nstart = 50,
             iter.max = 15)$tot.withinss})

tibble(n_clust = n_clust, wss = wss) %>% 
  ggplot(.) +
  geom_line(aes(x = n_clust, y = wss), col = "#006E4F") +
  scale_x_continuous(breaks = n_clust) +
  theme_bw()

K-means

クラスタ数は4つでK-meansをやってみる。その後は、Wide型のデータに対して、クラスタを付与して、可視化するためにLong型に変換しています。

set.seed(1989)
clusters <- kmeans(x = gapminder_wide %>% select(-country),
                   centers = 4)

gapminder_long <- gapminder_wide %>% 
  dplyr::bind_cols(tibble(cluster = clusters$cluster)) %>% 
  tidyr::pivot_longer(cols = c(-country, -cluster), names_to = "year", values_to = "lifeExp")

head(gapminder_long, 10)

## # A tibble: 10 × 4
##    country     cluster year  lifeExp
##    <fct>         <int> <chr>   <dbl>
##  1 Afghanistan       4 1952     28.8
##  2 Afghanistan       4 1957     30.3
##  3 Afghanistan       4 1962     32.0
##  4 Afghanistan       4 1967     34.0
##  5 Afghanistan       4 1972     36.1
##  6 Afghanistan       4 1977     38.4
##  7 Afghanistan       4 1982     39.9
##  8 Afghanistan       4 1987     40.8
##  9 Afghanistan       4 1992     41.7
## 10 Afghanistan       4 1997     41.8

クラスタの中心点もクラスタ毎に可視化できるようにしておく。

centers_long <- tibble::rownames_to_column(as_tibble(clusters$centers), "cluster") %>%
  tidyr::pivot_longer(cols = -cluster, names_to = "year", values_to = "lifeExp") 

head(centers_long, 10)

## # A tibble: 10 × 3
##    cluster year  lifeExp
##    <chr>   <chr>   <dbl>
##  1 1       1952     45.5
##  2 1       1957     48.8
##  3 1       1962     51.4
##  4 1       1967     54.4
##  5 1       1972     57.5
##  6 1       1977     60.3
##  7 1       1982     63.0
##  8 1       1987     65.7
##  9 1       1992     68.0
## 10 1       1997     69.8

可視化

可視化した結果を見る限り、考察としてはこんな感じだろうか。

ggplot() +
  geom_line(data = gapminder_long, aes(y = lifeExp, x = year, group = country), col = "gray") +
  geom_line(data = centers_long, aes(y = lifeExp, x = year, group = cluster), col = "#006E4F", size = 1) +
  facet_wrap( ~ cluster, nrow = 1) + 
  theme_bw()

クラスタ1は1952年~2007年にかけて相対的に急激に成長している国
クラスタ2は1952年~2007年にかけて相対的に緩やかに成長している国
クラスタ3は1952年~1990年にかけて相対的に急激に成長したが、下降した国
クラスタ4は1952年~2007年にかけて相対的に成長速度も遅く、下降していたりする国

うん…今回のデータであれば、クラスタ3と4でもう少し、上手く別れてくれるといいが、成長の仕方について、そこそこうまく分類できてそうなので、今回のデータであればよさそう。なので、とりあえず計算コストがそこまでかからないので、k-meansを時系列にやってみるのはありかも。

Memo〜時系列データをk-meansで分類〜

はじめに

サンプルデータ

データ加工

エルボー法でクラスタ数を考える

K-means

可視化