2021年9月7日火曜日

Googleモビリティレポートのデータを時系列クラスタリングする

Googleがコミュニティモビリティレポートを公開している。これはGoogleアカウントのロケーション履歴を有効にしているユーザー（サンプルベース）のデータに基づいた移動データで、2020 年 1 ⽉ 3 ⽇〜2 ⽉ 6 ⽇の 5 週間における該当曜⽇の中央値を基準とした増減のデータになっている。詳細はモビリティ レポートの CSV ドキュメントを参照。今回はこのデータを使って、Pythonで時系列クラスタリングをすると同様にPythonのtslearnを使って時系列クラスタリングを行う。


環境


WSL2(Ubuntu20.04)とJupyterLab。
 

その他必要なPythonライブラリをインストール。


モビリティレポートの可視化

モビリティレポートのcsvには世界中のデータが含まれるが、ここでは日本の都道府県ごとのデータのみを使う。まずはcsvをPandasのDataFrameとして読み込んで中身を確認。sub_region_1が都道府県のカラムで、retail_and_recreation_percent_change_from_baselineからresidential_percent_change_from_baselineが場所カテゴリごとの移動データになっている。

1
2
3
4
5
6
7
8
import pandas as pd
 
# Googleモビリティレポートのファイルパス
# 2020年1月3日から2月6日の5週間の曜日別中央値を基準とした値
MOBILITY_CSV_GOOGLE = 'Global_Mobility_Report.csv'
 
dfg = pd.read_csv(MOBILITY_CSV_GOOGLE)
dfg.info()



7月の都道府県ごとのtransit_stations_percent_change_from_baseline（公共交通機関の拠点）のデータを可視化してみる。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
%matplotlib inline
 
import matplotlib.pyplot as plt
import pandas as pd
 
# Googleモビリティレポートのファイルパス
# 2020年1月3日から2月6日の5週間の曜日別中央値を基準とした値
MOBILITY_CSV_GOOGLE = 'Global_Mobility_Report.csv'
 
# 対象とする場所のカテゴリ
# 公共交通機関の拠点（例: 地下鉄、バス、電車の駅）など
MOBILITY_CATEGORY = 'transit_stations_percent_change_from_baseline'
 
dfg = pd.read_csv(MOBILITY_CSV_GOOGLE)
# 都道府県ごとのデータを取得
dfg_jp = dfg[dfg['country_region_code']=='JP'].dropna(subset=['sub_region_1'])
 
# 都道府県のチェック
pref = dfg_jp['sub_region_1'].unique()
print(pref)
assert len(pref) == 47
 
# 日付をインデックス、都道府県をカラムとするDataFrameに変換
dfg_cat = dfg_jp.set_index(['date', 'sub_region_1'])[MOBILITY_CATEGORY].unstack()
 
# 2021年7月のデータを取得
dfg_cat_jul = dfg_cat['2021-07-01':'2021-07-31']
 
# データのプロット
ax = dfg_cat_jul.plot(figsize=(12,8), alpha=0.75, rot=90);
ax.legend(loc='center left', bbox_to_anchor=(1.0, 0.5))
plt.show()

47都道府県分をすべてプロットしているので見づらいグラフになっているが、4連休のあたりで都道府県ごとの違いが大きいように見える。


クラスター数を決める

データをk-means法で時系列クラスタリングするが、k-means法では事前にクラスター数を決める必要がある。クラスタリングの前にシルエット値からクラスター数を決定する。シルエット値はk-means法の最適なクラスター数を選択するでも使用しているが、今回はSelecting the number of clusters with silhouette analysis on KMeans clusteringを参考にシルエット値をプロットして確認してみる。なお、メトリックとしてユークリッド距離を使う。

以下はクラスター数２～５の範囲でシルエット値を算出してプロットするコード。scikit-learnにはデータごとのシルエット値を返してくれるsilhouette_samplesメソッドがあるがtslearnにはないので自作した。平均シルエット値が１に近く、プロット幅がなるべく均等になるクラスター数がよいと言われている。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
%matplotlib inline
 
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
import pandas as pd
 
from tslearn.clustering import TimeSeriesKMeans, silhouette_score
from tslearn.preprocessing import TimeSeriesScalerMinMax
from tslearn.utils import to_time_series_dataset
from scipy.spatial.distance import cdist
from sklearn.metrics import silhouette_samples as sklearn_silhouette_samples
 
# Googleモビリティレポートのファイルパス
# 2020年1月3日から2月6日の5週間の曜日別中央値を基準とした値
MOBILITY_CSV_GOOGLE = 'Global_Mobility_Report210814.csv'
 
