Rでコロナ陽性者数のクラスタリングをする。
はじめに
すっかりコロナも落ち着いた感じではあるが、厚生労働省はオープンデータとしてコロナの陽性者数を公開してくれている。
このデータを元に「コロナの陽性者数の変化に地理的相関があるのか？」を確かめてみる。つまり、「隣の県で陽性者数が増えたらその県でも増えるのではないか？」という話。目標としては以下のように日本地図を塗りつぶしてみたい。
なんとなく地理的相関はありそうなのが分かる。
この作図までをRで行おうと思う。クラスタリング自体はkmeans関数を呼び出すだけだが、tidyverseでの前処理が少し分かりにくいので、そこを中心に解説する。
参考文献
動作するコードは以下に置いてある。
日本地図を塗りつぶすのは以下のページを参照した。
データの前処理と整形
元データ
厚労省で公開されているデータは以下の様な感じで、県名と2020/1/16からの陽性者数が記載されたCSV形式のファイル。read_csvでDataframeとして読み込んで、後はtidyverseで整形する。
# A tibble: 1,173 × 49
Date ALL Hokkaido Aomori Iwate Miyagi Akita Yamagata Fukushima Ibaraki
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 2020/1/16 1 0 0 0 0 0 0 0 0
2 2020/1/17 0 0 0 0 0 0 0 0 0
3 2020/1/18 0 0 0 0 0 0 0 0 0
4 2020/1/19 0 0 0 0 0 0 0 0 0
日付のParse
まず、CSV(wide)を縦長(long)に変更。次に、日付が文字列になっているので、日付型に変換する。lubridateライブラリを使うと簡単にparseしてくれる。年/月/日の順番なので、ymdを使う。
pivot_longer(!Date, names_to = "prefecture", values_to = "number") %>%
mutate(date = ymd(Date)) %>%
余分なデータの削除
ALL（すべての県の陽性者数の合計）や、文字列形式の日付は不要なので削除する。
filter(prefecture != 'ALL') %>%
select(-Date) %>%
週ごとの平均値を求める
日次データだと変動が大きすぎるので、週毎に平均値を求める。floor_date(date, "week")でその週の最初の日付（日曜日）を取得できるので、dateを全部週初日に変更して、県名→週初日の順にグループ化する。その後summariseで(県名,週初日)のグループ毎に平均を求める。つまり、1週間の陽性者数の平均を県ごとに求める。
mutate(date = floor_date(date, "week")) %>%
group_by(prefecture, date) %>%
summarise(number = mean(number)) %>%
県ごとに陽性者数を正規化する
県によって人口は大きく異なり、陽性者数も大きく異なる。このままクラスタリングすると絶対値での影響が大きくなるので正規化する。つまり、期間中の陽性者数の合計が1になる様にする。summariseの引数に.groups = "drop_last"をつけることで最後のグループ化変数（週初日）のみを解除し、県名でのグループを残す。グループ毎の合計sum(number)で割ることで、正規化できる。最後にungroupですべてのグループ化を解除する。
mutate(date = floor_date(date, "week")) %>%
group_by(prefecture, date) %>%
- summarise(number = mean(number)) %>%
+ summarise(number = mean(number), .groups = "drop_last") %>%
+ mutate(normalized = number / sum(number)) %>%
+ ungroup() %>%
データの整形
最後にkmeans関数に渡しやすいように、整形する。
pivot_wider(names_from = "date", values_from = "normalized", values_fill = 0) %>%
column_to_rownames("prefecture")
クラスタリング
kmeans法
kmeans法を用いたクラスタリングはkmeans関数を使えばすぐにできる。結果をDataFrameにしてmapに対してinner_joinすれば、地図形式にクラスタ番号が振られうことになる。ただ、クラスタ数をあらかじめ決定しておく必要がある。
CLUSTER_NUM <- 7
km <- kmeans(data, CLUSTER_NUM)
cls <- data.frame(name=names(km$cluster), number=km$cluster)
map <- map_nipponmap %>% inner_join(cls, by="name")
クラスタ数の求め方
最適なクラスタ数を求めるにはエルボー法やシルエット法、ギャップ統計法を用いる必要がある。とりあえず全部をグラフ化してみる。出力されたグラフを見ると2か3かなぁという感じ。
library(cluster)
library(factoextra)
library(patchwork)
f1 <- fviz_nbclust(data, kmeans, method = "wss")
f2 <- fviz_nbclust(data, kmeans, method = "silhouette")
f3 <- fviz_nbclust(data, kmeans, method = "gap_stat")
f <- f1 + f2 + f3 + plot_layout(ncol = 1)
ggsave(file="cluster_num.pdf", plot=f , width=8, height=6)
クラスタリング結果
クラスタ数3で導出した結果。kmeansは乱数が入るので、毎回結果が違う。このため、参考程度にしかならないのは困ったものかも知れない。大きく、首都圏, 近畿, 愛知, 福岡でクラスタを形成し、あとは東日本と西日本に分かれている。県間の人の流れとか他の情報とのマッピングをとってみても面白いかもしれない。
Discussion