专栏名称: 连享会

连玉君老师团队分享，主页：lianxh.cn。白话计量，代码实操；学术路上，与君同行。

Sequence Analysis：生命历程中的序列分析

连享会 · 公众号 · · 2025-02-08 22:00

正文

👇 连享会 · 推文导航 | www.lianxh.cn

🍎 Stata： Stata基础 | Stata绘图 | Stata程序 | Stata新命令
📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享
💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归
⛳ 专题： SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析
🔃 因果： DID | RDD | 因果推断 | 合成控制法 | PSM-Matching
🔨 工具：工具软件 | Markdown | Python-R-Stata
🎧 课程：最新专题 | 计量专题 | 关于连享会

作者： Skylar B.Y. Sun

邮箱： [email protected]

1. Sequence 基本概念
2. sq.ado & sadi 命令安装
3. 数据导入
4. 自定义转换矩阵
5. Sequence data 基本数据生成
6. 计算 sequence distance
7. 生成不同聚类 (cluster)
8. Sequence 可视化

8.1 `sqparcoord`
8.2 `chronogram`
8.3 `sqindexplot`

9. 其他变量生成

编者按： 本文主要编译自以下两篇文章，特此感谢！

Source:

Halpin, B. (2017). SADI: Sequence analysis tools for Stata. The Stata Journal , 17 (3), 546-572. -PDF-

Brzinsky-Fay, C., Kohler, U., & Luniak, M. (2006). Sequence analysis with Stata. The Stata Journal , 6 (4), 435-460. -PDF-

温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。

生命历程调查 (Life History Analysis) 最近几年在社会学中越来越火热。与传统的只可以抓住特定时点的回归分析相比，生命历程调查最大的优势就是可以捕捉到一个人很长一段岁月中发生的故事。回归分析专注于找出样本的一个平均代表，而生命历程中的聚类分析则更倾向于通过细节找出一个群体中各种不同的小群体，旨在画出一幅包罗万象的丰富画卷。生命历程数据有很多研究方法，当下最火热的一种便是序列分析 (Sequence Analysis)。目前，常用的统计软件 R 和 Stata 都有序列分析相关的命令包，且 Stata 的命令相对更为丰富。因此，这篇就向大家介绍一下如何用 Stata 分析生命历程的数据。

1. Sequence 基本概念

Sequence 这个概念的缘起就是我们常说的 DNA 链条，后来逐渐发展成了社会学中的生命历程 sequence。Sequence 定义为遗传按次序排列的 元素 (element) ，比如说下图中的每个方格位置都是一个元素。每个元素都代表一个特定的 状态 (status) ，比如说工作状态、婚姻状态等。比如说下图中一共有五种状态，分别是 A B C E N。在一条 sequence 中，每个元素的都有其固定的位置。比如说 Position 7 上的元素是 A 。同一个元素连续出现了好几次便组成 元素片段 (episode / spell) 。比如说 P6 → P9 四个位置上的 A 就组成了一个元素片段。片段越长则代表对应的时间越长 (longer duration)。

Source: Cornwell, Benjamin. Social sequence analysis: Methods and applications . Vol. 37. Cambridge University Press, 2015. pp 59.

序列数据 (sequence data) 与截面时间序列数据 (cross-sectional time-series data) 及生存分析数据 (survival data) 有一些相似之处，比如说时间都是这种数据中的一个重要的组成部分。不过，与后两者不同的是，序列数据中的位置 (position) 一般指的是相对时间点，而不是一个绝对时间点。同时，每一条 sequence 都是被当成一个完整的个体去分析，且 sequence 中每一个元素都与这条 sequence 的特点息息相关。与生存分析数据还有一个不同，即序列数据中没有风险比的这个概念。

2. sq.ado & sadi 安装

Stata 目前有两套主要用于 sequence analysis 的命令。请按照如下方式安装。

net from http://teaching.sociology.ul.ie/sadi
net install sadi  /*导入最新版本的 sadi 命令*/

ssc install moremata
ssc install sq /*主要用来画 sequence index plot*/

3. 数据导入

在下面的例子中，我们将使用 McVicar 和 Anyadike-Danes (2002) 年的一套人物状态变换的数据。数据包括了连续72个月的数据点 / 位置 (state1 to state72)。每个位置有六种可能的元素，分别是受雇 (employment)、继续教育 (further education)、高等教育 (higher education)、高中教育 (secondary education)、培训 (training) 及失业 (unemployment)。

