Rで解析：グラフ理論に基づき時系列データの時間的非対称性を評価できる「topologyR」パッケージの紹介

数値データや時系列データの背後に潜む構造を把握するには、可視化やクラスタリングだけでなく、要素同士のつながり方そのものを調べる方法があります。しかし、グラフ理論や位相空間論に基づく解析を一から実装するには手間がかかります。

「topologyR」パッケージは、時系列データから可視性グラフを構築し、Nadaら(2018)の手法に基づいて要素間のつながりを表すトポロジー（位相構造）を生成する機能を提供します。生成したトポロジーの連結性や開集合の構成を調べることが可能です。また、時系列を前向き・後ろ向きの2方向から解析し、両者の位相構造の違いから時間的な不可逆性を評価することも可能です。

本パッケージの利用で、データに潜む構造的なつながりを位相空間の枠組みで調べられるのではないかと考えます。

パッケージバージョンは0.3.0。Windows 11 x64 (build 26200)のR version 4.6.1で確認しています。

R言語徹底解説 | Hadley Wickham, 石田基広, 市川太祐, 高柳慎一, 福島真太朗

Advanced Rの邦訳。オブジェクトの仕組みや関数型プログラミングなど、Rの言語仕様を深掘りする。

パッケージのインストール

下記コマンドを実行してください。

# パッケージのインストール
install.packages("topologyR")

# パッケージの読み込み
library("topologyR")

コマンド例

詳細はコメント、パッケージのヘルプを確認してください。

本パッケージでは、データの各要素を「点」とみなし、点同士の近傍関係から「開集合」の集まりであるトポロジー（位相構造）を作ります。近傍関係の作り方には、値の近さに基づく閾値ベースの方法と、時系列の可視性グラフを使うグラフ理論的な方法があり、後者はNadaら(2018)の手法に基づいています。

以下では、閾値の候補を探る関数から始め、水平可視性グラフの構築、そのグラフからのトポロジー生成、完全グラフやアレクサンドロフ位相との比較、そして時系列を前向き・後ろ向きの2方向から見る双方向解析までを順に紹介します。

IQR係数によるトポロジー特性の変化を分析する：analyze_topology_factorsコマンド

値の近さに基づく閾値ベースの近傍関係を作る際、閾値をIQR（四分位範囲）の何分の1にするか（IQR係数）によって、生成される基底の集合数や各集合の大きさがどう変わるかを調べます。理論的に厳密なグラフベースの方法（後述のgenerate_topologyコマンド）とは異なる、メトリック近似による構成である点に注意してください。

# 乱数の種を固定
set.seed(20260826)

# 宇治市内の観測地点における、50日分の最高気温の平年差(C)を模した架空データ
data <- round(rnorm(50, mean = 0, sd = 1.5), 2)

# IQR係数を変えながらトポロジーの特徴を分析
results <- analyze_topology_factors(data, plot = FALSE)

# 結果を表示
print(results)
  factor threshold base_size max_set_size min_set_size
1      1  1.790000       201           50            0
2      2  0.895000       166           50            0
3      4  0.447500        79           50            0
4      8  0.223750        96           50            0
5     16  0.111875        57           50            0

IQR係数を大きくするほど閾値（近傍とみなす距離）は狭くなりますが、基底集合数（base_size）は単調に減るとは限らず、データのばらつき方によって増減する場合があります。いずれの基底にも、大きさ0（空集合）と50（全体集合）が必ず含まれるため、min_set_sizeは常に0、max_set_sizeは常にデータ数と一致します。

複数の指標から閾値の候補を計算する：calculate_thresholdsコマンド

閾値ベースの近傍関係を作る際の閾値の候補を、平均・中央値・標準偏差・IQR・k近傍距離という5つの異なる考え方でまとめて計算します。

# 乱数の種を固定
set.seed(20260826)

# 宇治市内の別の観測地点における、100日分の降水量偏差(mm)を模した架空データ
data <- round(rnorm(100, mean = 0, sd = 2), 2)
# 複数の方法で閾値の候補を計算

calculate_thresholds(data)
$mean_diff
[1] 0.09121212

$median_diff
[1] 0.04

$sd
[1] 1.960688

$iqr
[1] 0.6

$dbscan
[1] 0.33

時系列データから水平可視性グラフを構築する：horizontal_visibility_graphコマンド

