Stataでシークエンス分析
Sequence Analysis のやり方を日本語で紹介したものがみあたらないので、Stataのパッケージマニュアルを参考に、紹介してみる。
Sequence Analysisとは？
Seaquence Analysisとは、並び順のパターンを発見する分析である。つまり、
11111
11311
11122
77777
のように4つの並び順があるとすると、上二つはかなり似ていて、四つ目はかなり違う。三番目はそこそこ似ている。みたいなのを発見する分析。手順としては以下の二段階を踏む。
- 配列同士の距離を計算する
- 距離からクラスター分析でいくつかのパターンに分類する
今回はStataを用いてSAを行う。
データの構築とパッケージの導入
データの構築
まずはパッケージの導入からデータの加工まで。とりあえず、定型データ（全てのサンプルが同じ長さ）形式にする^1。つまり、
11111
11311
221
のように３つ目だけ短い、みたいなのを避けて全部同じ長さのデータを作る。ひとまず、３つ目の空白に9とか代入して以下の形へ。
11111
11311
22199
パッケージの導入
なんはなくともパッケージの導入。
*** パッケージ導入 ssc install sq
データはロング型で与える。つまり、以下みたいな形でズラーッと長く^2。
| ID | order | みたい変数 |
|---|---|---|
| 1 | 1 | 1 |
| 1 | 2 | 1 |
| 1 | 3 | 1 |
| 1 | 4 | 1 |
| 1 | 5 | 1 |
| 2 | 1 | 1 |
| 2 | 2 | 1 |
| 2 | 3 | 3 |
| … | … | … |
つづいて、配列が入ってる変数と、各配列のid、配列の順序orderを指定して、sqsetの形を指定する（xtsetやstsetみたいにidや順序を指定したデータセット)。
*** データ形式指定
sqset みたい変数 id order, ltrim
分析
基礎分析
まずは基礎的な分析を行う。いわゆる記述統計的な分析。
| 関数 | 意味 | 例 |
|---|---|---|
| sqdes | 全要素をばらしたときの頻度を表しているはず^3 | |
| sqtab | 出現配列を出力。オプションは、rankは上位のみを出力。soは例えば11112と12222を同じものとしてカウント。seは12121と11112を同じものとしてカウント。 | |
| sqlength | 配列の長さ | 12121であれば5である。elementオプションは |
| sqelemcount | 配列内で異なる要素の数 | 11111であれば1、12121であれば2。 |
| sqepicount | 異なるエピソードの数。つまり数値の変化が何回起こったか-1 | 11111であれば1、12121であれば5。 |
*** 要素の頻度（長さ1の要素の頻度） sqdes *** 出現配列一覧 sqtab * top10のみ表示 sqtab,ranks(1/10) * 長さを考慮せずどういう順列か（ same order） sqtab,ranks(1/10) so * パターン構成要素が同じ配列の頻度（順不同でいかなる要素によって構成されるか same elements） sqtab,ranks(1/10) se * シークエンスの長さ egen length = sqlength() egen length1 = sqlength(),element(2) * 異なる要素の数 egen elemnum = sqelemcount() * エピソード単位のカウント egen epinum = sqepicount() * リスト表示で変数一覧 sqstatlist * 記述統計 sqstatsum * 1変数 sqstattab1 * クロス表 sqstattab2
図示
シークエンス分析の強みは何と言っても図示。
*** Sequence index plots
sqindexplot
シークエンス分析
以降では、いよいよ配列間の距離を計算しクラスタリングする。ここで重要なのは以下の点。
- 配列間の距離をいかに決定するか
- いかにクラスタリングするか
距離を求める
まず、配列間の距離は Levenshtein距離 を用いる。Levenshtein距離は、 indel-cost と substitution cost の2つの合計である。indel-cost(insert delete)とは、配列のある要素を挿入するないし、消すことである。例えば1231→131(消去)や、1231→12231(挿入)を指す。一方 substitution cost はある要素を他の要素に変更することである。例えば1231→1241（3を4に変更）。これらの操作それぞれにコストを決め（例えばindel-costを1、substitution costを2）、「コスト × 行われた手順」が Levenshtein距離となる。
ただし、これらの手順が一通りとは限らない。そこで、最小距離を求めるアルゴリズムが Needleman-Wunsch アルゴリズムである。
まずsqomで距離を求める。全ての距離を計算するオプションがfull。結構時間がかかる。名前を指定しなければSQDistというマトリックスに距離行列がしまわれる。
よくわからないが、データの形式上 sqclusterdat をして特殊な画面？に移動する。そこでクラスター分析を行い、sqclusterdat,returnで元の画面に戻って来る。クラスター分析の際、何クラスで割り当てるかはgroupの中で指定する。indexplotをクラスターごとに出力したところで今日はここまで。
*** 距離を計算（OM） sqom, full *** クラスター分析 sqclusterdat * 階層的クラスター分析（今回はウォード法で行う） * クラスターを計算 clustermat wardslinkage SQdist, name(wards) add * デンドログラムの出力 cluster tree wards, cutnumber(20) * カテゴリー値の割り当て cluster generate cl3=group(3), name(wards) cluster generate cl4=group(4), name(wards) cluster generate cl5=group(5), name(wards) cluster generate cl6=group(6), name(wards) sqclusterdat, return sqindexplot, by(cl3)