use http://teaching.sociology.ul.ie/bhalpin/mvad   /*数据导入*/

describe

Contains data from mvad.dta
  obs:           712                          
 vars:            86                          20 Dec 2020 20:53
 size:       244,928                          
----------------------------------------------------------
              storage   display    value
variable name   type    format     label      variable label
----------------------------------------------------------
id              float   %9.0g                 
state1          float   %9.0g      state      1 state
state2          float   %9.0g      state      2 state
state3          float   %9.0g      state      3 state
state4          float   %9.0g      state      4 state
state5          float   %9.0g      state      5 state
(output omitted)
state72         float   %9.0g      state      72 state
(output omitted)          
livboth         float   %9.0g                 
--------------------------------------------------------
Sorted by: id

4. 自定义转换矩阵

虽然之后我们也会介绍一种通过数据信息自己生成的转换矩阵，但在不少研究中，很多学者依旧倾向于自定义转换矩阵。

定义转换矩阵的原则主要是基于作者对于不同状态之间差别的理解。比如说，同一个状态之间的转换都是 0 (下图中的对角线都是 0) ，即从 Employment → Employment 之间并不需要转换。比如我们可能认为从 Employment → Unemployment 之间的差距是最大的，故定义这个值为 3。与此同时，我们可能觉得从 Employment → Further Education 或者 Employment → Higher Education 之间虽然有区别，但区别也不大，就定义这两个值为 1 。注意，一般来讲，我们会认为两种状态之间的转变，不论是哪个方向的转变 (比如说从 Employment → Unemployment 对比从 Unemployment → Employment ) 的差距都是一样的，故下面的矩阵通常是对称的。我目前读的文献还未见过用不同值代表两个方向的转变，如果大家有兴趣的话也可以自己试试。

	Employment	Further Education	Higher Education	Secondary Education	Training	Unemployment
Employment	0	1	1	2	1	3
Further Education	1	0	1	2	1	3
Higher Education	1	1	0	2	1	2
Secondary Education	2	2	2	0	1	1
Training	1	1	1	1	0	2
Unemployment	3	3	2	1	2	0

5. Sequence data 生成

以下命令用于生成基本的 sequence analysis 数据。

先来 reshape 一下数据。

*-先把数据转换成 long format
. reshape long state, i(id) j(order) 

(note: j = 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)

Data                               wide   ->   long   
-----------------------------------------------------
Number of obs.                      712   ->   51264
Number of variables                  86   ->      16
j variable (72 values)                    ->   order
xij variables:
              state1 state2 ... state72   ->   state
-----------------------------------------------------

对数据进行 sqset ，将其设定为序列分析的数据专有格式，其目的和作用与面板分析中的 xtset ，以及时间序列分析中的 tsset 命令相似。

. sqset st id order   

       element variable:  state, 1 to 6 /*6个state*/
       identifier variable:  id, 1 to 712
       order variable:  order, 1 to 72

看一看各种 sequence 样式的分布。

*-列表呈现
. sqtab, ranks (1/10)   /*列出前10个最常见的 sequence样式*/

Sequence-Patte |
            rn |      Freq.     Percent        Cum.
---------------+-----------------------------------
          1:72 |         40       39.22       39.22
     4:26 3:46 |         13       12.75       51.96
     5:24 1:48 |         11       10.78       62.75
     2:24 1:48 |          7        6.86       69.61
     4:27 3:45 |          6        5.88       75.49
 1:2 2:10 1:60 |          5        4.90       80.39
 1:2 2:22 1:48 |          5        4.90       85.29
     4:39 3:33 |          5        4.90       90.20
      6:2 1:70 |          5        4.90       95.10
 6:2 5:22 1:48 |          5        4.90      100.00
---------------+-----------------------------------
         Total |        102      100.00

*-描述性统计分析
. sqdes	 /*sequence concentration直观表示*/

   # of observed sequences: 712
overall # of obs. elements: 6
       max sequence length: 72
 # of producible sequences: 1.064e+56