時系列の各時点を「点」とし、途中の時点にさえぎられずに互いを見通せる（間の値がどちらよりも小さい）2時点の間にエッジを張ることで、水平可視性グラフ（HVG）を構築します。得られた隣接構造は、次のgenerate_topologyコマンドでトポロジーを生成する際の部分基底として使われます。

# 宇治市の茶畑にある観測点で記録した、8日間の日較差(最高気温-最低気温、C)
series <- c(8.2, 6.1, 9.4, 5.8, 10.3, 7.5, 6.9, 11.2)

# 水平可視性グラフを構築
g <- horizontal_visibility_graph(series)

# エッジの総数を表示
g$n_edges
[1] 11

# エッジの一覧を表示
g$edges
   from to
1     1  2
2     2  3
3     1  3
4     3  4
5     4  5
6     3  5
7     5  6
8     6  7
9     7  8
10    6  8
11    5  8

# 有向グラフとしても構築
gd <- horizontal_visibility_graph(series, directed = TRUE)

# 1日目から前方に見える観測点を表示
gd$out_adjacency[[1]]
[1] 2 3

可視性グラフの隣接構造からNadaトポロジーを生成する：generate_topologyコマンド

horizontal_visibility_graphコマンドなどで得た隣接構造（近傍関係）を部分基底とし、共通部分による閉包で基底を求め、特殊化前順序に基づいて連結性を厳密に判定します。閾値ベースの方法と異なり、値の大小関係そのものから構造を捉える、理論的に厳密な構成です。

# 嵐山地区の観測点における、5日間の日較差(C)
series2 <- c(9.4, 5.8, 10.3, 7.5, 6.9)

# 水平可視性グラフの隣接構造を取得
g2 <- horizontal_visibility_graph(series2)

# 隣接構造からNadaトポロジーを生成
topo <- generate_topology(g2$adjacency, g2$n)
Warning message:
Enumerating topology of a space already determined to be disconnected (2 components). Consider using max_open_sets = 0 to skip enumeration when only connectivity is needed.

# 連結かどうかを確認
topo$connected
[1] FALSE

# 連結成分の一覧を表示
topo$components
[[1]]
[1] 1 2 3

[[2]]
[1] 4 5

この例では、1〜3日目のグループと4〜5日目のグループに分かれ、位相的には連結していないという結果になりました。

Rによるデータサイエンス(第2版):データ解析の基礎から最新手法まで | 金明哲

主成分分析やクラスター分析から、サポートベクターマシン・ニューラルネットワーク・生存分析まで、幅広い手法をRで一望できる。

＜おすすめのRに関する書籍です＞

基礎から学ぶ統計学実践編〜PCを使った計算とカテゴリカルデータ解析

PCで実際に計算しながら進める実践編。クロス集計などカテゴリカルデータの解析を扱う。

www.amazon.co.jp

完全グラフから自明なトポロジーを構築する：complete_topologyコマンド

すべての点が互いに隣接する完全グラフから、generate_topologyコマンドと同じ手続きでトポロジーを構築します。完全グラフでは各点の近傍が「自分以外の全要素」という自明な集合になるため、生成されるトポロジーは値の大小に関係なく常に自明な（空集合と全体集合だけからなる）ものになります。値そのものの構造を調べたい場合は、可視性グラフを使うgenerate_topologyコマンドを使ってください。

# 京都市内5地点で観測した降水量(mm)
data <- c(12.4, 8.9, 15.2, 3.1, 9.7)

# 完全グラフからトポロジーを構築
result <- complete_topology(data)

# 開集合の数を表示
result$n_open_sets
[1] 2

# 連結かどうかを確認
result$connected
[1] TRUE

有向可視性グラフからアレクサンドロフ・トポロジーを生成する：generate_alexandrov_topologyコマンド

有向グラフの到達可能性（ある点から辿り着ける点の集合）を開集合とみなす、アレクサンドロフ・トポロジーを生成します。可視性グラフでは連続する観測点同士は必ず辺で結ばれるため、到達可能性は時刻順の全順序となり、生成されるアレクサンドロフ・トポロジーは常に連結し、基底の要素数は常にデータ数nに一致します。

# 宇治市の茶畑の観測点における、8日間の日較差(C)（水平可視性グラフの節と同じ系列）
series <- c(8.2, 6.1, 9.4, 5.8, 10.3, 7.5, 6.9, 11.2)

