R qgraphでspring-embeddedレイアウトのグラフを描く

[R](https://www.r-project.org/)の[qgraph](http://sachaepskamp.com/qgraph/reference/qgraph.html)を使って、下図のようなspring-embeddedレイアウトのネットワーク図を書く方法の備忘録です。
![qgraph-spearman_sample_w400.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/469745/aab5bf2b-9d45-2b1d-92b1-babe29debed7.png)
題材は、自分がこれまでコーディングしたCSSファイルを[analyze-css](https://github.com/macbre/analyze-css)で解析したメトリクスデータを集めたもの。メトリクスの相関の強さをqgraphで表してみました。

**動作環境および必要なソフトウェア**
バージョンは本記事執筆時点で使用しているもの

* macOS Mojave 10.14.6
* analyze-css 0.12.7 [https://github.com/macbre/analyze-css](https://github.com/macbre/analyze-css)
* Node.js v12.13.0 [https://nodejs.org/en/](https://nodejs.org/en/)
* R 3.6.1 [https://www.r-project.org/](https://www.r-project.org/)
* qgraph 1.6.3 [http://sachaepskamp.com/qgraph/](http://sachaepskamp.com/qgraph/)

## 第1章　分析用のデータを集める
[analyze-css](https://github.com/macbre/analyze-css)は、CSSファイルを解析してルール記述やセレクタの個数などの情報を出力してくれるコマンドラインツールです。ver.0.12.7現在、37種類の[メトリクス](https://github.com/macbre/analyze-css#metrics)がサポートされています。
### analyze-cssの使い方
npmでインストールします。

```bash
$ npm install --global analyze-css
```

計測したいCSSファイルのパス名を引数に指定して`analyze-css`を実行すると、解析結果がJSON形式で標準出力に書き出されるので、リダイレクトして適当なファイルに保存します。引数`--pretty`を指定すると、項目ごとに改行とインデントが付いて読みやすい形式になります。

```bash
$ analyze-css --pretty --file examples/style.css > style.json
```

### JSONをCSV形式に変換
Rで処理しやすくするため、ファイルのフォーマットをJSON形式からCSV形式に変換します。pythonを使うと簡単です。抽出する項目はキー"metrics"が保持するデータだけです。

```python:css-metrics.py
import json
import sys

args = sys.argv
if len(args) < 2:
	sys.exit('usage: python %s json_file' % args[0])

with open(args[1], 'r') as f:
	jsn = json.load(f)
	for key in jsn['metrics']:
		print(key + ',' + str(jsn['metrics'][key]))
```

jsonファイルを１つのフォルダに集め、一括してCSV形式に変換します。

```bash
$ ls json/
project1.json	project2.json	project3.json
```

以下のシェルスクリプトを実行すると、同じフォルダに拡張子.csvのファイルが出来ます。

```bash:mkcsv.sh
#!/bin/bash

for i in $(find json -name '*.json' -print)
do
  echo $i
  outf=${i%.*}.csv
  echo "metrics,value" > $outf
  python css-metrics.py $i | sort >> $outf
done
```

```bash
$ ls -p
css-metrics.py	json/		mkcsv.sh
$ ./mkcsv.sh 
$ ls json
project1.csv	project2.csv	project3.csv
project1.json	project2.json	project3.json
```

バラバラのCSVファイルを１個のCSVファイルに統合します。その際、１行目に列見出し、2行目以降にデータ行が追加されるようにしたいので、Pythonを使って元のデータ（metrics,valueの２列構成）の行と列を転置して出力します。また、１列目には、ファイル名から抽出したプロジェクト名を挿入し、見出しをprojectに書き換えます。

```python:csv-transpose.py
import csv
import sys
import os

hdr = False
args = sys.argv
usage = 'usage: python %s csv_file [ H ]' % args[0]
if len(args) < 2:
	sys.exit(usage)
if len(args) > 2:
	if args[2] == 'H':
		hdr = True
	else:
		sys.exit(usage)
	 
basename = os.path.basename(args[1])
tpl = os.path.splitext(basename)
projname = tpl[0]
key_list = []
val_list = []

with open(args[1], 'r') as f:
	reader = csv.reader(f)
	header = next(reader)
	key_list.append('project')
	val_list.append(projname);
	for row in reader:
		key_list.append(row[0])
		val_list.append(row[1])
	if hdr == True:
		print(','.join(key_list))
	print(','.join(val_list))
```

以下のシェルスクリプトを走らせてフォルダ内のCSVファイルを一括して変換し、結果の標準出力をファイルに保存します。

```bash:merge-list.sh
#!/bin/bash

hdr=H
for i in $(find json -name '*.csv' -print)
do
  python csv-transpose.py $i $hdr
  hdr=
done
```

```bash
$ ls -p
css-metrics.py		json/			mkcsv.sh
csv-transpose.py	merge-list.sh
$ ./merge-list.sh > css-metrics-data.csv
```

保存したCSVファイルを表計算アプリで開いて、正しく変換できたか確認します。これでRに読み込ませるデータが準備できました！
![csv-convert_w1000.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/469745/840682c4-3540-37ef-4a10-0ac06e43897e.png)

## 第2章　Rにデータを読み込んでqgraphでグラフを描く
以下、**[R]**で始まる見出しは「RのTips」、**[qgraph]**で始まる見出しは「qgraphのTips」のつもりで書いてあります。Tipsとして読めるよう、分析の詳細については（退屈ですし）端折ってます。
### [R] データフレームの作成： CSVファイルの読み込み
Rを起動し、Rコンソールから以下のコマンドを入力して、前節で作ったCSVファイルをデータフレームに読み込ませます。

```R
df <- read.table("css-metrics-data.csv", header=T, sep=",")
```

### [R] データフレームから列を削除（その１）
データの相関関係の分析には[スピアマンの順位相関係数](https://en.wikipedia.org/wiki/Spearman%27s_rank_correlation_coefficient)を使います。dfには、相関係数の算出に向かない列も含まれているため、事前に取り除いておきます。たとえば以下のような列は不要なので削除します。

1. １列目の見出し（プロジェクト名）
1. 値が全部ゼロ
1. 値が全部同じ（ゼロ以外）
1. 値が0と1しかない（値のばらつきがほとんどない）

項目1〜3については以下のコマンドで対応します。`df[,-1]`は１列目を除く列、`colSums(df) != 0`は列の和が0でない（オール0でない※）列のみTRUE、`apply(df, 2, var) != 0`は分散が0でない（値が均一でない）列のみTRUEになることを利用して、`df`の該当列のみ抽出しています。
※測定値が0以上であることが前提

```R
df <- df[,-1]
df <- df[colSums(df) != 0]
df <- df[apply(df, 2, var) != 0]
```

項目4については、以下のコードのように関数`is.bin()`を定義し、引数に与えたベクトル`x`の要素の集合(重複なし）が{0, 1}に等しいかどうかを判定します。データフレーム `df`から`is.bin()`がTRUEになる列を除外して、`df`を置き換えます。

```R
is.bin <- function(x) {
   return(setequal(unique(x), c(0,1)))
}
df <- df[!apply(df, 2, is.bin)]
```

### [R] スピアマンの順位相関係数の算出
相関係数の算出には関数`cor()`を使います。`df`のすべての列の組み合わせに対して、スピアマンの順位相関係数を計算した結果を変数`v`に保存します。

```R
v <- cor(df, method="spearman")
```

### [qgraph] springレイアウトのグラフ作成
ここまでの状態で一旦グラフを描いてみましょう。まずqgraphパッケージをインストールします。依存するパッケージが芋づる式に多数インストールされるため少々時間がかかります。

```R
install.packages("qgraph")
```

qgraphパッケージをロードしてみて、正しくインストールできたことを確認します。

```R
> library(qgraph)
Registered S3 methods overwritten by 'huge':
  method    from   
  plot.sim  BDgraph
  print.sim BDgraph
```

qgraphパッケージをロードした後、関数`qgraph()`に順位相関係数のベクトル`v`を渡して、springレイアウトでグラフを描かせます（下図）。

```R
qgraph(v, minimum=0.7, layout="spring")
```

![qgraph-spearman_01_w1000.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/469745/0b01c517-8347-27b6-dd89-5bdb96e9a9c1.png)
グラフの各ノードが、変数`v`の各要素つまりメトリクスを表します。ノードのラベルはメトリクス名を３文字に省略した文字列です。
例：sBC = selectorsByClass
ラベル名と元のメトリクス名をグラフの凡例に表示する方法は後述します。

ノードどうしを結ぶ辺は相関の強さを表します。線が太く・色が濃いほど相関が強いことを示しています。`qgraph()`は、デフォルトでは全ノードのすべての組み合わせに対する相関関係の辺を描画しますが、変数が多くビジーすぎるので、`qgraph()`に引数`minimum=0.7`を指定することにより、相関係数が0.7以上の辺のみをグラフに表示させました。

### [R] データフレームから列を範囲で抽出
sBT (selectorsByTag) とsTT (specificityTagTotal) など、CSSセレクタの[詳細度](https://specificity.keegan.st/)に関するメトリクスを表すノードと、セレクタの個数を表すノードはどちらも同じ要素を計数しているため、両者の増減が連動し相関が強く出るのは当然なので、両方をグラフに入れると冗長だと判りました。そこで、データフレームから詳細度に関するメトリクスspecificity* の列を削除します。

データフレームの列名を確認するには関数`colnames()`を使います。

```R
> colnames(df)
 [1] "colors"                       "comments"                     "commentsLength"                   
 .....    
[19] "selectorsById"                "selectorsByPseudo"            "selectorsByTag"              
[22] "specificityClassAvg"          "specificityClassTotal"        "specificityIdAvg"            
[25] "specificityIdTotal"           "specificityTagAvg"            "specificityTagTotal"     
```

specificity* の列はインデックス22以降に集まってますので、列1〜21までを抽出して`df`を書き換えます。

```R
df <- df[1:21]
```

### [R] データフレームから列を削除（その２）
このほか、ノードimp (imports) と emR (emptyRules) から出ている辺も、元のデータの偏りから生じたノイズだと判ったため（CSSの作りからして本来無関係）、これらもデータフレームから削除することにします。

データフレームから列を名指しで除外するには、列名の集合に対して関数`setdiff()`を利用します。

```R
df <- df[setdiff(colnames(df), c("emptyRules", "imports"))]
```

### [qgraph] グラフの書式を調整
グラフに凡例を追加し、ノード間の辺に相関係数の値を表示し、文字の大きさを調整します。

```R
v <- cor(df, method="spearman")
qgraph(v, minimum=0.7, layout="spring", repulsion=0.6,
    legend=T, nodeNames=names(v[,1]), legend.cex=0.5,
    edge.labels=T, edge.label.cex=0.5,
    label.cex=0.8)
```
![qgraph-spearman_04_final_w1000.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/469745/a4541aef-a6ad-76d0-ac5d-554f3623fa76.png)

`qgraph()`の引数に追加した項目は以下のとおりです。

`repulsion=0.6`
1より小さくすると、ノードの“島”どうしの間隔が近くなり、個々の島の占める領域が広がります。ノードが密集している部分のノード間の辺が長く表示されて、結合状態がよく見えるようになります。

`legend=T, nodeNames=names(v[,1]), legend.cex=0.5`
グラフの右側に凡例を表示します。`nodeNames`は、ノードのラベルに対するフルネームのベクトルを指定します。この例では、関数`cor()`の戻り値を保存した変数`v`に入っている列名のリストを`v[,1]`により取り出しています。凡例の文字列の拡大率を0.5に設定します。

`edge.labels=T, edge.label.cex=0.5`
辺にラベル(相関係数の値）を表示し、ラベル文字列の拡大率を0.5に設定します。

`label.cex=0.8`
ノードのラベル文字列の拡大率を0.8に設定します。

## [qgraph] ノードをグループ化して色分けする
グラフのノードを以下の４つのグループに分けて、それぞれノードに色をつけてみましょう。

| グループ | 色 | 意味 | ノード |
|:-:|:-:|:--|:--|
| A  | `#FF4081` | CSSの保守性に関わる指標   | lng, sLA  |
| B  | `#A1887F` | CSSの構成要素の計数 | clr, dcl, mdQ, rls, slc, sBC, sBI, sBP, sBT   |
| C  | `#2196F3` | 複雑な構造の計数| cmS, dpP, dpS, mcS, qlS, rCN  |
| D  | `#BDBDBD` | その他 | cmm, cmL  |

まず、`qgraph()`の引数に渡すデータを準備します。ノードのラベルに対応する列名を`colnames(df)`から調べ、そのインデックスをノード番号として`c()`を使ってグループごとのベクタを作成します。4個のグループに対するベクタを要素とするリストを`list()`を使って作成し、変数`groups`に保存します。

グループの配色は、RGB値`#...`や色名の配列で指定し、変数`color`に保存します。配列の要素1が`groups`リストの要素1に対応します。

```R
groups <- list(c(8,14), c(1,5,9,13,15,16,17,18,19), c(4,6,7,10,11,12), c(2,3))
color <- c("#FF4081", "#A1887F", "#2196F3", "#BDBDBD")
```

`qgraph()`に引数`groups=groups`、`color=color`を渡すとグラフのノードが塗り分けられます。以下のコードでは`label.color`でラベルの文字列の色を指定し、`borders=F`でノードの境界線を非表示にしました。ただし、凡例左側の小さな円の境界線は`borders=F`では消えません。対策として、`border.color="white"`も指定して凡例側の円を若干見やすくしました。

```R
qgraph(v, minimum=0.7, layout="spring", repulsion=0.6,
       legend=T, nodeNames=names(v[,1]), legend.cex=0.5,
       edge.labels=T, edge.label.cex=0.5, label.cex=0.8,
       groups=groups, color=color,
       label.color="white", borders=F, border.color="white")
```
![qgraph-spearman_05_groups_w1000.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/469745/8183a961-6ef0-9108-4b51-0c38b78e29b4.png)
これでグラフが完成しました！

`#FF4081`のノードと`#2196F3`のノードの関係を見て、私が書いていたCSSはセレクタの指定が複雑になる傾向にあることが判りました。確かに、`analyze-css`の出力jsonファイルには、`.foo .bar .baz p`みたいにdivの階層辿るようなセレクタを使ってる箇所が多く報告されていました[^1]。このようなセレクタは詳細度が無駄に上がってしまい、media queryの中などでCSSプロパティを上書きしたい時、セレクタを元と同じ形式かそれ以上の詳細度で指定しなければならず、メンテナンスしにくいCSSになる要因ですね。
[^1]: コーディングを始めた頃、既存サイトに下層ページを追加する案件を担当していて、手本にした既存のCSSが`#product .foo ul.bar li {...}`みたいに冗長なセレクタを多数使用していたのに倣ったのがいけなかった。

## [R] 箱ひげ図が付いた散布図の描き方
おまけとして、散布図のX軸とY軸の外側に箱ひげ図が付いた下図のようなグラフをRで描く方法です。このグラフは、CSSのlng (length)  と sLA (selectorLengthAvg) のデータの分布をプロットしたものです。２本の赤い線はそれぞれlngとsLAの中央値を示します。
![scatterplot_06_w450.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/469745/277be7b4-25f7-dbd9-330f-a728fff5d96e.png)

箱ひげ図が付いた散布図は、carパッケージの関数`scatterplot()`で簡単に描けます。

```R
install.packages("car")
```

以下のコードではcarパッケージをロード後、`par()`でグラフィックスパラメータを退避・更新してから、`scatterplot()`を呼び出します。さらに、関数`abline()`でX, Y軸に平行な直線を引き、最後にグラフィックパラメータを元に戻しています。グラフィックスパラメータを更新した理由は、箱ひげ図の外側に余白を入れたいので、`par(oma=...)`を呼んで作図領域の底辺と左辺にマージンを設定するためです。

```R
library(car)

oldpar <- par(no.readonly=T)
oma <- c(2,2,0,0)
par(oma=oma)

scatterplot(df$length, df$selectorLengthAvg,
	smooth=F, regLine=F, grid=F,
	xlab="length", ylab="selectorLengthAvg",
	pch=16, cex=1.5, reset.par=F)

abline(v=median(df$length), col="red")
abline(h=median(df$selectorLengthAvg), col="red")

par(oldpar)
```

lngとsLAが共に大きいグループ＝赤い線で区切られた右上の領域＝にプロットされたプロジェクトは、ページ制作当時、ハマりまくったものばかりなので、生産性の悪さ加減がグラフに現れてるなと納得しました（苦笑）