----------------------------------------------------------
Observations |     Sequences  % of observed           Cum.
-------------+--------------------------------------------
           1 |           505       70.92696       70.92696  /*有505个sequencey样式只出现了一次*/
           2 |            28       3.932584       74.85955
           3 |             7       .9831461        75.8427
           4 |             7       .9831461       76.82584
           5 |             5       .7022472       77.52808
           6 |             1       .1404494       77.66853
           7 |             1       .1404494       77.80898
          11 |             1       .1404494       77.94943
          13 |             1       .1404494       78.08988
          40 |             1       .1404494       78.23033  /*有1种sequencey样式出现了40次*/
             | 
       Total |           557       78.23034               
----------------------------------------------------------

进一步，我们使用 egen 命令生成一些有用的变量。其他变量生成命令包括 sqelemcount() sqepicount() sqgapcount() 等等，可以通过 help sqegen 或者 help sqstat 找到更详细的说明。

/* () 中可以留空，这里是生成每条 sequence 的长度*/
. egen length = sqlength()  
 
/*不过我们这里的所有 sequence 长度都是一样的*/
. tab length   

  Length of |
   sequence |      Freq.     Percent        Cum.
------------+-----------------------------------
         72 |     51,264      100.00      100.00
------------+-----------------------------------
      Total |     51,264      100.00

/*生成每条 sequence 中含有状态1(employment)的长度*/
. egen length1 = sqlength(), element(1)  

/*id1的sequence里面，employment这个state出现了68次*/
. tab length1 if id == 1  

  Length of |
episodes of |
  element 1 |      Freq.     Percent        Cum.
------------+-----------------------------------
         68 |         72      100.00      100.00
------------+-----------------------------------
      Total |         72      100.00

6. 计算 sequence distance

每个 sequence 的组成都不同，计算两个 sequence 之间的不同的大小可以决定哪些 sequence 可以大致被分为一类。我这里讲两种距离计算的方法，一个用到了上面我们自定义的转换矩阵，另一个用的是用数据生成的矩阵。具体哪种方法更好呢，那就见仁见智啦。

先介绍第一种方法，用的就是我们自定义的转换矩阵。

. reshape wide state, i(id) j(order) /*先把数据 reshape 回 wide 格式*/

*-输入上述的自定义转换矩阵
#delimit ;
matrix mvdanes = (0,1,1,2,1,3 \
                  1,0,1,2,1,3 \
                  1,1,0,2,1,2 \
                  2,2,2,0,1,1 \
                  1,1,1,1,0,2 \
                  3,3,2,1,2,0 );
#delimit cr

set matsize 4000

/*用定义的 matrix 及 oma 算法，subsmat中就是我们定义的矩阵的名字*/
. oma state1-state72, subsmat(mvdanes) pwd(omd) length(72) indel(1.5)

下面用数据来生成转换矩阵，主要用到 trans2subs 这个命令。

/*用数据自己生成的转换矩阵来计算距离*/
. preserve
. reshape long state, i(id) j(m)
. trans2subs   state, id(id) subs(tpr1)
. matrix list tpr1

symmetric tpr1[6,6]
          c1        c2        c3        c4        c5        c6
r1         0
r2  1.147539         0
r3  1.064734  1.849958         0
r4  1.643575  1.757525  1.671111         0
r5  1.182927  1.844291      1.96   1.90181         0
r6  1.207729  1.525335  1.831594  1.803575  1.608297         0


. trans2subs state, id(id) subs(tpr2) diag

. matrix list tpr2

symmetric tpr2[6,6]
          c1        c2        c3        c4        c5        c6
r1         0
r2  1.967601         0
r3   1.98727  1.993341         0
r4  1.984684  1.987531  1.982969         0
r5  1.959993  1.992045  1.999488  1.994867         0
r6  1.951231   1.96336  1.996033  1.985649  1.972029         0

. restore

/*再次使用 oma 算法及生成的矩阵 tpr 来计算距离*/
oma        state1-state72, subsmat(tpr1) pwd(tpr) length(72) indel(1.5)

两个矩阵算出来的距离有什么区别吗？用 corrsqm 这个命令。

/*看看两种方法生成的 distance 是否相同*/
. corrsqm omd tpr

VECH correlation between omd and omv: 0.7726 /*correlation 还可以，不错！*/

7. 生成不同聚类

简单起见，我们就选择用 omd 这个距离来对所有人进行 sequence 分组。具体分成几组好呢，我们可以根据理论来看，也可以画出 index plot 来看。

分组生成命令 clustermat 及 cluster generate 如下：