# 対象とする場所のカテゴリ
# 公共交通機関の拠点（例: 地下鉄、バス、電車の駅）など
MOBILITY_CATEGORY = 'transit_stations_percent_change_from_baseline'
 
 
def silhouette_samples(X, labels, **kwds):
    """シルエット値の計算
    ユークリッド距離の計算方法は以下tslernのコードを参照
    https://github.com/tslearn-team/tslearn/blob/42a56cc/tslearn/clustering/utils.py#L66-L197
    """
 
    X_ = X.reshape((X.shape[0], -1))
    sklearn_X = cdist(X_, X_, metric='euclidean')
 
    # 計算済みのデータを渡すのでmetric='precomputed'にする
    return sklearn_silhouette_samples(X=sklearn_X,
                                  labels=labels,
                                  metric='precomputed',
                                  **kwds)
 
 
def plot_silhouette_scores(ts_data, silhouette_values, n_clusters, km):
    """シルエット値のプロット
    以下のコードをほぼコピー
    https://scikit-learn.org/stable/auto_examples/cluster/plot_kmeans_silhouette_analysis.html#sphx-glr-auto-examples-cluster-plot-kmeans-silhouette-analysis-py
    """  
    cluster_labels = km.labels_
    centers = km.cluster_centers_
     
    # Create a subplot with 1 row and 2 columns
    fig, (ax1, ax2) = plt.subplots(1, 2)
    fig.set_size_inches(18, 7)
 
    ax1.set_xlim([-0.1, 1])
    # The (n_clusters+1)*10 is for inserting blank space between silhouette
    # plots of individual clusters, to demarcate them clearly.
    ax1.set_ylim([0, len(ts_data) + (n_clusters + 1) * 10])
 
    y_lower = 10
    for i in range(n_clusters):
        # Aggregate the silhouette scores for samples belonging to
        # cluster i, and sort them
        ith_cluster_silhouette_values = \
            silhouette_values[cluster_labels == i]
 
        ith_cluster_silhouette_values.sort()
 
        size_cluster_i = ith_cluster_silhouette_values.shape[0]
        y_upper = y_lower + size_cluster_i
 
        color = cm.nipy_spectral(float(i) / n_clusters)
        ax1.fill_betweenx(np.arange(y_lower, y_upper),
                          0, ith_cluster_silhouette_values,
                          facecolor=color, edgecolor=color, alpha=0.7)
 
        # Label the silhouette plots with their cluster numbers at the middle
        ax1.text(-0.05, y_lower + 0.5 * size_cluster_i, str(i))
 
        # Compute the new y_lower for next plot
        y_lower = y_upper + 10  # 10 for the 0 samples
 
    ax1.set_title("The silhouette plot for the various clusters.")
    ax1.set_xlabel("The silhouette coefficient values")
    ax1.set_ylabel("Cluster label")
 
    # The vertical line for average silhouette score of all the values
    ax1.axvline(x=np.mean(silhouette_values), color="red", linestyle="--")
 
    ax1.set_yticks([])  # Clear the yaxis labels / ticks
    ax1.set_xticks([-0.1, 0, 0.2, 0.4, 0.6, 0.8, 1])
 
    # 2nd Plot showing the actual clusters formed
    colors = cm.nipy_spectral(cluster_labels.astype(float) / n_clusters)
    ax2.scatter(ts_data[:, 0], ts_data[:, 1], marker='.', s=30, lw=0, alpha=0.7,
                c=colors, edgecolor='k')
 
    # Labeling the clusters
    # centers = clusterer.cluster_centers_
    # Draw white circles at cluster centers
    ax2.scatter(centers[:, 0], centers[:, 1], marker='o',
                c="white", alpha=1, s=200, edgecolor='k')
 
    for i, c in enumerate(centers):
        ax2.scatter(c[0], c[1], marker='$%d$' % i, alpha=1,
                    s=50, edgecolor='k')
 
    ax2.set_title("The visualization of the clustered data.")
    ax2.set_xlabel("Feature space for the 1st feature")
    ax2.set_ylabel("Feature space for the 2nd feature")
 
    plt.suptitle(("Silhouette analysis for KMeans clustering on sample data "
                  "with n_clusters = %d" % n_clusters),
                 fontsize=14, fontweight='bold')
 
    plt.show()
 
 
def show_silhouette_scores(ts_data, dates):
    # クラスター数2-5を確認
    for n_clusters in range(2, 6):
        print('クラスター数 =', n_clusters)
  
        km = TimeSeriesKMeans(n_clusters=n_clusters, metric='euclidean', verbose=False, random_state=1).fit(ts_data)
        assert len(km.labels_) == 47
        print('クラスタリング結果 =', km.labels_)
 
        silhouette_values = silhouette_samples(ts_data, km.labels_)
        silhouette_avg = np.mean(silhouette_values)
        print('シルエット値平均 =', silhouette_avg)
 
        plot_silhouette_scores(ts_data, silhouette_values, n_clusters, km)
 
 
