PythonでRの`cld`関数みたいなことをする

時間のない方は，[結果](#結果)だけお読みください．

# 要約

Compact letter display (CLD)は，多重検定による多対多の有意関係を，文字列で表すことで，有意関係を少ないスペースで示すことができる．
任意の有意関係からCLDを導出することは処理区が増えてくると困難になり，自動化が求められる．
R言語では，`cld`関数を用いることでCLDの導出が可能だが，Pythonではそのようなパッケージが存在しなかった．
本記事では，グラフ理論を用いることで，任意の有意関係を示す行列から，どの処理区にどの文字列を対応させるかを導出する手法を示す．
これにより，Pythonによるデータビジュアリゼーションをより容易に行うことができる．

# イントロダクション

## 処理区が増えるとCLDを導出するのが困難になる

[Tukey-Kramer法](https://ja.wikipedia.org/wiki/テューキーの範囲検定)などを用いて多重検定したとき，3~4つくらいの処理区ならば，有意差のある対応関係を図示することは容易である．

たとえば，処理区A~Dに関して多重検定を行い，以下の表に示した*p*値が得られたとすると，

|   *p* |    A   |    B   |    C   |    D   |
|------:|:------:|:------:|:------:|:------:|
| **A** |   -1   | 0.7238 | 0.0427 | 0.0005 |
| **B** | 0.7238 |   -1   | 0.2550 | 0.0081 |
| **C** | 0.0427 | 0.2550 |   -1   | 0.9574 |
| **D** |  0.005 | 0.0081 | 0.9574 |   -1   |

有意差のあった対応関係は，以下のように図示できる(\*, *p* < 0.05; \*\*, *p* < 0.01; \*\*\*, *p* < 0.001)．
<img width="238" alt="Untitled.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/214914/20337d20-11ca-2985-1815-383be720907d.png">

しかし，処理区が5~6つを越えてくると，有意差のある処理区の関係性が複雑になり，上記に示したような，線による図示が困難となる．
そういう場合には，有意差の有無を文字によって示すことがある．
さきほどの図であれば，
<img width="238" alt="Untitled_2.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/214914/14024cab-5a3d-8d6b-2974-12e1c2615224.png">
のように示し，同じ文字を含む処理区では有意差なし(*p* >= 0.05)を意味する．

この方法は，[Letter-based display of all-pairwise comparison](https://www.google.co.jp/search?q=Letter-based%20display%20of%20all-pairwise%20comparison)や[Compact letter display](https://www.google.co.jp/search?q=compact%20letter%20display) (CLD)などと呼ばれ，複雑な有意差の関係性を少ないスペースで示すことができる．
しかし，CLDを有意関係から取得するのは非常に難しい．
たとえば，以下の表では，10の処理区に関して多重検定を行い，*p* < 0.05ならば1，*p* >= 0.05ならば0として示したが，この表からCLDを矛盾なく導くのは大変困難である．

|       |  A |  B |  C |  D |  E |  F |  G |  H |  I |  J |
|------:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|
| **A** | -1 |  1 |  0 |  0 |  0 |  0 |  1 |  1 |  0 |  1 |
| **B** |  1 | -1 |  0 |  1 |  0 |  0 |  0 |  0 |  0 |  0 |
| **C** |  0 |  0 | -1 |  0 |  1 |  0 |  0 |  0 |  1 |  0 |
| **D** |  0 |  1 |  0 | -1 |  1 |  1 |  0 |  0 |  0 |  1 |
| **E** |  0 |  0 |  1 |  1 | -1 |  0 |  0 |  0 |  0 |  0 |
| **F** |  0 |  0 |  0 |  1 |  0 | -1 |  1 |  0 |  0 |  0 |
| **G** |  1 |  0 |  0 |  0 |  0 |  1 | -1 |  1 |  1 |  0 |
| **H** |  1 |  0 |  0 |  0 |  0 |  0 |  1 | -1 |  0 |  0 |
| **I** |  0 |  0 |  1 |  0 |  0 |  0 |  1 |  0 | -1 |  0 |
| **J** |  1 |  0 |  0 |  1 |  0 |  0 |  0 |  0 |  0 | -1 |

ちなみに，以下の表のようなCLDが，有意関係に矛盾が無い．

|       |   CLD   |
|------:|:-------:|
| **A** |   adej  |
| **B** |   bchi  |
| **C** |  defghi |
| **D** |  dfgjk  |
| **E** |   abc   |
| **F** |   abei  |
| **G** |   cfh   |
| **H** |   bgik  |
| **I** |   abjk  |
| **J** |   bchi  |

R言語では，[`cld`関数](https://www.rdocumentation.org/packages/lsmeans/versions/2.27-62/topics/cld)を用いることで，CLDを導出することができる．
R言語の`cld`のドキュメントでは，Piephoら(2004)[^piepho]の論文が引用されており，当該論文の手法でも矛盾のないCLDが可能である．
しかし，Piephoらの方法では冗長な文字列が生じ，Piephoらは，必要十分なCLDを可能にするアルゴリズムを示していない．
また，Piephoらの方法は，計算にかかるコストが高く，より高速で簡便なアルゴリズムが必要である．

## グラフ理論によってCLDが導出できる

Pythonには，[scikit-posthocs](https://github.com/maximtrp/scikit-posthocs)のような，優秀な多重検定パッケージが存在するが，CLDの導出には対応していない．
本記事では，[グラフ理論](https://ja.wikipedia.org/wiki/グラフ理論)を用いたCLDの導出法を示す．

上記の10の処理区からなる有意関係は，有意差のない処理区同士がつながった無向グラフであるとみなすことができる．
<img width="353" alt="Unknown.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/214914/c07ce5e4-8af9-7d73-453f-804a6f11af1b.png">

有意関係をグラフで表現したとき，次のことが分かる:

1. 任意の2つの処理区に関して，距離が2以上のとき，その2つの処理区間は有意差がある
2. 1から，同じ[クリーク](https://ja.wikipedia.org/wiki/クリーク_(グラフ理論))に属する処理区の間では有意差が無い
3. 2から，[クリーク辺被覆問題 (Clique Edge Cover Problem)](https://www.google.co.jp/search?q=clique%20edge%20cover%20problem)を解き，得られたクリーク集合の各要素に含まれる処理区に同一の文字列を与えることで，必要最小限の文字列からなるCLDが得られる．

よって，
1. 有意差の無い処理区をつなげる隣接行列を得る
2. 1からグラフを作る
3. 全辺を被覆するクリークの集合を得る
4. 各クリークに対して`a`から始まる文字列を決める
5. 各クリークに含まれる全ての処理区に，4で決めた文字列を加える
というアルゴリズムで，CLDを得ることができると考えられる．

# 結果

有意関係の表から，有意差のない関係を1とする行列を作成すると，有意差のない処理区がつながった無向グラフを表現する[隣接行列](https://ja.wikipedia.org/wiki/隣接行列)が得られる．
この隣接行列から，`networkx`を用いてグラフを作成し，`find_cliques()`によってグラフ被覆に十分なクリークを得る．
このクリークごとに文字列を割り振ることで，必要最小限の文字列からなるCLDを導出できる．

```python
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt

significance = np.matrix(
    [[-1,  1,  0,  0,  0,  0,  1,  1,  0,  1],
     [ 1, -1,  0,  1,  0,  0,  0,  0,  0,  0],
     [ 0,  0, -1,  0,  1,  0,  0,  0,  1,  0],
     [ 0,  1,  0, -1,  1,  1,  0,  0,  0,  1],
     [ 0,  0,  1,  1, -1,  0,  0,  0,  0,  0],
     [ 0,  0,  0,  1,  0, -1,  1,  0,  0,  0],
     [ 1,  0,  0,  0,  0,  1, -1,  1,  1,  0],
     [ 1,  0,  0,  0,  0,  0,  1, -1,  0,  0],
     [ 0,  0,  1,  0,  0,  0,  1,  0, -1,  0],
     [ 1,  0,  0,  1,  0,  0,  0,  0,  0, -1]])

labels = ('A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J',)

# 有意差の無いもの(0)を1とする隣接行列を作る
mat = (significance == 0)

V = nx.from_numpy_matrix(mat)
G = nx.Graph(V)

# クリーク辺被覆問題を解き，同じ文字を与えるべきクリーク集合を得る
cliques = nx.find_cliques(G)

# 各クリーク集合の要素について，同一の文字列ラベルを与える
cld_dict = {i: [] for i in labels}
for i, clique in enumerate(cliques):
    letter = chr(i + 97)
    for j in clique:
        cld_dict[j].append(letter)

for i in cld_dict:
    cld_dict[i] = ''.join(cld_dict[i])

print(cld_dict)

# {'A': 'bejk'
#  'B': 'fghi'
#  'C': 'cdegik'
#  'D': 'abcde'
#  'E': 'fhj'
#  'F': 'fgjk'
#  'G': 'chi'
#  'H': 'adfg'
#  'I': 'abfj'
#  'J': 'fghi'}
```

[^piepho]: [Piepho, H.P. (2004) "An algorithm for a letter-based representation of all pairwise comparisons" *Journal of Computational and Graphical Statistics*, **13**, 2:456-466](https://www.tandfonline.com/doi/abs/10.1198/1061860043515)