数据科学:R语言实战

978-7-115-43590-3
作者: 【美】Dan Toomey(丹·图米)
译者: 刘丽君李成华卢青峰
编辑: 王峰松
分类: R语言

图书目录:

详情

这本书旨在探索数据科学家感兴趣的核心话题。它收集了各种数据源,并使用公开可用的R函数和R包评估这些数据。同时,在本书中,读者还将了解业界经常使用的分析技术。最后,读者将学会如何用R语言来实现一系列数据科学技术。

图书摘要

版权信息

书名:数据科学:R语言实战

ISBN:978-7-115-43590-3

本书由人民邮电出版社发行数字版。版权所有,侵权必究。

您购买的人民邮电出版社电子书仅供您个人使用,未经授权,不得以任何方式复制和传播本书内容。

我们愿意相信读者具有这样的良知和觉悟,与我们共同保护知识产权。

如果购买者有侵权行为,我们可能对该用户实施包括但不限于关闭该帐号等维权措施,并可能追究法律责任。

• 著    [美] Dan Toomey

  译    刘丽君 李成华 卢青峰

  责任编辑 王峰松

• 人民邮电出版社出版发行  北京市丰台区成寿寺路11号

  邮编 100164  电子邮件 315@ptpress.com.cn

  网址 http://www.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无法与相当稀疏的数据或带有许多异常值的数据一起工作。此外,如果集群的线性形状不好,就会出现问题。图示应证明数据是否适合此算法。

1.用法

用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 包含了专家诊断

2.示例

首先,在正态分布中生成了一百组随机数,并且将随机数分配给矩阵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重新开始。

1.用法

用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函数返回的结果,这样难以说明结果;或可绘图表示结果,这样能够直观地了解结果。

2.示例

使用带有两个(直观上)清晰集群的一组简单数据,如下表所示,此数据存储在名为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编程中进行分层聚类。

1.用法

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,即相异度矩阵

2.示例

开始,我们通过使用下列代码在正态分布上生成某些随机数据:

> 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算法是基于模型的聚类、分类及密度估计,包括贝叶斯正则化。

1.用法

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 包含了发出哪个警告,默认为无

2.模型名称列表

当Mclust函数试图决定哪个项目属于某一个集群时,此函数就会使用模型。单变量、多变量及单一分量数据集有不同的模型名称。目的在于选取说明数据的模型,例如:VII将用于数据,并且用每个集群的等体积取代进行球状取代。

模型 数据集类型
单变量混合
E 等方差(一维)
V 变量方差(一维)
多变量混合
EII 球形,同等体积
VII 球形,不等体积
EEI 对角线,同等体积和形状
VEI 对角线,不同体积、同等形状
EVI 对角线,同等体积、不同形状
VVI 对角线,不同体积和形状
EEE 椭圆,同等体积、形状和方位
EEV 椭圆,同等体积和同等形状
VEV 椭圆,等形状
VVV 椭圆,不同体积、形状和方位
单一分量
X 单变量正态
XII 球形多变量正态
XXI 对角线多变量正态
XXX 椭圆多变量正态

3.示例

首先,我们必须加载包括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.1X1.4会产生距离最近的集群,如下图所示。

第三幅图表示不同选择对聚类迭代的影响,通过消除图中的任意点来突出主要集群,此图用于主要集群,如下所示。

第四幅图是每个集群的轨道图,相对于每个集群的中心,突出显示点可能会出现在哪个地方,如下所示。

密度估计是评估观测组中群组的概率密度函数的过程。密度估计过程会采用观测值,将其分散到若干数据点中,运行FF转换以便确定核心,然后运行线性近似来评估密度。

密度估计对难以察觉的总体分布函数进行预估。用于生成密度估计的方法如下所示。

通过R编程中的density函数进行密度估计。R中其他用于密度估计的函数如下所示。

函数 描述
DBSCAN 此函数用于确定固定点集群的聚类
OPTICS 此函数用于确定广泛分布集群的聚类

1.用法

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的收敛公差

2.示例

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编程中,箱线图是散点图的一种。

1.示例1

在此示例中,我们需要生成100个随机数,然后将盒中的点绘制成图。

然后,我们用第一个异常值的标识符来对其进行标记:

> y <- rnorm(100)
> boxplot(y)
> identify(rep(1, length(y)), y, labels = seq_along(y))

注意图中接近异常值的0。

2.示例2

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")

3.另一个有关异常检测的示例

当有两个维度时,我们同样也可以使用箱线图的异常值检测。注意:我们通过使用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编程带有可以让您完全控制异常的机制:编写能够用于做决策的函数。

1.用法

我们可使用name函数创建异常,如下所示:

name <- function(parameters,…) {
  # determine what constitutes an anomaly
  return(df)
}

这里,参数是我们需要在函数中使用的数值。假设我们将函数返回为数据框。利用此函数可以完成任何工作。

2.示例1

我们会在此示例中使用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

为了获得预期的结果,可以通过向函数传送不同的参数值灵活地对准则进行轻微调整。

3.示例2

另一个受欢迎的功能包是DMwR。它包括同样可以用于定位异常值的lofactor函数。通过使用以下指令可对DMwR功能包进行安装:

> install.packages("DMwR")
> library(DMwR)

我们需要从数据上移除“种类”列,原因在于其是根据数据进行分类的。可通过使用以下指令完成移除:

> nospecies <- data[,1:4]

现在,我们确定框中的异常值:

> scores <- lofactor(nospecies, k=3)

然后,我们查看异常值的分布:

> plot(density(scores))

一个兴趣点在于:在若干异常值中(即密度约为4)是否存在非常接近的等式。

