第2章 序列的数据挖掘
数据科学：R语言实战
数据挖掘技术一般用于检测数据中的序列或模式。本章中，我们将试图使数据遵循一种模式，在此模式中，一个或一系列事件可以通过一致的方式预测另一个数据点。
本章描述了在数据集中查找模式的不同方法：
- 查找模式
- 在数据中查找模式
- 约束条件
我们可以在很多较大的数据集中查找模式。这涵盖了很多区域，比如混合人口的变化、使用手机的频率、高速公路质量衰退、年龄因素造成的事故等。不过我们能明确地感受到，有很多模式和序列正等待我们去发现。
我们可以通过使用R编程中的一些工具找到这些模式。大多数模式因约束条件而在一定程度上受到限制，如序列的有用时间。
2.1 模式
我们来回顾一下确定数据中模式的方法：
模型类别 模型如何工作
2.1.1 Eclat
Eclat算法用于频繁项集的挖掘。这种情况下，我们寻找行为相似的模式，与之相对的是寻找不规则模式（与处理其他数据挖掘的方法类似）。
Algorithm通过数据中的交集来估算同时频繁出现事件候选项（如购物车项目）的支持度。然后通过对频繁候选项进行测试来证实数据集中的模式。
1．用法
在R编程中使用Eclat就是使用arules功能包中的eclat函数。使用Eclat算法的R编程遵循了此处提出的约定：
> eclat(data,
parameter = NULL,
control = NULL)
下表对eclat函数的不同参数进行了说明：
常见的ECParameter如下所示：
ECControl常见值如下所示：
调用eclat函数返回数据中出现的频繁项集。
Eclat的实施包括成年人数据集。成年人数据集包括人口统计局数据中约50000行的数据。
2．使用eclat找到成年人行为的相似点
使用下列代码找到成年人行为的相似点：
> library("arules")
> data("Adult")
> dim(Adult)
[1] 48842 115
> summary(Adult)
transactions as itemMatrix in sparse format with
48842 rows (elements/itemsets/transactions) and
115 columns (items) and a density of 0.1089939
most frequent items:
capital-loss=None capital-gain=None
46560 44807
native-country=United-States race=White
43832 41762
workclass=Private (Other)
33906 401333
element (itemset/transaction) length distribution:
sizes
9 10 11 12 13
19 971 2067 15623 30162
Min. 1st Qu. Median Mean 3rd Qu. Max
9.00 12.00 13.00 12.53 13.00 13.00
includes extended item information - examples:
labels variables levels
1 age=Young age Young
2 age=Middle-aged age Middle-aged
3 age=Senior age Senior
includes extended transaction information - examples:
transactionID
1 1
2 2
3 3
检查最终结果时，我们会注意到以下细节：
摘要共48842行，115列。
已列出常见项目：白种人。
有很多描述符，如age=Young。
3．查找数据集中的频繁项目
处理数据集时，通过下列代码挖掘出现的频繁项集：
> data("Adult")
> itemsets <- eclat(Adult)
parameter specification:
tidLists support minlenmaxlen target ext
FALSE 0.1 1 10 frequent itemsets FALSE
algorithmic control:
sparse sort verbose
7 -2 TRUE
eclat - find frequent item sets with the eclat algorithm
version 2.6 (2004.08.16) (c) 2002-2004 Christian Borgelt
createitemset ...
set transactions ...[115 item(s), 48842 transaction(s)] done [0.03s].
sorting and recoding items ... [31 item(s)] done [0.00s].
creating bit matrix ... [31 row(s), 48842 column(s)] done [0.02s].
writing ... [2616 set(s)] done [0.00s].
Creating S4 object ... done [0.00s].
默认值已发现2616个频繁集合。如果我们寻找前五个集合，将会看到下列输出数据：
> itemsets.sorted <- sort(itemsets)
> itemsets.sorted[1:5]
items support
1 {capital-loss=None} 0.9532779
2 {capital-gain=None} 0.9173867
3 {native-country=United-States} 0.8974243
4 {capital-gain=None,
capital-loss=None} 0.8706646
5 {race=White} 0.8550428
以下是对之前输出数据的研究所得：
- 普查数据中的大多数人未要求资本损失或资本利得（这种财政税收事件并非正常状态）。
- 大多数人来自美国。
- 大多数是白种人。
4．集中于最高频率的示例
为了进一步证实数据，我们可以将范围缩减至数据集中出现的最高频率（可以通过调节minlen参数直至处理完一项集合来实现操作）：
> itemsets <- eclat(Adult, parameter=list(minlen=9))
> inspect(itemsets)
items support
1 {age=Middle-aged,
workclass=Private,
marital-status=Married-civ-spouse,
relationship=Husband,
race=White,
sex=Male,
capital-gain=None,
capital-loss=None,
native-country=United-States} 0.1056673
按照预期，由一位美国本土且拥有工作的已婚男士填写普查数据表格。
2.1.2 arulesNBMiner
在R中，arulesNBMiner是一个功能包，用于寻找一个集合中两个或两个以上项目的共现。底层模型，即负二项式模型，允许高度偏态次数分配，否则会很难确定最小项集容量。我们在正被挖掘的较大数据集中寻找频繁数据集。当确定使用arulesNBMiner时，您应该看到一些迹象：项目集频率正出现在数据子集合中。
1．用法
将arulesNBMiner作为功能包进行操作，并且必须将此功能包安装于您的R编程环境中。可以通过使用任意数据集来学习如何使用模型/函数中包含的工具，如下所示：
> results <-NBMiner(data, parameter, control = NULL)
下表对NBMiner函数的不同参数进行了说明：
NBMinerParameters的数值如下所示：
NBMinerParameters(data, trim = 0.01, pi = 0.99,
theta = 0.5, minlen = 1, maxlen = 5, rules = FALSE,
plot = FALSE, verbose = FALSE, getdata = FALSE)
功能包中的Agrawal数据可以直接使用。注意：Agrawal数据是为集中事务通过特别合成而生成的。代码如下所示：
> data(Agrawal)
> summary(Agrawal.db)
transactions as itemMatrix in sparse format with
20000 rows (elements/itemsets/transactions) and
1000 columns (items) and a density of 0.00997795
most frequent items:
item540 item155 item803 item741 item399 (Other)
1848 1477 1332 1295 1264 192343
element (itemset/transaction) length distribution:
sizes
1 2 3 4 5 6 7 8 9 10 11 12 13
15 88 204 413 737 1233 1802 2217 2452 2444 2304 1858 1492
14 15 16 17 18 19 20 21 22 23 24 25
1072 706 431 233 138 83 46 19 10 1 1 1
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.000 8.000 10.000 9.978 12.000 25.000
includes extended item information - examples:
labels
1 item1
2 item2
3 item3
includes extended transaction information - examples:
transactionID
1 trans1
2 trans2
3 trans3
> summary(Agrawal.pat)
set of 2000 itemsets
most frequent items:
item399 item475 item756 item594 item293 (Other)
29 29 29 28 26 3960
element (itemset/transaction) length distribution:sizes
1 2 3 4 5 6
702 733 385 134 34 12
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.00 1.00 2.00 2.05 3.00 6.00
summary of quality measures:
pWeightspCorrupts
Min.:2.100e-08 Min. :0.0000
1st Qu.:1.426e-04 1st Qu. :0.2885
Median :3.431e-04 Median :0.5129
Mean :5.000e-04 Mean :0.5061
3rd Qu.:6.861e-04 3rd Qu. :0.7232
Max. :3.898e-03 Max. :1.0000
includes transaction ID lists: FALSE
以下是对之前输出数据的研究所得：
共20000行，1000列。
所有纵列都用项目399、项目475等命名。
2000个子集在少数例子中具有偏态（如容量1有702、容量2有733等）。
2．为频繁集挖掘Agrawal数据
如果以Agrawal数据为例，可以得到下列输出数据：
> mynbparameters <- NBMinerParameters(Agrawal.db)
> mynbminer <- NBMiner(Agrawal.db, parameter = mynbparameters)
> summary(mynbminer)
set of 3332 itemsets
most frequent items:
item540 item615 item258 item594 item293 (Other)
69 57 55 50 46 6813
element (itemset/transaction) length distribution:sizes
1 2 3 4 5
1000 1287 725 259 61
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.000 1.000 2.000 2.128 3.000 5.000
summary of quality measures:
precision
Min.:0.9901
1st Qu. :1.0000
Median :1.0000
Mean :0.9997
3rd Qu. :1.0000
Max. :1.0000
以下是对之前输出数据的研究所得：
项目近乎均匀分布。
项集长度1或2有较大偏斜。
2.1.3 Apriori
Apriori是可以帮助了解关联规则的分类算法。与事务的实施方式相对。这种算法尝试找到数据集中常见的子集合，必须满足最小阈值以便核实关联。
Apriori的支持度和置信度概念十分有趣。Apriori方法会从您的数据集中返回有趣的关联，如当出现Y时，会返回X。支持度是包含X和Y的事务的百分比。置信度是同时包含X和Y的事务的百分比。支持度的默认值为10，置信度的默认值为80。
1．用法
apriori方法的使用规则如下所示：
apriori(data, parameter = NULL, appearance = NULL, control = NULL)
下表对apriori函数的不同参数进行了说明：
2．评估购物篮中的关联
我们正寻找食品超市中典型购物篮内购买的项目之间的关联。为此，我们将按下列步骤进行操作。
（1）下载下列arules功能包：
> install.packages("arules")
> library(arules)
（2）下载事务，即比利时杂货店数据：
> tr <- read.transactions("http://fimi.ua.ac.be/data/retail.dat",
format="basket")
（3）大概了解数据：
>summary(tr)
transactions as itemMatrix in sparse format with
88162 rows (elements/itemsets/transactions) and
16470 colunsis (items) and a density of 0.0006257289
most frequent items:
39 48 38 32 41 (Other)
50675 42135 15596 15167 14945 770058
element (itemset/transaction) length distribution:
sizes
1 2 3 4 5 6 7 8 9 10 11 12 13
3016 5516 6919 7210 6814 6163 5746 5143 4660 4086 3751 3285 2866
14 15 16 17 18 19 20 21 22 23 24 25 26
2620 2310 2115 1874 1645 1469 1290 1205 981 887 819 684 586
27 28 29 30 31 32 33 34 35 36 37 38 39
582 472 480 355 310 303 272 234 194 136 153 123 115
40 41 42 43 44 45 46 47 48 49 50 51 52
112 76 66 71 60 50 44 37 37 33 22 24 21
53 54 55 56 57 58 59 60 61 62 63 64 65
21 10 11 10 9 11 4 9 7 4 5 2 2
66 67 68 71 73 74 76
5 3 3 1 1 1 1
Min. 1st Qu. Median Mean 3rd QU. Max.
1.00 4.00 8.00 10.31 14.00 76.00
includes extended item information - examples:
labels
1 0
2 1
3 10
以下是对之前输出数据的研究所得：
共88162个购物篮，对应16470个项目。
成对项目很受欢迎（项目39有50675个）。
（4）一起看一下最频繁的项目：
> itemFrequencyPlot(tr, support=0.1)
我们可以再次看到少数频率比平常频率更高的项目。
（5）现在，为合适的关联构建一些规则：
> rules <- apriori(tr, parameter=list(supp=0.5,conf=0.5))
parameter specification:
confidenceminvalsmaxaremavaloriginalSupport support minlen
0.5 0.1 1 none FALSE TRUE 0.5 1
maxlen target ext
10 rules FALSE
algorithmic control:
filter tree heap memopt load sort verbose
0.1 TRUE TRUE FALSE TRUE 2 TRUE
apriori - find association rules with the apriori algorithm
version 4.21 (2004.05.09) (c) 1996-2004 Christian Borgelt
set item appearances ...[0 item(s)] done [0.00s].
set transactions ...[16470 item(s), 88162 transaction(s)] done
[0.13s].
sorting and recoding items ... [1 item(s)] done [0.01s].
creating transaction tree ... done [0.02s].
checking subsets of size 1 done [0.00s].
writing ... [1 rule(s)] done [0.00s].
creating S4 object ... done [0.01s].
（6）然后用一条规则作为结束。规则摘要如下：
> summary(rules)
set of 1 rules
rule length distribution (lhs + rhs):sizes
1
1
Min. 1st Qu. Median Mean 3rd Qu. Max.
1 1 1 1 1 1
summary of quality measures:
support confidence lift
Min.:0.5748 Min. : 0.5748 Min. :1
1st Qu. :0.5748 1st Qu. :0.5748 1st Qu. :1
Median :0.5748 Median :0.5748 Median :1
Mean :0.5748 Mean :0.5748 Mean :1
3rd Qu. :0.5748 3rd Qu. :0.5748 3rd Qu. :l
Max. :0.5748 Max. :0.5748 Max. :1
mining info:
datantransactions support confidence
tr 88162 0.5 0.5
规则的支持度有力，置信度较低。
（7）具体规则：
>inspect(rules)
lhsrhs support confidence lift
1{} => {39} 0.5747941 0.5747941 1
正如我们猜想的那样，大多数人将项目39放入购物篮。
（8）我们可以寻找更多与规则相关的信息，以便全面了解其产生的影响。
>interestMeasure(rules, c("support", "chiSquare", "confidence",
"conviction", "cosine", "leverage", "lift", "oddsRatio"), tr)
sapply(method, FUN = function(m) interestMeasure(x, m,
transactions, reuse, ...))
support 0.5747941
chiSquareNaN
confidence 0.5747941
conviction 1.0000000
cosine 0.7581518
leverage 0.0000000
lift 1.0000000
oddsRatioNaN
一条派生规则将这些措施的置信度展现得非常全面。
2.1.4 用TraMineR确定序列
TraMineR功能包用于挖掘序列，并将其可视化，其思想是发现序列。可以将序列分布、序列频率及湍流等绘图的图解设备构建到功能包中。此外，还有一些自然出现的项目，其中的数据有重复的序列，如在一些社会科学场地，数据会自然地循环项目。
通过此文件，我将带您大概了解TraMineR，以便生成一系列用于发现序列的工具。在挖掘操作中选择何种工具取决于您自己。
您可以同时使用TraMineR功能包及一对内置数据集：
1．用法
seqdef函数用于确定数据中出现的序列：
seqdef(data, var=NULL, informat="STS", stsep=NULL,
alphabet=NULL, states=NULL, id=NULL, weights=NULL,
start=1, left=NA, right="DEL", gaps=NA,
missing=NA, void="%", nr="*", cnames=NULL,
xtstep=1, cpal=NULL, missing.color="darkgrey",
labels=NULL, ...)
下表对seqdef函数的不同参数进行了说明：
2．确定训练和职业中的序列
在这一示例中，我们将看到人们生活中从训练到工作的进程中时间的序列。我们期望看到从失业未经训练的状态至经过训练并最终成为全职员工的进程。
TraMineR功能包中的一些函数有助于序列分析。我们使用seqdef来创建数据对象，以便其他函数也可以使用。可以用其他方法设置或保留参数，如下所示：
seqdef(data, var=NULL, informat="STS", stsep=NULL,
alphabet=NULL, states=NULL, id=NULL, weights=NULL, start=1,
left=NA, right="DEL", gaps=NA, missing=NA, void="%", nr="*",
cnames=NULL, xtstep=1, cpal=NULL, missing.color="darkgrey",
labels=NULL, ...)
大多数参数可以在默认值下使用。
如您所见，seq数据对象是plot函数的第一个参数。您可以用实际想用的plot函数（如下面编码中所用的seqiplot）代替XXX。
seqXXXplot(seqdata, group=NULL, type="i", title=NULL,
cpal=NULL, missing.color=NULL,
ylab=NULL, yaxis=TRUE, axes="all", xtlab=NULL, cex.plot=1,
withlegend="auto", ltext=NULL, cex.legend=1,
use.layout=(!is.null(group) | withlegend!=FALSE),
legend.prop=NA, rows=NA, cols=NA, ...)
多数参数是您在plot中所需的标准图像的强化形式，如ylab是y轴的标记。
首先，我们必须用下列编码将TraMineR加载入您的环境中。
> install.packages("TraMineR")
> library ("TraMineR")
我们将使用TraMineR功能包中内置的mvad数据集。mvad数据集追踪了712个个体在20世纪90年代自训练至工作的进程。我们可以按下列形式使用mvad数据集。
> data(mvad)
数据摘要如下所示：
> summary(mvad)
id weight male catholic Belfast
Min.: 1.0 Min.:0.1300 no :342 no :368 no :624
1st Qu. :178.8 1st Qu. :0.4500 yes:370 yes:344 yes: 88
Median :356.5 Median :0.6900
Mean :356.5 Mean :0.9994
3rd Qu. :534.2 3rd Qu. :1.0700
Max. :712.0 Max. :4.4600
N.EasternSouthern S.Eastern Western Grammar funemp
no :503 no :497 no :629 no :595 no :583 no :595
yes:209 yes:215 yes: 83 yes:117 yes:129 yes:117
gcse5eqfmprlivboth Jul.93
no :452 no :537 no :261 school :135
yes:260 yes:175 yes:451 FE : 97
employment :173
training :122
joblessness:185
HE : 0
我们可以查看标准标识符来了解体重、性别、宗教等信息。
截取序列数据（我们正通过86使用17列，因为这适用于人们在数据调查不同点的状态），并将数据的这部分应用于序列确定函数，如下所示：
> myseq <- seqdef(mvad, 17:86)
[>] 6 distinct states appear in the data:
1 = employment
2 = FE
3 = HE
4 = joblessness
5 = school
6 = training
[>] state coding:
[alphabet] [label] [long label]
1 employmentemploymentemployment
2 FEFEFE
3 HEHEHE
4 joblessnessjoblessnessjoblessness
5 schoolschoolschool
6 trainingtrainingtraining
[>] 712 sequences in the data set
[>] min/max sequence length: 70/70
这样看来是正确的，我们可以参照相关状态（失业、上学、训练及工作）来获取所需的行序列数据。
我们可以使用一些内置图表来将已确定的序列可视化。序列如下所示。
- seqiplot：指数图表。
- seqfplot：频率图表。
- seqdplot：分布图表。
以指数图表为例：
> seqiplot(myseq)
通过参照个人不同状态间界定的转换期，您会发现连续几个月都有训练。您应进行核实，以便数据显示的信息与您对序列数据的理解相一致。
现在，以频率图表为例：
> seqfplot(myseq)
现在我们来看序列在不同时间的频率。多次观看后我们会看到同一序列的人群集，如经过一段时间的训练后会有工作。
现在，以分布图表为例：
> seqdplot(myseq)
我们来看看序列状态在不同时期的分布情况。通常情况下，人们在上学或训练后开始工作。明白了吧！
通过以下指令我们可以看到序列的熵：
> seqHtplot(myseq)
熵在不同时期的变化特点：明显降低后会出现细微的上升。这与不同人群会在最初做出不同选择的情况一致（很多状态），如上学或训练，然后进行工作，成为劳动力（一种状态）。
有一个有趣的想法为数据湍流。湍流传达出一个信息，即从数据中可见的某个特定事例可以推导出多少不同的后续序列。我们可以用seqST函数将湍流可视化。seqST函数将序列数据作为参数，并返还湍流数据。让我们继续以示例进行说明：
> myturbulence <- seqSt(myseq)
> hist(myturbulence)
我们可以看到带有长尾数的近乎标准化分布。大多数状态分为少量后续状态以及少数状态或多或少的异常值。
2.1.5 序列相似点
TraMineR功能包也可以确定序列的度量，如不同序列间的相异点。
- 最长公共前缀（LCP）：我们可以通过比较相同的最长序列前缀来确定相似点。
- 最长公共序列（LCS）：我们也可以通过查看两个序列之间的相同部分，根据其内部的最长序列来确定相似点。
- 最佳匹配（OM）距离：指生成一个不同序列的最佳编辑距离，在此距离下，插入及删除的成本最小。
使用TraMineR中的seqdist函数可以实现所有这些功能。
1．序列度量
我们可以用seqdist计算LCP。
2．用法
seqdist函数使用规则如下：
seqdist(seqdata, method, refseq=NULL, norm=FALSE,
indel=1, sm=NA, with.missing=FALSE, full.matrix=TRUE)
下表对seqdist函数的不同参数进行了说明：
3．示例
seqdist函数用法示例如下：
（1）使用被构建到功能包的famform序列：
> data(famform)
（2）界定可用的序列对象：
> seq <- seqdef(famform)
[>] found missing values ('NA') in sequence data
[>] preparing 5 sequences
[>] coding void elements with '%' and missing values with '*'
[>] 5 distinct states appear in the data:
1 = M
2 = MC
3 = S
4 = SC
5 = U
[>] state coding:
[alphabet] [label] [long label]
1 MMM
2 MCMCMC
3 SSS
4 SCSCSC
5 UUU
[>] 5 sequences in the data set
[>] min/max sequence length: 2/5
> seq
Sequence
[1] S-U
[2] S-U-M
[3] S-U-M-MC
[4] S-U-M-MC-SC
[5] U-M-MC
（3）确定使用序列3和序列4的LCP：
> seqLLCP(seq[3,],seq[4,])
[1] 4
我们可以得到四个前置匹配（S-U-M-MC与S-U-M-MC-SC相比）。
（4）我们可以直接计算LCS度量：
> seqLLCS(seq[1,],seq[2,])
[1] 2
我们可以在2中找到常用序列。
（5）也可以直接确定OMD：
[>] cost <- seqsubm(seq, method-"CONSTANT". Cval=2)
[>] creating 5x5 substitution-cost matrix using 2 as constant
value
> cost
M-> MC-> S-> SC-> U->
M-> 0 2 2 2 2
MC-> 2 0 2 2 2
S-> 2 2 0 2 2
SC-> 2 2 2 0 2
U-> 2 2 2 2 0
OMD仅为2（仅用较小序列来诠释概念）。