dfg = pd.read_csv(MOBILITY_CSV_GOOGLE)
 
# 都道府県ごとのデータを取得
dfg_jp = dfg[dfg['country_region_code']=='JP'].dropna(subset=['sub_region_1'])
 
# 都道府県のチェック
pref = dfg_jp['sub_region_1'].unique()
print(pref)
assert len(pref) == 47
 
# 日付をインデックス、都道府県をカラムとするDataFrameに変換
dfg_cat = dfg_jp.set_index(['date', 'sub_region_1'])[MOBILITY_CATEGORY].unstack()
 
# 2021年7月のデータを取得
dfg_cat_jul = dfg_cat['2021-07-01':'2021-07-31']
 
# データのプロット
ax = dfg_cat_jul.plot(figsize=(12,8), alpha=0.75, rot=90);
ax.legend(loc='center left', bbox_to_anchor=(1.0, 0.5))
plt.show()
 
# 正規化
ts_data = TimeSeriesScalerMinMax().fit_transform(dfg_cat_jul.values.transpose())
 
prefs = dfg_cat_jul.columns
dates = dfg_cat_jul.index
 
show_silhouette_scores(ts_data, dates) 

実際にはクラスター数２から５の結果が表示されるが、以下はクラスター数4の結果。平均シルエット値だけをみると他のクラスター数がよさそうだが、ここではクラスター数4を採用する。


時系列クラスタリングの結果を可視化

クラスター数が決まったので、時系列クラスタリングして、クラスターごとの平均を可視化してみる。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
%matplotlib inline
 
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
import pandas as pd
 
from tslearn.clustering import TimeSeriesKMeans, silhouette_score
from tslearn.preprocessing import TimeSeriesScalerMinMax
from tslearn.utils import to_time_series_dataset
from scipy.spatial.distance import cdist
from sklearn.metrics import silhouette_samples as sklearn_silhouette_samples
 
# Googleモビリティレポートのファイルパス
# 2020年1月3日から2月6日の5週間の曜日別中央値を基準とした値
MOBILITY_CSV_GOOGLE = 'Global_Mobility_Report210814.csv'
 
# 対象とする場所のカテゴリ
# 公共交通機関の拠点（例: 地下鉄、バス、電車の駅）など
MOBILITY_CATEGORY = 'transit_stations_percent_change_from_baseline'
 
 
## シルエット値プロット用の関数は省略 ##
 
 
def tsclustering(ts_data, prefs, n_clusters):
 
    km = TimeSeriesKMeans(n_clusters=n_clusters, metric='euclidean', verbose=False, random_state=1).fit(ts_data)
  
    # クラスタリングの結果
    result = pd.DataFrame({'clust': km.labels_}, index=prefs)
    for clust, g in result.groupby(by='clust'):
        print('Cluster: ', clust)
        print(g.index)
 
    # クラスターごとの中心をプロット
    for i, c in enumerate(km.cluster_centers_):
        plt.plot(c.T[0], label=i)
  
    plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
    plt.show()
 
 
dfg = pd.read_csv(MOBILITY_CSV_GOOGLE)
# 都道府県ごとのデータを取得
dfg_jp = dfg[dfg['country_region_code']=='JP'].dropna(subset=['sub_region_1'])
 
# 都道府県のチェック
pref = dfg_jp['sub_region_1'].unique()
print(pref)
assert len(pref) == 47
 
# 日付をインデックス、都道府県をカラムとするDataFrameに変換
dfg_cat = dfg_jp.set_index(['date', 'sub_region_1'])[MOBILITY_CATEGORY].unstack()
 
# 2021年7月のデータを取得
dfg_cat_jul = dfg_cat['2021-07-01':'2021-07-31']
 
# データのプロット
ax = dfg_cat_jul.plot(figsize=(12,8), alpha=0.75, rot=90);
ax.legend(loc='center left', bbox_to_anchor=(1.0, 0.5))
plt.show()
 
# 正規化
ts_data = TimeSeriesScalerMinMax().fit_transform(dfg_cat_jul.values.transpose())
 
prefs = dfg_cat_jul.columns
dates = dfg_cat_jul.index
 
# show_silhouette_scores(ts_data, dates)
 
n_clusters = 4
tsclustering(ts_data, prefs, n_clusters=n_clusters)

時系列クラスタリングの結果は以下の通り。4連休での差が大きい。4連休で最も移動が減少しているクラスター2には東京、大阪、愛知などの大都市圏が含まれる。逆に4連休に移動が上昇したクラスター１，３には地方の都道府県が含まれている。クラスター２に沖縄県が入っているが、これは台風で天候が悪かった影響と考えられる。

投稿者 時刻:
ラベル: , , ,

0 件のコメント:

コメントを投稿

登録: コメントの投稿 (Atom)