关联规则说明了两个数据集之间的关联。此规则常用于购物篮分析。一组事务中的每个事务(购物袋)可能包含多个不同项目,那么如何能够让产品销售有关联呢?常见关联如下所示。

在关联规则中,R中广泛使用的工具是apriori。

1.用法

可调用apriori规则的程序库,如下所示:

apriori(data, parameter = NULL, appearance = NULL, control = NULL)

下表对apriori程序库的不同参数进行了说明。

参数 描述
data 事务数据
parameter 存储了挖掘的默认行为,其中支持度为0.1、置信度为0.8、最大长度为10,可相应地改变参数值
appearance 用于限制规则中出现的项目
control 用于调整所用算法的性能

2.示例

需要加载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通过数据中的交集来估算同时频繁出现事件候选项(如购物车项目)的支持度。然后通过对频繁候选项进行测试来证实数据集中的模式。

1.用法

在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行的数据。

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

检查最终结果时,我们会注意到以下细节:

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

按照预期,由一位美国本土且拥有工作的已婚男士填写普查数据表格。

在R中,arulesNBMiner是一个功能包,用于寻找一个集合中两个或两个以上项目的共现。底层模型,即负二项式模型,允许高度偏态次数分配,否则会很难确定最小项集容量。我们在正被挖掘的较大数据集中寻找频繁数据集。当确定使用arulesNBMiner时,您应该看到一些迹象:项目集频率正出现在数据子集合中。

1.用法

将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

以下是对之前输出数据的研究所得:

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

以下是对之前输出数据的研究所得:

Apriori是可以帮助了解关联规则的分类算法。与事务的实施方式相对。这种算法尝试找到数据集中常见的子集合,必须满足最小阈值以便核实关联。

Apriori的支持度和置信度概念十分有趣。Apriori方法会从您的数据集中返回有趣的关联,如当出现Y时,会返回X。支持度是包含X和Y的事务的百分比。置信度是同时包含X和Y的事务的百分比。支持度的默认值为10,置信度的默认值为80。

1.用法

apriori方法的使用规则如下所示:

apriori(data, parameter = NULL, appearance = NULL, control = NULL)

下表对apriori函数的不同参数进行了说明:

参数 描述
data 这是所需的数据集
parameter 这是用于控制处理方式的参数列表。支持度的默认值为0.1,置信度默认值为0.8,maxlen默认值为10
appearance 控制了所用的数据值
control 控制了算法,特别是分类的效能

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

以下是对之前输出数据的研究所得:

(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 此数据集包含了每月活动单个的状态数据

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函数的不同参数进行了说明:

参数 描述
data 矩阵
var 会有一个纵列列表,其中包含序列或代表所有纵列都存在的“NULL”
informat 包含了原始数据的格式,可以是下列格式中的任意一种:
● STS
● SPS
● SPELL
stsep 分隔符
alphabet 所有可能状态的列表
states 包含了短时状态的标记

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(myseq)

通过参照个人不同状态间界定的转换期,您会发现连续几个月都有训练。您应进行核实,以便数据显示的信息与您对序列数据的理解相一致。

现在,以频率图表为例:

> seqfplot(myseq)

现在我们来看序列在不同时间的频率。多次观看后我们会看到同一序列的人群集,如经过一段时间的训练后会有工作。

现在,以分布图表为例:

> seqdplot(myseq)

我们来看看序列状态在不同时期的分布情况。通常情况下,人们在上学或训练后开始工作。明白了吧!

通过以下指令我们可以看到序列的熵:

> seqHtplot(myseq)

熵在不同时期的变化特点:明显降低后会出现细微的上升。这与不同人群会在最初做出不同选择的情况一致(很多状态),如上学或训练,然后进行工作,成为劳动力(一种状态)。

有一个有趣的想法为数据湍流。湍流传达出一个信息,即从数据中可见的某个特定事例可以推导出多少不同的后续序列。我们可以用seqST函数将湍流可视化。seqST函数将序列数据作为参数,并返还湍流数据。让我们继续以示例进行说明:

> myturbulence <- seqSt(myseq)
> hist(myturbulence)

我们可以看到带有长尾数的近乎标准化分布。大多数状态分为少量后续状态以及少数状态或多或少的异常值。

TraMineR功能包也可以确定序列的度量,如不同序列间的相异点。

使用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函数的不同参数进行了说明:

参数 描述
seqdata 这是状态序列(用seqdef界定)
method 包含了待用的LCP方法
refseq 这是可选参考序列
norm 将距离标准化
indel 仅用于OM
sm 这是替代矩阵(忽略LCP)
with.missing 如果出现缺失间隙,数值为“TRUE”
full.matrix 如果数值为“TRUE”,返回全矩阵

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(仅用较小序列来诠释概念)。

实际问题
什么时候做、怎样做以及为什么这样做?
挑战

本章,我们探讨了确定数据序列的不同模式。通过使用eclat函数查找数据集模式,以便寻找人口中的相似模式。使用TraMineR查找购物篮中的项目频集。使用apriori规则确定购物篮中的项目关联。使用TraMineR确定成年人职业转换期的序列,并通过序列数据可用的大量图形特征将其可视化。最后,用seqdist检查序列之间的相似点和不同点。

下一章我们将探讨文本挖掘或基于文本的数据集,而非数值形式或分类属性。


相关图书

科研论文配图绘制指南——基于R语言
科研论文配图绘制指南——基于R语言
R语言编程:基于tidyverse
R语言编程:基于tidyverse
R语言医学多元统计分析
R语言医学多元统计分析
Python与R语言数据科学实践
Python与R语言数据科学实践
R数据挖掘实战
R数据挖掘实战
R语言机器学习实战
R语言机器学习实战

相关文章

相关课程