clustermat wards omd, name(oma) add
cluster generate o=groups(8 12)  /*依次分成8组及12组；之后再看怎样分结果最好*/

tab	o8

	o8	Freq.	Percent	Cum.
				
	1	93	13.06	13.06
	2	139	19.52	32.58
	3	62	8.71	41.29
	4	146	20.51	61.80
	5	93	13.06	74.86
	6	30	4.21	79.07
	7	47	6.60	85.67
	8	102	14.33	100.00
				
	Total	712	100.00

tab o12

        o12 |      Freq.     Percent        Cum.
------------+-----------------------------------
          1 |         93       13.06       13.06
          2 |         42        5.90       18.96
          3 |         97       13.62       32.58  /*看来这里把上面的第2组劈开成了两组*/
          4 |         62        8.71       41.29
          5 |         89       12.50       53.79
          6 |         57        8.01       61.80  /*这里把上面的第4组劈开成了两组*/
          7 |         24        3.37       65.17
          8 |         35        4.92       70.08
          9 |         34        4.78       74.86  /*这里把上面的第5组劈开成了三组*/
         10 |         30        4.21       79.07
         11 |         47        6.60       85.67
         12 |        102       14.33      100.00
------------+-----------------------------------
      Total |        712      100.00

以分成八组为例的话，我们看看八组中每组的代表 sequence 长什么样子，即生成每组的 medoid sequence ，主要用到 discrepancy 这个命令。

stripe state1-state72, gen(seqstr) symbols("EFHSTU")  /*先生成seqstr变量，里面放入sequence的简化信息*/
list seqstr in 1/5, clean  /*列出前5个sequence的样子*/
                                                                         seqstr  
  1.   TTEEEETTEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE  
  2.   UUFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH  
  3.   UUTTTTTTTTTTTTTTTTTTTTTTTTFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEEEEEEEEEEUU  
  4.   TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTEEEEEEEEEEEEEEUUUUUUUUU  
  5.   UUFFFFFFFFFFFFFFFFFFFFFFFFFHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH  
  

discrepancy o8, dist(omd) id(id) dcg(dx) niter(1)
sort o8 dx
by o8: gen medoid = _n==1
list o8 dx seqstr if medoid, clean

       o8         dx                                                                     seqstr  
  1.    1   1.241993   EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE  
 94.    2   3.936442   TTTTTTTTTTTTTTTTTTTTTTTTEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE  
233.    3   3.302549   FFFFFFFFFFFFFFFFFFFFFFFFFFFHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH  
295.    4   6.059251   SSFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE  
441.    5   26.78911   TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTEEEEEEEEEEEEEEUUUUUUUUU  
534.    6   8.887777   TTTTTTTTTTTTTTTTTTTTUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUU  
564.    7   7.227705   SSSSSSSSSSSSSSSSSSSSSSSSSSEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE  
611.    8   3.599865   SSSSSSSSSSSSSSSSSSSSSSSSSSHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH

可以看出，如果分成八组的话，第一组的代表 sequence jiushi从头到尾一直是处于 employment 状态的人。

8. Sequence 可视化

Sequence 图像主要有三大类， sqparcoord ， chronogram ， sqindexplot ，及 transition pattern graph。 sqparcoord 可以画出出现频率比较高的 sequence。Chronogram 主要为了看出在每种 state 里面的持续时间。Indexplot 把每条 sequence 都化成了一条横线，然后叠加在一起看规律。Transition pattern graph 主要可以看出整体而言，每种 state 之间的转化情况。

8.1 `sqparcoord`

先来看看几条出现频率最高的 sequence 的样子吧！

preserve
reshape long state, i(id) j(order)  /*sqparcoord只可以在 long format下使用，因此 reshape to long*/
sqparcoord, so offset(0.5) wlines(1) ylabel(1(1)6)  
restore

可以看出呢，不少人一直处于 employment 状态（从头至尾都是状态 1）。也可以看出不少人从状态5 (Training) 转换到状态1。我随便改了改一些线条的颜色，因为我觉得彩色一点比较好看。

8.2 `chronogram`

chronogram state*, by(o8, legend(off)) name(chronogram, replace)  /*用之前的分成的八组来生成 chronogram*/

我把 legend(off) 了……不过还是可以看出，第一组的话，蓝色代表 employment。