书名:数据科学:R语言实战
ISBN:978-7-115-43590-3
本书由人民邮电出版社发行数字版。版权所有,侵权必究。
您购买的人民邮电出版社电子书仅供您个人使用,未经授权,不得以任何方式复制和传播本书内容。
我们愿意相信读者具有这样的良知和觉悟,与我们共同保护知识产权。
如果购买者有侵权行为,我们可能对该用户实施包括但不限于关闭该帐号等维权措施,并可能追究法律责任。
• 著 [美] Dan Toomey
译 刘丽君 李成华 卢青峰
责任编辑 王峰松
• 人民邮电出版社出版发行 北京市丰台区成寿寺路11号
邮编 100164 电子邮件 315@ptpress.com.cn
• 读者服务热线:(010)81055410
反盗版热线:(010)81055315
Copyright © Packt Publishing 2014. First published in the English language under the title R for Data Science,ISBN 978-1-78439-086-0. All rights reserved.
本书中文简体字版由Packt Publishing公司授权人民邮电出版社出版。未经出版者书面许可,对本书的任何部分不得以任何方式或任何手段复制和传播。
版权所有,侵权必究。
本书讲述的是R语言在数据科学中的应用,目标读者是从事不同行业的数据分析师、数据挖掘工程师、机器学习工程师、自然语言处理工程师、数据科学家,以及从事大数据和人工智能领域的工作者、学生、老师等。
本书最大的优点在于其通俗易懂、容易上手,每一个实例都有现成的数据和源代码,读者不仅能理解整个案例的来龙去脉,还可以直接编译本书提供的所有源代码,从而了解怎么从实际问题转变成可实现的代码,感受R语言的魅力,让数据产生价值。这种学习和实践相结合的方式非常适合初学者和有一定经验的数据分析师。
本书的内容涵盖了基于数据挖掘的常用模型,包括分类、聚类、关联分析、预测、异常检测等,还包括机器学习的常用算法和自然语言处理、数据可视化等内容。本书内容全面,做到了易读、易用、易理解、易实现、易上手,是不可多得的R语言书籍。
Dan Toomey具有20多年开发应用程序方面的经验,曾在多个行业及公司担任不同的职位,包括唯一投稿人、副总裁及首席技术官。近10年,Dan一直在美国马萨诸塞州东部地区的公司工作。Dan以Dan Toomey软件公司的名义,成为这些领域的开发承包商。
刘丽君,韩国国立全北大学博士,加拿大圣西维尔大学博士后,一直从事物联网、工业大数据等方面的数据分析、市场分析等工作,目前任武汉泰迪智慧科技有限公司CEO,对数据敏感,并对数据怎么转变成价值、数据与商业的关系有独到见解。
李成华,数据挖掘与机器学习方向博士,约克大学博士后,麻省理工学院访问科学家,曾任海信集团数据挖掘专家,京东深度神经网络实验室首席科学家,长期从事数据挖掘、机器学习、深度学习和自然语言处理等方面的研究和工作,擅长自动问答以及基于自然语言的人机交互。
卢青峰,硕士毕业于美国威斯康辛州立大学,毕业后从事数据分析、挖掘等相关工作至今。曾先后在敦煌网、百度、京东等行业领先的公司从事数据挖掘、用户行为分析、推荐等工作。
在本书翻译的过程中,来自“统计之都”的核心成员,包括九峰医疗的数据研究总监李舰、雪晴数据网的创始人陈堰平、京东搜索与大数据平台部的刘思喆等,对本书的中文译稿进行了仔细的审阅,并提出了很多有价值的意见和建议。在此对他们的付出和努力表示感谢。
Amar Gondaliya是数据科学家,任职于一家著名的医疗机构。他是大数据及数据挖掘的爱好者,也是Pingax(www.pingax.com)的顾问,专注于使用数据挖掘技术构建预测模型。他热衷于研究大数据技术,并提出大数据技术的解决方案。
他在数据科学领域工作了2年多,是Pingax的投稿人,撰写关于机器学习及其在R中实现的文章。他一直致力于机器学习新技术及大数据分析的研究。
Mohammad Rafi是一名软件工程师,热衷数据分析、编程以及修修补补。他致力于R、Python、Hadoop及JavaScript等技术的研究。白天,他是工程师;晚上,他是游戏发烧友。自撰写本书起,他热衷于Raspberry Pi超级小电脑。
他在应用开发、数据处理、搜索专家及网络数据分析方面有6年以上的专业经验。最初,他在一家名为Position2的网络营销公司工作。此后,先后在《印度斯坦时报》、谷歌及InMobi等公司工作。
Tengfei Yin获得了中国南开大学生命科学与技术学院的理学学士学位,并且在美国爱荷华州立大学完成了分子、细胞和发育生物学(MCDB)的学习,重点研究计算生物学和生物信息学。他的研究领域包括信息可视化、高通量生物数据分析、数据挖掘、机器学习及应用统计遗传学。他开发并且维护了R及Bioconductor中的多种软件包。
R是为数据操作及统计计算提供语言及环境的软件包,同样也能够用图表表示产生的统计数据。
R具有以下特性:
第1章:模式的数据挖掘,包括R中的数据挖掘。在此示例中,我们会寻找数据集中的模式。本章探讨通过使用多种工具进行聚类分析的示例,同样也包括了异常检测和关联规则的使用。
第2章:序列的数据挖掘,探讨了R中可以让读者发现数据中序列的方法。多种可用的R功能包能够帮助读者确定序列并通过图形描绘做进一步分析。
第3章:文本挖掘,介绍了多种挖掘R中文本的方法。我们会着眼于能够帮读者处理并分析源中文本或文字的工具,同样也会研究XML处理能力。
第4章:数据分析——回归分析,探讨了对数据进行回归分析的不同方法。本章运用多种方法来运行简易回归和多元回归以及后续的显示。
第5章:数据分析——相关性,探讨了多种相关功能包。本章不仅运用皮尔森、多分格、四分、异构及部分相关性,还运用基本的相关性及协方差对数据进行分析。
第6章:数据分析——聚类,探讨了各种聚类分析的参考文献。本章包括K-means、PAM和若干其他聚类技术。R程序员可利用所有技术。
第7章:数据可视化——R图形,探讨了各种使数据可视化的方法。我们会着眼于数据的全部范围,包括典型的类显示、第三方工具的相互作用和地图的使用。
第8章:数据可视化——绘图,探讨了绘制R中数据的不同方法。本章不仅提供了可用于绘制数据的自定义显示的示例,还提供了带有标准化显示的简易绘图的示例。
第9章:数据可视化——三维,作为直接从R创建数据三维显示的指南,我们同样也会考虑使用更大数据集的三维显示。
第10章:机器学习实战,探讨了如何使用R进行机器学习。本章包括把数据集分成训练数据及测试数据,从训练数据中构建模型并对模型进行试验,与试验数据进行比较。
第11章:用机器学习预测事件,使用了时间序列数据集。本章包括将数据转化成R时间序列,然后将季节性、趋势及不规则分量分离出来。其目标在于模拟或预测未来事件。
第12章:监督学习和无监督学习,阐释了如何使用监督和无监督学习来构建模型。本章包括多种有关监督和无监督学习的方法。
关于本书,读者需要将R安装在机器(或将要运行脚本的机器)上。许多平台都可使用R。本书并非局限于R的现有特定版本。
为了更好地使用本书,读者需要开发R项目的交互工具。主要工具是R Studio,是许多平台均可使用的完全交互式且自包含程序,可以让读者输入脚本、显示数据并显示图形结果。各种R的安装总会有R命令行工具。
本书的目标读者是牢固掌握先进数据分析技术的数据分析师。同样也要求基本了解R语言及某些数据科学话题。本书假设读者可以使用R环境,并且对涉及的统计有所了解。
欢迎来自读者的反馈,以让我们了解读者对本书的看法——读者喜欢或不喜欢的内容。读者反馈是重要的部分,能够帮助我们更好地阐述本书的主题,以便读者以后能够更好地掌握其内容。
将一般反馈发送至邮箱:feedback@packtpub.com,需在邮件主题中标明本书的标题。
如果读者擅长某一话题,并且对写书或撰稿感兴趣,登录www.packtpub.com/authors,可查看作者指南。
尽管我们已格外小心以确保内容的准确性,仍然无法避免书中会出现疏漏。如果读者发现在我们出版的书籍中存在错误——可能是文字或代码错误——请告知我们,我们将感激不尽。疏漏改正后,可以避免影响其他读者的阅读体验,并且有助于提高本书的后续版本。如果读者发现任何勘误,请按照以下方法告知:登录http://www.packtpub.com/ submit-errata,选择书籍,单击“勘误提交表格”链接,并输入勘误内容。一旦勘误通过核实,将会接收提交,并且勘误会上传至我们的网站或添加至此标题下“勘误”一节中现有的勘误表。
登录https://www.packtpub.com/books/content/support并在“搜索栏”输入书名,可查看以前提交的勘误。“勘误”一节中会提供必要信息。
非法翻印网上有版权的材料是所有媒体中存在已久的问题。在Packt,我们很重视对自身版权及授权的保护。如果读者在网上看见以任何形式盗版我们著作的行为,请立即向我们提供位置地址或网站名称,以便我们采取补救措施。
请将带有涉嫌盗版内容链接的邮件发送至copyright@packtpub.com与我们联系。
感谢读者对保护作者及我们能为读者提供有价值的内容而提供的帮助。
如果读者对本书留有疑问,可发送邮件至questions@packtpub.com与我们联系,我们会竭尽所能解答读者的疑问。
数据挖掘常用于检测数据中的模式或规则。
兴趣点在于仅能够通过使用大数据集进行检测的不明显模式。一段时间内可以检测更简易的模式,如用于购买关联或时间选择的购物篮分析。我们对R编程的兴趣在于检测意外的关联,这能够带来新的机会。
某些模式本质上是有序的,例如,基于以往结果预测系统中的故障,通过使用大数据集,以往结果会更加明确。下一章会探讨相关内容。
本章探讨使用R来发现数据集不同方法中的模式。
通过运用各种算法(如下表列举的一些算法)可进行聚类分析。
模型类别 | 模型如何工作 |
---|---|
连通性 | 此模型计算了点间的距离,并且基于距离远近对点进行组织 |
分割 | 此模型将数据分割成集群,并且将每个数据点与一个集群结合起来,最主要的集群是K-means |
分布模型 | 此模型使用了统计分布以便确定集群 |
密度 | 此模型确定了数据点的紧密度以便到达分布的密集区。DBSCAN常用于密集的分布,OPTICS常用于更加稀疏的分布 |
在算法中,同样也有更高级别的粒度,其中包括硬聚类或软聚类等。
在R编程中,聚类工具用于:
K-means聚类是将数据集分割成k集群的方法。需要预先确定将数据集分成几个集群。K-means算法的步骤如下所示。
(1)从数据中选取k随机行(质心)(可使用预先确定的集群数量)。
(2)我们使用Lloyd's算法(系统默认)来确定集群。
(3)根据与质心的距离而对每个数据点进行分配。
(4)将每个质心重新计算为与其相关的所有点的平均值。
(5)对与质心距离最近的数据点进行重新分配。
(6)继续执行步骤3和步骤4,直到不再分配数据点或循环次数已为最大值为止。
这是启发式算法,最好多次运行流程。通常,此算法会在R中快速运行,原因在于每个步骤都不难执行。其目的是通过不断完善术语,将平方和降到最小。
预先确定集群数量可能会出现问题。用图表表示数据(或平方等数值)应直观地呈现数据的逻辑分组。通过确定截止进行选择的步骤来进行迭代,以此来确定分组大小(本章后面会使用此步骤)。同样也有其他企图计算选择截止的R功能包。完成后,应核实选取集群的拟合。
使用平均值(步骤3中)表示K-means无法与相当稀疏的数据或带有许多异常值的数据一起工作。此外,如果集群的线性形状不好,就会出现问题。图示应证明数据是否适合此算法。
用kmeans函数在R编程中进行K-means聚类。使用k-means聚类的R编程遵循了此处提出的约定(注意:读者可以通过使用内联help函数确定函数的约定以便获取此信息,如?kmeans):
kmeans(x,
centers,
iter.max = 10,
nstart = 1,
algorithm = c("Hartigan-Wong",
"Lloyd",
"Forgy",
"MacQueen"),
trace=FALSE)
下表对不同参数进行了说明。
参数 | 描述 |
---|---|
x | 这是待分析的数据矩阵 |
centers | 这是集群数量 |
iter.max | 这是最大迭代次数(除非停止重赋值) |
nstart | 这是随机集的使用次数 |
algorithm | 可以是下列类型之一:Hartigan-Wong、Lloyd、Forgy或MacQueen算法 |
trace | 随着算法的不断发展,这提供了当前的跟踪信息 |
调用kmeans函数能够返回带有下列属性的kmeans对象。
属性 | 描述 |
---|---|
cluster | 包含了集群分配 |
centers | 包含了集群中心 |
totss | 提供了总平方和 |
withinss | 这是每个聚类内部平方和的向量 |
tot.withinss | 包含了距离平方和总量 |
betweenss | 包含了聚类组间的平方和 |
size | 包含了每个聚类内的数据点数量 |
iter | 包含了执行迭代的次数 |
ault | 包含了专家诊断 |
首先,在正态分布中生成了一百组随机数,并且将随机数分配给矩阵x,如下所示:
>x <- rbind(matrix(rnorm(100, sd = 0.3), ncol = 2),
matrix(rnorm(100, mean = 1, sd = 0.3), ncol = 2))
我们可以显示生成的数值,如下所示:
>x
[,1] [,2]
[1,] 0.4679569701 -0.269074028
[2,] -0.5030944919 -0.393382748
[3,] -0.3645075552 -0.304474590
…
[98,] 1.1121388866 0.975150551
[99,] 1.1818402912 1.512040138
[100,] 1.7643166039 1.339428999
生成的kmeans对象值可确定并显示如下(使用10个集群):
> fit <- kmeans(x,10)
> fit
K-means clustering with 10 clusters of sizes 4, 12, 10, 7, 13, 16, 8,
13, 8, 9
Cluster means:
[,1] [,2]
1 0.59611989 0.77213527
2 1.09064550 1.02456563
3 -0.01095292 0.41255130
4 0.07613688 -0.48816360
5 1.04043914 0.78864770
6 0.04167769 -0.05023832
7 0.47920281 -0.05528244
8 1.03305030 1.28488358
9 1.47791031 0.90185427
10 -0.28881626 -0.26002816
Clustering vector:
[1] 7 10 10 6 7 6 3 3 7 10 4 7 4 7 6 7 6 6 4 3 10
4 3 6 10 6 6 3 6 10 3 6 4 3 6 3 6 6 6 7 3 4 6 7 6
10 4 10 3 10 5 2 9 2
[55] 9 5 5 2 5 8 9 8 1 2 5 9 5 2 5 8 1 5 8 2 8
8 5 5 8 1 1 5 8 9 9 8 5 2 5 8 2 2 9 2 8 2 8 2 8
9
Within cluster sum of squares by cluster:
[1] 0.09842712 0.23620192 0.47286373 0.30604945 0.21233870 0.47824982
0.36380678 0.58063931 0.67803464 0.28407093
(between_SS / total_SS = 94.6 %)
Available components:
[1] "cluster" "centers" "totss" "withinss" "tot.
withinss" "betweenss" "size" "iter" "ifault"
如果我们查看结果,就会发现某些有意思的数据点:
我们任意选取数量为10的集群,但是应核实此数字确实好用。通过使用一系列的集群大小多次运行kmeans函数,最终会得到一张与下列示例中图表相似的图表。
例如,如果我们运行下列代码,并且记录结果,输出数据如下所示:
results <- matrix(nrow=14, ncol=2, dimnames=list(2:15,c("clusters","s
umsquares")))
for(i in 2:15) {
fit <- kmeans(x,i)
results[i-1,1] <- i
results[i-1,2] <- fit$totss
}
plot(results)
数据分布得越多,就越能够明确划分集群的最大数量,正如深层聚类的平方和不会有所增加。然而,由于测试中使用了平滑数据,因此可以增加集群数量。
一旦确定了集群,您就能够收集到可视化表示形式,如下图所示:
K-medoids聚类是确定数据集中集群的另一种方法。medoid是数据集的实体。数据集表示将medoid嵌入的群组。K-means与质心一起工作。质心是人为创造的概念,用以表示集群。所以,实际上,medoid是数据集的一部分,质心是派生数量。
当围绕medoids进行分割时,务必确保注意下列事项。
算法的两个阶段包括多个步骤,具体如下。
(1)随机选取成为medoids的k实体(可向算法提供k实体)。
(2)计算相异度矩阵(计算数据集中观测值之间的所有成对相异点(距离))以便发现距离。
(3)将每个实体分配给距离最近的medoid。
(1)搜索实体的每个集群,实体大大降低了平均相异系数,因此成为集群的medoid。
(2)如果medoid出现变化,请从构建阶段的步骤3重新开始。
用pam函数在R编程中进行K-medoid聚类:
pam(x, k, diss, metric, medoids, stand, cluster.only, do.swap, keep.
diss, keep.data, trace.lev)
下表对pam函数的不同参数进行了说明。
参数 | 描述 |
---|---|
x | 这是数据矩阵或相异度矩阵(基于diss标记) |
k | 这是集群数量,其实0比k小,而k比实体数量小 |
diss | 数值如下所示: ● 如果x是矩阵,就为“FALSE” ● 如果x是相异度矩阵,就为“TRUE” |
metric | 这是用于计算相异度矩阵的字符串度量,可以是下列类型之一: ● 用于欧几里得距离的euclidean ● 用于曼哈顿距离的manhattan |
medoids | 如果分配到了“NULL”,就需要开发一组medoids;否则,这就是一组初步medoids |
stand | 如果x是数据矩阵,在计算相异度矩阵之前,会使x的度量标准化 |
cluster.only | 如果数值集为“TRUE”,仅需要计算并返回聚类 |
do. swap | 包含了布尔值,用于决定是否应该进行交换 |
keep. diss | 包含了布尔值,用于决定相异点是否应该保存在结果中 |
keep. data | 包含了布尔值,用于决定数据是否应该在结果中保留 |
trace. lev | 包含了诊断的整体跟踪级别,其中0表示无跟踪信息 |
可显示从pam函数返回的结果,这样难以说明结果;或可绘图表示结果,这样能够直观地了解结果。
使用带有两个(直观上)清晰集群的一组简单数据,如下表所示,此数据存储在名为medoids. csv的文件中:
对象 | x | y |
---|---|---|
1 | 1 | 10 |
2 | 2 | 11 |
3 | 1 | 10 |
4 | 2 | 12 |
5 | 1 | 4 |
6 | 3 | 5 |
7 | 2 | 6 |
8 | 2 | 5 |
9 | 3 | 6 |
使用文件medoids. csv上的pam函数,如下所示:
# load pam function
> library(cluster)
#load the table from a file
> x <- read.table("medoids.csv", header=TRUE, sep=",")
#execute the pam algorithm with the dataset created for the example
> result <- pam(x, 2, FALSE, "euclidean")
Looking at the result directly we get:
> result
Medoids:
ID Object x y
[1,] 2 2 2 11
[2,] 7 7 2 6
Clustering vector:
[1] 1 1 1 1 2 2 2 2 2
Objective function:
build swap
1.564722 1.564722
Available components:
[1] "medoids" "id.med" "clustering" "objective" "isolation"
[6] "clusinfo" "silinfo" "diss" "call" "data"
对结果进行评估:
让我们通过摘要函数来进行更清楚的观察,可见下列结果:
> summary(result)
Medoids:
ID Object x y
[1,] 2 2 2 11
[2,] 7 7 2 6
Clustering vector:
[1] 1 1 1 1 2 2 2 2 2
Objective function:
build swap
1.564722 1.564722
Numerical information per cluster:
sizemax_dissav_diss diameter separation
[1,] 4 2.236068 1.425042 3.741657 5.744563
[2,] 5 3.000000 1.676466 4.898979 5.744563
Isolated clusters:
L-clusters: character(0)
L*-clusters: [1] 1 2
Silhouette plot information:
cluster neighbor sil_width
2 1 2 0.7575089
3 1 2 0.6864544
1 1 2 0.6859661
4 1 2 0.6315196
8 2 1 0.7310922
7 2 1 0.6872724
6 2 1 0.6595811
9 2 1 0.6374808
5 2 1 0.5342637
Average silhouette width per cluster:
[1] 0.6903623 0.6499381
Average silhouette width of total data set:
[1] 0.6679044
36 dissimilarities, summarized :
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.4142 2.3961 6.2445 5.2746 7.3822 9.1652
Metric : euclidean
Number of objects : 9
Available components:
[1] "medoids" "id.med" "clustering" "objective" "isolation"
[6] "clusinfo" "silinfo" "diss" "call" "data"
下载示例代码
凡购买过Packt出版的书籍的用户,均可访问http://www.packtpub.com,登录账户,并下载书籍的示例代码文件。如果读者是在别处购买的本书,可登录http://www.packtpub.com/support,注册账户,以便通过邮箱直接获取文件。
摘要提供了有关medoid的更多详情,并介绍了如何选取medoid。然而,同样需注意相异度。
通过用图表表示数据,可见下列输出数据:
#plot a graphic showing the clusters and the medoids of each cluster
> plot(result$data, col = result$clustering)
生成图和预期一样。从图中能够清楚地看到分裂为两个medoid的数据,各有各的空间,并且用颜色区分二者。
分层聚类是确定数据集中某一层集群的方法。
通过使用分层聚类,我们试图创建集群的层次,有如下两种方法。
通常,通过使用系统树图的树形模型/图模型,显示生成的分层。
用hclust函数在R编程中进行分层聚类。
hclust函数调用如下所示:
hclust(d, method = "complete", members = NULL)
下表对hclust函数的不同参数进行了说明。
参数 | 描述 |
---|---|
d | 这是矩阵 |
method | 这是使用的附聚法。应是以下方法(不同的缩写)之一:ward.D、ward.D2、single、complete、average(= UPGMA)、mcquitty(= WPGMA)、median(= WPGMC)或centroid(= UPGMC) |
members | 可为“NULL”或d,即相异度矩阵 |
开始,我们通过使用下列代码在正态分布上生成某些随机数据:
> dat <- matrix(rnorm(100), nrow=10, ncol=10)
> dat
[,1] [,2] [,3] [,4] [,5]
[,6]
[1,] 1.4811953 -1.0882253 -0.47659922 0.22344983 -0.74227899
0.2835530
[2,] -0.6414931 -1.0103688 -0.55213606 -0.48812235 1.41763706
0.8337524
[3,] 0.2638638 0.2535630 -0.53310519 2.27778665 -0.09526058
1.9579652
[4,] -0.50307726 -0.3873578 -1.54407287 -0.1503834
Then, we calculate the hierarchical distribution for our data as
follows:
> hc <- hclust(dist(dat))
> hc
Call:
hclust(d = dist(dat))
Cluster method : complete
Distance : euclidean
Number of objects : 10
生成的数据对象未提供任何信息。我们通过使用系统树图显示分层集群,如下所示:
>plot(hc)
系统树图是预期的形状。我发现这些图表有些不易了解,但是,如果仔细研究,可得出如下推论。
期望最大化(EM)是评估统计模型中参数的过程。
对于给定的模型,参数如下所示。
执行期望最大化的步骤如下。
(1)将未知参数(T)初始化为随机值。
(2)通过使用新的参数值计算最优缺失值(Z)。
(3)使用刚刚计算出的最优缺失值(Z)以便更好地评估未知参数(T)。
(4)迭代采用步骤2和步骤3,直至收敛。
此版本的算法生成了硬参数值(Z)。实际上,可能对软数值感兴趣,因为将概率分配给了参数的各种数值(Z)。通过硬值,我们可以选取特定的Z值。相反,我们可能会使用软值,因为通过某些概率分布Z会出现变化。
利用mclust程序库里的Mclust函数在R编程中使用EM。Mclust的完整描述是通过EM算法拟合的正态混合模型,EM算法是基于模型的聚类、分类及密度估计,包括贝叶斯正则化。
Mclust函数如下所示:
Mclust(data, G = NULL, modelNames = NULL,
prior = NULL, control = emControl(),
initialization = NULL, warn = FALSE, ...)
下表对Mclust函数的不同参数进行了说明。
参数 | 描述 |
---|---|
data | 包含了矩阵 |
G | 包含了使用的集群数量的向量,用于计算BIC,默认值为1:9 |
modelNames | 包含了使用的模型名称的向量 |
prior | 包含了平均值的可选共轭先验 |
control | 包含了EM的控制参数列表,默认值为List |
initialization | 包含了“NULL”或包含一个或更多下列分量的列表: ● hcPairs:用于合并各组 ● subset:在初始化期间使用 ● noise:初次推测噪声 |
warn | 包含了发出哪个警告,默认为无 |
当Mclust函数试图决定哪个项目属于某一个集群时,此函数就会使用模型。单变量、多变量及单一分量数据集有不同的模型名称。目的在于选取说明数据的模型,例如:VII将用于数据,并且用每个集群的等体积取代进行球状取代。
模型 | 数据集类型 |
---|---|
单变量混合 | |
E | 等方差(一维) |
V | 变量方差(一维) |
多变量混合 | |
EII | 球形,同等体积 |
VII | 球形,不等体积 |
EEI | 对角线,同等体积和形状 |
VEI | 对角线,不同体积、同等形状 |
EVI | 对角线,同等体积、不同形状 |
VVI | 对角线,不同体积和形状 |
EEE | 椭圆,同等体积、形状和方位 |
EEV | 椭圆,同等体积和同等形状 |
VEV | 椭圆,等形状 |
VVV | 椭圆,不同体积、形状和方位 |
单一分量 | |
X | 单变量正态 |
XII | 球形多变量正态 |
XXI | 对角线多变量正态 |
XXX | 椭圆多变量正态 |
首先,我们必须加载包括mclust函数的程序库(我们需要将程序库安装在本地环境中),如下所示:
> install.packages("mclust")
> library(mclust)
我们会在此示例中使用iris数据,如下所示:
> data <- read.csv("http://archive.ics.uci.edu/ml/machine-learningdatabases/
iris/iris.data")
现在,我们可以通过EM计算最优匹配(注意:大写Mclust),如下所示:
> fit <- Mclust(data)
结果显示如下所示:
> fit
'Mclust' model object:
best model: ellipsoidal, equal shape (VEV) with 2 components
> summary(fit)
----------------------------------------------------
Gaussian finite mixture model fitted by EM algorithm
----------------------------------------------------
Mclust VEV (ellipsoidal, equal shape) model with 2 components:
log.likelihood n df BIC ICL
-121.1459 149 37 -427.4378 -427.4385
Clustering table:
1 2
49 100
简单地显示fit数据对象不会提供太多信息,仅表示使用什么计算数据集的密度。
summary指令显示了有关结果的更多详细信息,如下所列:
为了能够直观地确认结果,可用图表表示结果,如下所示:
> plot(fit)
注意:EM 的plot指令生成了下列四个图。
下图是密度图。
第一幅图用不同的模型名称说明BIC范围和分量的数量。在这种情况下,我们可能不应该使用VEV,例如:
第二幅图表示对使用数据馈送的每个分量与数据馈送的每隔一个分量进行比较以便确定产生的聚类。选取提供数据最优聚类的分量。熟悉数据是选取聚类的合适数据点的关键,这只是其中一例。
在此示例中,我认为选取X5.1和X1.4会产生距离最近的集群,如下图所示。
第三幅图表示不同选择对聚类迭代的影响,通过消除图中的任意点来突出主要集群,此图用于主要集群,如下所示。
第四幅图是每个集群的轨道图,相对于每个集群的中心,突出显示点可能会出现在哪个地方,如下所示。
密度估计是评估观测组中群组的概率密度函数的过程。密度估计过程会采用观测值,将其分散到若干数据点中,运行FF转换以便确定核心,然后运行线性近似来评估密度。
密度估计对难以察觉的总体分布函数进行预估。用于生成密度估计的方法如下所示。
通过R编程中的density函数进行密度估计。R中其他用于密度估计的函数如下所示。
函数 | 描述 |
---|---|
DBSCAN | 此函数用于确定固定点集群的聚类 |
OPTICS | 此函数用于确定广泛分布集群的聚类 |
density函数调用如下所示:
density(x, bw = "nrd0", adjust = 1,
kernel = c("gaussian", "epanechnikov",
"rectangular",
"triangular", "biweight",
"cosine", "optcosine"),
weights = NULL, window = kernel, width,
give.Rkern = FALSE,
n = 512, from, to, na.rm = FALSE, ...)
下表对density函数的不同参数进行了说明。
参数 | 描述 |
---|---|
x | 矩阵 |
bw | 使用的平滑带宽 |
adjust | 这是倍增器,用于调节带宽 |
kernel | 使用的是更加平滑的核心,必须是下列核心之一: ● gaussian ● rectangular ● triangular ● epanechnikov ● biweight ● cosine ● optcosine |
weights | 与x具有相同长度的观测权的向量 |
window | 使用的核心 |
width | S兼容性参数 |
give.Rkern | 如果此参数的值为“TRUE”,就表示未预估密度 |
N | 预估的密度点数 |
from, to | 使用的最左边和最右边的点 |
na.rm | 如果此参数的值为“TRUE”,就表示移除缺失值 |
通过使用以下指令可发现可用带宽:
bw.nrd0(x)
bw.nrd(x)
bw.ucv(x, nb = 1000, lower = 0.1 * hmax, upper = hmax, tol = 0.1 *
lower)
bw.bcv(x, nb = 1000, lower = 0.1 * hmax, upper = hmax, tol = 0.1 *
lower)
bw.SJ(x, nb = 1000, lower = 0.1 * hmax, upper = hmax, method =
c("ste", "dpi"), tol = 0.1 * lower)
下表对bw函数的不同参数进行了说明。
参数 | 描述 |
---|---|
x | 数据集 |
nb | 二元数量 |
lower, upper | 需要最小化的带宽范围 |
method | ste方法用于解方程,或者dpi方法用于直接插件 |
tol | ste的收敛公差 |
iris数据集可使用如下:
> data <- read.csv("http://archive.ics.uci.edu/ml/machine-learningdatabases/
iris/iris.data")
The density of the X5.1 series (sepal length) can be computed as
follows:
> d <- density(data$X5.1)
> d
Call:
density.default(x = data$X5.1)
Data: data$X5.1 (149 obs.); Bandwidth 'bw' = 0.2741
x y
Min.:3.478 Min. :0.0001504
1st Qu.:4.789 1st Qu.:0.0342542
Median :6.100 Median :0.1538908
Mean :6.100 Mean :0.1904755
3rd Qu.:7.411 3rd Qu.:0.3765078
Max. :8.722 Max. :0.3987472
密度值可绘制如下:
> plot(d)
此图表示大部分数据都是出现在5和7之间。所以,萼片长度平均低于6。
我们可以使用R编程来检测数据集中的异常。异常检测可用于入侵检测、欺诈检测、系统健康状态等不同领域。在R编程中,这些被称为异常值。R编程允许用多种方法对异常值进行检测:
R编程存在可以显示异常值的函数:identify (in boxplot)。
boxplot函数生成了一个盒须图(请看下图)。boxplot函数有若干图形选项,比如此示例,我们无需进行任何设置。
identify函数是便于标记散点图中点的方法。在R编程中,箱线图是散点图的一种。
在此示例中,我们需要生成100个随机数,然后将盒中的点绘制成图。
然后,我们用第一个异常值的标识符来对其进行标记:
> y <- rnorm(100)
> boxplot(y)
> identify(rep(1, length(y)), y, labels = seq_along(y))
注意图中接近异常值的0。
boxplot函数同样也会自动计算数据集的异常值。
首先,我们会生成100个随机数(注意:此数据是随机生成的,所以结果可能不同):
> x <- rnorm(100)
可通过使用下列代码查看摘要信息:
> summary(x)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-2.12000 -0.74790 -0.20060 -0.01711 0.49930 2.43200
现在,我们可以通过使用下列代码显示异常值:
> boxplot.stats(x)$out
[1] 2.420850 2.432033
下列代码会用图表表示数据集,并且突出显示异常值:
> boxplot(x)
注意图中接近异常值的0。
我们可以通过使用汽车的内置数据生成含有更常见数据的箱线图,这些数据与异常值存在相同的问题,如下所示:
boxplot(mpg~cyl,data=mtcars, xlab="Cylinders", ylab="MPG")
当有两个维度时,我们同样也可以使用箱线图的异常值检测。注意:我们通过使用x和y中异常值的并集而非交集解决问题。此示例就是要显示这样的点。代码如下所示:
> x <- rnorm(1000)
> y <- rnorm(1000)
> f <- data.frame(x,y)
> a <- boxplot.stats(x)$out
> b <- boxplot.stats(y)$out
> list <- union(a,b)
> plot(f)
> px <- f[f$x %in% a,]
> py <- f[f$y %in% b,]
> p <- rbind(px,py)
> par(new=TRUE)
> plot(p$x, p$y,cex=2,col=2)
虽然R确实按照我们的要求去做了,但是此图看起来不对。我们完全是在编造数据,在真正的use case中,需要利用专业知识以便确定这些异常值是否正确。
考虑到构成异常的多样性,R编程带有可以让您完全控制异常的机制:编写能够用于做决策的函数。
我们可使用name函数创建异常,如下所示:
name <- function(parameters,…) {
# determine what constitutes an anomaly
return(df)
}
这里,参数是我们需要在函数中使用的数值。假设我们将函数返回为数据框。利用此函数可以完成任何工作。
我们会在此示例中使用iris数据,如下所示:
> data <- read.csv("http://archive.ics.uci.edu/ml/machine-learningdatabases/
iris/iris.data")
如果我们决定当萼片低于4.5或高于7.5时存在异常,就使用下列函数:
> outliers <- function(data, low, high) {
> outs <- subset(data, data$X5.1 < low | data$X5.1 > high)
> return(outs)
>}
然后,我们会得出下列输出数据:
> outliers(data, 4.5, 7.5)
X5.1 X3.5 X1.4 X0.2 Iris.setosa
8 4.4 2.9 1.4 0.2 Iris-setosa
13 4.3 3.0 1.1 0.1 Iris-setosa
38 4.4 3.0 1.3 0.2 Iris-setosa
42 4.4 3.2 1.3 0.2 Iris-setosa
105 7.6 3.0 6.6 2.1 Iris-virginica
117 7.7 3.8 6.7 2.2 Iris-virginica
118 7.7 2.6 6.9 2.3 Iris-virginica
122 7.7 2.8 6.7 2.0 Iris-virginica
131 7.9 3.8 6.4 2.0 Iris-virginica
135 7.7 3.0 6.1 2.3 Iris-virginica
为了获得预期的结果,可以通过向函数传送不同的参数值灵活地对准则进行轻微调整。
另一个受欢迎的功能包是DMwR。它包括同样可以用于定位异常值的lofactor函数。通过使用以下指令可对DMwR功能包进行安装:
> install.packages("DMwR")
> library(DMwR)
我们需要从数据上移除“种类”列,原因在于其是根据数据进行分类的。可通过使用以下指令完成移除:
> nospecies <- data[,1:4]
现在,我们确定框中的异常值:
> scores <- lofactor(nospecies, k=3)
然后,我们查看异常值的分布:
> plot(density(scores))
一个兴趣点在于:在若干异常值中(即密度约为4)是否存在非常接近的等式。
关联规则说明了两个数据集之间的关联。此规则常用于购物篮分析。一组事务中的每个事务(购物袋)可能包含多个不同项目,那么如何能够让产品销售有关联呢?常见关联如下所示。
在关联规则中,R中广泛使用的工具是apriori。
可调用apriori规则的程序库,如下所示:
apriori(data, parameter = NULL, appearance = NULL, control = NULL)
下表对apriori程序库的不同参数进行了说明。
参数 | 描述 |
---|---|
data | 事务数据 |
parameter | 存储了挖掘的默认行为,其中支持度为0.1、置信度为0.8、最大长度为10,可相应地改变参数值 |
appearance | 用于限制规则中出现的项目 |
control | 用于调整所用算法的性能 |
需要加载apriori规则的程序库,如下所示:
> install.packages("arules")
> library(arules)
加载购物篮数据:
> data <- read.csv("http://www.salemmarafi.com/wp-content/
uploads/2014/03/groceries.csv")
然后,我们可以从数据中生成规则:
> rules <- apriori(data)
parameter specification:
confidenceminvalsmaxaremavaloriginalSupport support minlenmaxlen
target
0.8 0.1 1 none FALSE TRUE 0.1 1
10 rules
ext
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 ...[655 item(s), 15295 transaction(s)] done [0.00s].
sorting and recoding items ... [3 item(s)] done [0.00s].
creating transaction tree ... done [0.00s].
checking subsets of size 1 2 3 done [0.00s].
writing ... [5 rule(s)] done [0.00s].
creating S4 object ... done [0.00s].
在结果中有若干个突出的点:
我们可以对已生成的规则进行检查,如下所示:
> rules
set of 5 rules
> inspect(rules)
lhsrhs support confidence lift
1 {semi.finished.bread=} => {margarine=} 0.2278522 1
2.501226
2 {semi.finished.bread=} => {ready.soups=} 0.2278522 1
1.861385
3 {margarine=} => {ready.soups=} 0.3998039 1
1.861385
4 {semi.finished.bread=,
margarine=} => {ready.soups=} 0.2278522 1
1.861385
5 {semi.finished.bread=,
ready.soups=} => {margarine=} 0.2278522 1
2.501226
为了方便阅读,对代码的格式做了微调整。
仔细观察这些规则,发现在市场上购买面包、汤及人造黄油之间确实存在关系,至少在市场收集数据的时候是这样。
如果我们改变了计算中使用的参数(阈值),就会得到一组不同的规则。例如,核实下列代码:
> rules <- apriori(data, parameter = list(supp = 0.001, conf = 0.8))
此代码生成了500多个规则,但是这些规则的含义值得质疑,原因在于,目前规则的置信度为0.001。
本章探讨了聚类分析、异常检测及关联规则。“聚类分析”部分使用了K-means聚类、K-medoids聚类、分层聚类、期望最大化及密度估计。“异常检测”部分通过使用内置R函数发现了异常值,并开发出自己的专用R函数。“关联规则”部分使用了apriori功能包,用以确定数据集中的关联。
下一章我们会涉及序列的数据挖掘。
数据挖掘技术一般用于检测数据中的序列或模式。本章中,我们将试图使数据遵循一种模式,在此模式中,一个或一系列事件可以通过一致的方式预测另一个数据点。
本章描述了在数据集中查找模式的不同方法:
我们可以在很多较大的数据集中查找模式。这涵盖了很多区域,比如混合人口的变化、使用手机的频率、高速公路质量衰退、年龄因素造成的事故等。不过我们能明确地感受到,有很多模式和序列正等待我们去发现。
我们可以通过使用R编程中的一些工具找到这些模式。大多数模式因约束条件而在一定程度上受到限制,如序列的有用时间。
我们来回顾一下确定数据中模式的方法:
模型类别 | 模型如何工作 |
---|---|
eclat | 此模型用于项集模式检测,其中购物车最为常见 |
arules | 此模型确定数据集中的项目共现 |
apriori | 此模型学习数据集中的关联规则 |
TraMineR | 这是一个用于挖掘序列的R功能包 |
Eclat算法用于频繁项集的挖掘。这种情况下,我们寻找行为相似的模式,与之相对的是寻找不规则模式(与处理其他数据挖掘的方法类似)。
Algorithm通过数据中的交集来估算同时频繁出现事件候选项(如购物车项目)的支持度。然后通过对频繁候选项进行测试来证实数据集中的模式。
在R编程中使用Eclat就是使用arules功能包中的eclat函数。使用Eclat算法的R编程遵循了此处提出的约定:
> eclat(data,
parameter = NULL,
control = NULL)
下表对eclat函数的不同参数进行了说明:
参数 | 描述 |
---|---|
data | 待分析的数据矩阵 |
parameter | ECParameter或列表的对象 |
control | ECControl或列表的对象 |
常见的ECParameter如下所示:
参数 | 描述 |
---|---|
support | 此参数界定了一个项集的最小支持度(默认值为0.1) |
minlen | 此参数包含了一个项集的最小容量(默认值为1) |
maxlen | 此参数包含了一个项集的最大容量(默认值为10) |
target | 此参数界定了待挖掘关联项集的类型: ● 频繁项集 ● 最频繁项集 ● 频繁闭项集 |
ECControl常见值如下所示:
参数 | 描述 |
---|---|
sort | 此参数为下列数值之一: ● 1代表升序 ● -1代表降序 ● 0代表未分类 ● 2代表升序 ● -2代表事务容量和的降序 |
verbose | 此参数显示进度信息 |
调用eclat函数返回数据中出现的频繁项集。
Eclat的实施包括成年人数据集。成年人数据集包括人口统计局数据中约50000行的数据。
使用下列代码找到成年人行为的相似点:
> 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
检查最终结果时,我们会注意到以下细节:
处理数据集时,通过下列代码挖掘出现的频繁项集:
> 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
以下是对之前输出数据的研究所得:
为了进一步证实数据,我们可以将范围缩减至数据集中出现的最高频率(可以通过调节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
按照预期,由一位美国本土且拥有工作的已婚男士填写普查数据表格。
在R中,arulesNBMiner是一个功能包,用于寻找一个集合中两个或两个以上项目的共现。底层模型,即负二项式模型,允许高度偏态次数分配,否则会很难确定最小项集容量。我们在正被挖掘的较大数据集中寻找频繁数据集。当确定使用arulesNBMiner时,您应该看到一些迹象:项目集频率正出现在数据子集合中。
将arulesNBMiner作为功能包进行操作,并且必须将此功能包安装于您的R编程环境中。可以通过使用任意数据集来学习如何使用模型/函数中包含的工具,如下所示:
> results <-NBMiner(data, parameter, control = NULL)
下表对NBMiner函数的不同参数进行了说明:
参数 | 描述 |
---|---|
data | 待分析的数据矩阵 |
parameter | 参数列表(自动转换为NBMinerParameters的对象) |
control | 使用的控制列表(自动转换为NBMinerControl的对象),目前仅verbose和调试逻辑可用 |
NBMinerParameters是用于调用NBMiner的参数块,其架构如下所示:
NBMinerParameters(data, trim = 0.01, pi = 0.99,
theta = 0.5, minlen = 1, maxlen = 5, rules = FALSE,
plot = FALSE, verbose = FALSE, getdata = FALSE)
NBMinerParameters的数值如下所示:
参数 | 描述 |
---|---|
data | 事务 |
trim | 从数据频率分布尾数修剪的分数 |
pi | 精度阈值π |
theta | 剪枝参数θ |
minlen | 项集中所发现项目的最小数(默认值为1) |
maxlen | 项集中所发现项目的最大数(默认值为5) |
rules | 包含了布尔值,用于确定是否挖掘NB精确规则而非NB频繁项集 |
plot | 包含了布尔值,用于确定是否为模型绘图 |
verbose | verbose输出参数,用于估算程序 |
getdata | 用于获取研究及估算的计数 |
功能包中的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
以下是对之前输出数据的研究所得:
如果以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
以下是对之前输出数据的研究所得:
Apriori是可以帮助了解关联规则的分类算法。与事务的实施方式相对。这种算法尝试找到数据集中常见的子集合,必须满足最小阈值以便核实关联。
Apriori的支持度和置信度概念十分有趣。Apriori方法会从您的数据集中返回有趣的关联,如当出现Y时,会返回X。支持度是包含X和Y的事务的百分比。置信度是同时包含X和Y的事务的百分比。支持度的默认值为10,置信度的默认值为80。
apriori方法的使用规则如下所示:
apriori(data, parameter = NULL, appearance = NULL, control = NULL)
下表对apriori函数的不同参数进行了说明:
参数 | 描述 |
---|---|
data | 这是所需的数据集 |
parameter | 这是用于控制处理方式的参数列表。支持度的默认值为0.1,置信度默认值为0.8,maxlen默认值为10 |
appearance | 控制了所用的数据值 |
control | 控制了算法,特别是分类的效能 |
我们正寻找食品超市中典型购物篮内购买的项目之间的关联。为此,我们将按下列步骤进行操作。
(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
以下是对之前输出数据的研究所得:
(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
一条派生规则将这些措施的置信度展现得非常全面。
TraMineR功能包用于挖掘序列,并将其可视化,其思想是发现序列。可以将序列分布、序列频率及湍流等绘图的图解设备构建到功能包中。此外,还有一些自然出现的项目,其中的数据有重复的序列,如在一些社会科学场地,数据会自然地循环项目。
通过此文件,我将带您大概了解TraMineR,以便生成一系列用于发现序列的工具。在挖掘操作中选择何种工具取决于您自己。
您可以同时使用TraMineR功能包及一对内置数据集:
数据集 | 描述 |
---|---|
actcal | 此数据集包含了2000年以来每月活动单个的状态符号 |
biofam | 此数据集包含了15岁至30岁期间单个的家庭生活状态 |
mvad | 此数据集包含了每月活动单个的状态数据 |
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函数的不同参数进行了说明:
参数 | 描述 |
---|---|
data | 矩阵 |
var | 会有一个纵列列表,其中包含序列或代表所有纵列都存在的“NULL” |
informat | 包含了原始数据的格式,可以是下列格式中的任意一种: ● STS ● SPS ● SPELL |
stsep | 分隔符 |
alphabet | 所有可能状态的列表 |
states | 包含了短时状态的标记 |
在这一示例中,我们将看到人们生活中从训练到工作的进程中时间的序列。我们期望看到从失业未经训练的状态至经过训练并最终成为全职员工的进程。
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(myseq)
通过参照个人不同状态间界定的转换期,您会发现连续几个月都有训练。您应进行核实,以便数据显示的信息与您对序列数据的理解相一致。
现在,以频率图表为例:
> seqfplot(myseq)
现在我们来看序列在不同时间的频率。多次观看后我们会看到同一序列的人群集,如经过一段时间的训练后会有工作。
现在,以分布图表为例:
> seqdplot(myseq)
我们来看看序列状态在不同时期的分布情况。通常情况下,人们在上学或训练后开始工作。明白了吧!
通过以下指令我们可以看到序列的熵:
> seqHtplot(myseq)
熵在不同时期的变化特点:明显降低后会出现细微的上升。这与不同人群会在最初做出不同选择的情况一致(很多状态),如上学或训练,然后进行工作,成为劳动力(一种状态)。
有一个有趣的想法为数据湍流。湍流传达出一个信息,即从数据中可见的某个特定事例可以推导出多少不同的后续序列。我们可以用seqST函数将湍流可视化。seqST函数将序列数据作为参数,并返还湍流数据。让我们继续以示例进行说明:
> myturbulence <- seqSt(myseq)
> hist(myturbulence)
我们可以看到带有长尾数的近乎标准化分布。大多数状态分为少量后续状态以及少数状态或多或少的异常值。
TraMineR功能包也可以确定序列的度量,如不同序列间的相异点。
使用TraMineR中的seqdist函数可以实现所有这些功能。
我们可以用seqdist计算LCP。
seqdist函数使用规则如下:
seqdist(seqdata, method, refseq=NULL, norm=FALSE,
indel=1, sm=NA, with.missing=FALSE, full.matrix=TRUE)
下表对seqdist函数的不同参数进行了说明:
参数 | 描述 |
---|---|
seqdata | 这是状态序列(用seqdef界定) |
method | 包含了待用的LCP方法 |
refseq | 这是可选参考序列 |
norm | 将距离标准化 |
indel | 仅用于OM |
sm | 这是替代矩阵(忽略LCP) |
with.missing | 如果出现缺失间隙,数值为“TRUE” |
full.matrix | 如果数值为“TRUE”,返回全矩阵 |
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(仅用较小序列来诠释概念)。
本章,我们探讨了确定数据序列的不同模式。通过使用eclat函数查找数据集模式,以便寻找人口中的相似模式。使用TraMineR查找购物篮中的项目频集。使用apriori规则确定购物篮中的项目关联。使用TraMineR确定成年人职业转换期的序列,并通过序列数据可用的大量图形特征将其可视化。最后,用seqdist检查序列之间的相似点和不同点。
下一章我们将探讨文本挖掘或基于文本的数据集,而非数值形式或分类属性。