# 有向の水平可視性グラフを構築
gd <- horizontal_visibility_graph(series, directed = TRUE)

# 到達可能性に基づくアレクサンドロフ・トポロジーを生成
alex <- generate_alexandrov_topology(gd$out_adjacency, gd$n)

# 連結かどうかを確認
alex$connected
[1] TRUE

# 基底集合の数を表示
length(alex$base)
[1] 8

トポロジーの連結性を素朴な方法で判定する：is_topology_connectedコマンド

開集合の一覧（トポロジー）を渡し、同じ開集合に含まれる要素同士を辺でつないだグラフが連結かどうかを判定します。generate_topologyコマンドの厳密な判定（特殊化前順序による判定）とは異なり、必要条件ではあっても十分条件ではない簡易な判定である点に注意してください。厳密な判定にはis_topology_connected_exactコマンドを使います。

# 宇治茶の産地(1:宇治市 2:宇治田原町 3:和束町 4:南山城村 5:笠置町)を
# 出荷組合ごとにグループ分けした結果を、重なりのある開集合として表現
topology <- list(c(1, 2), c(2, 3, 4), c(4, 5))

# グラフとしての連結性を判定
is_topology_connected(topology)
[1] TRUE

双方向の可視性グラフから時間的不可逆性を評価する：generate_bitopologyコマンド

時系列から前向き・後ろ向き2種類の有向可視性グラフを作り、それぞれについてNadaトポロジー（前向きトポロジーτ+、後ろ向きトポロジーτ-）を生成したうえで、両者の違いを数値化します。可逆な過程であればτ+とτ-はよく似た構造になりますが、非対称な過程（例えば緩やかな上昇と急な下降を繰り返すような過程）では両者が乖離します。戻り値のinvariantsに、bitopology_invariantsコマンドが計算する非可逆性の指標がまとまっています。

まず、水平可視性グラフの節で使った8日分の日較差データで、戻り値の構造を確認します。

# 宇治市の茶畑の観測点における、8日間の日較差(C)
series <- c(8.2, 6.1, 9.4, 5.8, 10.3, 7.5, 6.9, 11.2)

# 双方向のトポロジー解析を実行
bt <- generate_bitopology(series)

# 前向きトポロジーの連結成分数を表示
bt$invariants$forward_components
[1] 3

# 後ろ向きトポロジーの連結成分数を表示
bt$invariants$backward_components
[1] 3

# 成分数に基づく非可逆性の指標を表示
bt$invariants$irreversibility_components
[1] 0

前向き・後ろ向きの連結成分数が3で一致しており、irreversibility_componentsも0になりました。この8日分のデータには、少なくとも連結成分数で見る限り、明確な時間的非対称性は見られません。次に、緩やかに増加していく傾向を持つ、より長い系列で確認します。

# 乱数の種を固定
set.seed(20260826)

# ある観光地点の来訪者数が緩やかに増加していく様子を模した、40日分の架空の指標
n <- 40

# 全体として緩やかに上昇するトレンドを作成
trend <- seq(0, 4, length.out = n)

# 平均0、標準偏差0.6の正規分布ノイズを重ねる
noise <- rnorm(n, mean = 0, sd = 0.6)
series_trend <- round(trend + noise, 2)

# 双方向のトポロジー解析を実行
bt2 <- generate_bitopology(series_trend)

# 前向き・後ろ向きの連結成分数を表示
bt2$invariants$forward_components
[1] 23

bt2$invariants$backward_components
[1] 21

# 成分数に基づく非可逆性の指標を表示
bt2$invariants$irreversibility_components
[1] 0.08695652

# 非対称性の向き(正なら前向きの方がより連結)を表示
bt2$invariants$asymmetry_direction
[1] -2

この例では前向きトポロジーの連結成分数（23）の方が後ろ向き（21）よりわずかに多く、asymmetry_directionは-2となりました。これは後ろ向きトポロジーの方がわずかに連結性が高いことを示しており、必ずしも「なだらかな上昇・急な下降」という典型例と同じ向きの非対称性が出るとは限らない、という点も含めて、実際のデータで試してみる価値があります。

Rではじめるデータサイエンス 第2版 | Hadley Wickham, Mine Çetinkaya-Rundel, Garrett Grolemund, 大橋 真也

tidyverseの事実上の公式教科書。データの取り込み・整形・可視化・モデリング・レポート化までを一本の流れとして学べる。

この記事が誰かの役に立ちますように。