生信算法生信笔记

层次聚类分析案例(二)

2019-11-14  本文已影响0人  11的雾

之前的笔记:
聚类介绍:点这里
层次聚类分析案例(一)

案例二:亚马逊雨林烧毁情况

1999~2010年,33000平方英里(85500平方公里),即2.8%的亚马逊雨林被烧毁。这一结果是被NASA领导的研究项目发现的。该研究的主要目的是衡量森林树冠下暗火的蔓延程度。该研究发现火灾烧毁的森林比用于农耕而砍伐的森林面积大很多。然而,森林烧毁情况和火灾之间没有建立起联系。

如何建立火灾和森林烧毁情况之间的联系,需要基于NASA的Aqua卫星上的大气红外探测仪(AIRS)设备的湿度数据。火灾频率与夜间的低湿度相吻合,低湿度使得地表的低强度火灾能够持续燃烧。

准备工作

为了进行层次聚类,我们应该使用采集于亚马逊雨林(1999~2010年)的数据集。

第1步:收集和描述数据

该任务使用名为NASAUnderstory的数据集。该数据以标准格式存储在名为NASAUnderstory.csv的CSV格式的文件中。其中包含64行数据和32个变量。数据在这里

非数值型数据如下:Overstory Species ;Labels 。其余为数值型变量。
具体实施步骤
以下为实现细节。

第2步:探索数据

让我们探索数据并理解变量间的关系。我们从导入名为NASAUnderstory.csv的文件开始,将该数据存储到NASA数据框中:

NASA = read.csv("Data/NASAUnderstory.csv",header = T)

下一步,我们应该获取每一个物种列标签的长版本:

NASA.lab = NASA$Labels

之后,输出NASA.lab数据框。这包含了每一个物种的完整名字。

结果如下:


接下来把整个数据内容传递给NASA数据框:

NASA = NASA[,-32]

输出NASA数据框可以把整体的数据内容显示出来。

NASA

结果如下:


第3步:转换数据

接下来进行数据归一化。scale()函数将中心化并拉伸所有先前提到的数值型变量列:

NASAscale <- scale(NASA[,3:31])

这会拉伸缩NASA数据框中第3~31列的所有数值型变量。

输出NASAscale数据框可以显示所有的拉伸和中心化后的数值。

NASAscale

结果如下:


为了将一个向量编码为一个因子,需要使用factor函数。如果自变量排序是TRUE,因子等级将被排序。这里,我们将OverstorySpecies列作为一个值传给factor函数:

rownames(NASAscale) = as.factor(NASA$Overstory.Species)

as.factor()返回列名的数据框。

输出数据框rownames(NASAscale)将显示OverstorySpecies列的所有值:

rownames(NASAscale)

结果如下:


第4步:训练和评估模型

下一步是训练模型。首先计算距离矩阵。dist()函数用来计算并返回距离矩阵,使用特定的距离度量来计算数据矩阵中各行间的距离。这里使用的距离度量可以是欧式距离、最大距离、曼哈顿距离、堪培拉距离、二进制距离,或闵可夫斯基距离。这里的距离度量使用欧式距离。欧式距离计算两个向量间的距离为sqrt(sum((x_i-y_i)^2))。结果存储在一个新的数据框dist1中。

dist1 <- dist(NASAscale, method = "euclidean")

下一步是使用Ward方法进行聚类。hclust()函数对一组不同的n个对象进行聚类分析。第一阶段,每个对象被指派给它自己的簇。之后每个阶段,算法迭代聚合两个最相似的簇。这个过程不断持续直到只剩一个簇。hclust()函数要求我们以距离矩阵的形式提供数据。dist1数据框被传入。默认使用全链接进行聚类。此外还可以使用不同的聚集方法,包括ward.D、ward.D2、single、complete和average。

clust1 <- hclust(dist1,method = "ward.D")

clust1

调用clust1,结果显示了聚集方法、计算距离的方法,以及对象的数量。结果如下:


第5步:绘制模型

plot()函数是绘制R对象的通用函数。这里,plot()函数用来绘制系统树图:

plot(clust1,labels=NASA[,2], cex=0.5,

     xlab="",ylab="Distance",main="Clustering for NASA Understory Data")

结果如下:


rect.hclust()函数会围绕系统树图的某些枝干绘制矩形以强调对应的簇。系统树图首先在某个等级上被剪切,之后在选定的枝干上绘制长方形。

clust1作为一个对象传入函数,同时传入的还有需要形成的簇的数量。

rect.hclust(clust1,k=2)

结果如下:


cuts()函数基于期望的簇数量或者切割高度切割这棵树中的节点到不同的组。这里,clust1被当作一个对象传入该函数,同时传入期望的簇数量。

cuts = cutree(clust1,k=2)

cuts

结果如下:


第6步:改进模型

首先需要载入以下包:

library(vegan)

vegan包最初被社会学和植被生态学家广泛使用。它包括排序方法、多样性分析以及其他功能。其中一些流行的工具包括多样性分析、物种丰富度模型、物种丰富度分析、差异分析等。

下一步是通过使用jaccard距离方法训练模型来改进模型。第一步是使用vegdist()函数计算距离矩阵。该函数计算成对元素的距离,并把结果存储在一个新的数据框dist1中。jaccard系数衡量有限样本集的相似性,该系数通过两集合的交集大小除以两集合的并集大小计算得来。

dist1 <- vegdist(NASA[,3:31],method="jaccard", upper=T)

下一步是使用Ward方法进行聚类。使用hclust()函数:

clust1 <- hclust(dist1,method = "ward.D")

clust1

调用clust1显示了使用的聚集方法、距离度量方法和对象的个数。结果如下:


plot()函数是绘制R语言对象的通用函数:

plot(clust1,labels=NASA[,2],cex=0.5,

     xlab=" ",ylab="Distance",main="Clustering for NASA Understory Data")

clust1数据框作为对象传入该函数。cex给出了缩放数值,通过该值可以把文字和符号相对放大到默认值。
结果如下:


clust1对象被作为一个对象传入函数,同时传入的还有聚类的数量:

rect.hclust(clust1,k=2)

结果如下:


cuts()函数能够基于期望的组数量或者切割的高度将树切割到不同的组:

cuts = cutree(clust1,k=2)

cuts

结果如下:


使用判别函数绘制两个类的解。

clusplot()函数能够绘制二维的聚类图。这里,把NASA数据框作为对象传入。

clusplot(pam(NASAscale,2))

结果如下:


使用判别函数绘制两个类的解。

为了区分给定的类别,plotcluster()函数使用映射函数进行绘图。不同的映射方法包括经典判别坐标、展现平均值和协方差结构区别的方法、非对称方法(从混合类中分离出同质类)、本地基于近邻的方法和基于鲁棒协方差矩阵的方法。

clusplot()函数可以画出二维的聚类图。这里把NASA的数据框作为对象传入:

clusplot(NASA,cuts,color = TRUE,shade = TRUE, labels = 2,lines = 0,

         main = "NASA Two Cluster Plot, Ward's Method, First two PC")

结果如下:


接下来,对NASAscale数据框用t()函数进行转换:

library(fpc)

NASAtrans=t(NASAscale)

下一步是使用闵可夫斯基距离方法改进模型。第一步是计算距离矩阵。这里使用dist()函数。

闵可夫斯基距离经常在变量为比例尺度以及绝对零值的情况下使用。

dist1 <- dist(NASAtrans, method="minkowski",p=3)

下一步是使用Ward方法进行聚类。使用hclust()方法。

clust1 <- hclust(dist1,method="ward.D")

clust1

调用clust1函数显示了使用的聚集方法、距离度量方法和对象的个数,结果如下:


plot()函数是绘制R对象的通用函数。这里,plot()函数被用来绘制系统树图:

plot(clust1, labels=NASA.lab[1:29], cex=1,
     xlab="",ylab="Distance",main="Clustering for NASA Understory Data")

结果如下:


rect.hclust()函数会围绕系统树图的某些枝干绘制矩形以强调对应的簇。系统树图首先在某个等级上被剪切,之后在选定的枝干上绘制长方形。

clust1对象被当作一个对象传入该函数,同时被传入的还有聚类的个数:

rect.hclust(clust1,k=3)

结果如下:


cuts()函数能够基于期望的簇数量或者切割的高度将树切分到三个不同的组。这里传入clust1对象和聚类的个数。

cuts = cutree(clust1,k=3)

cuts

结果如下:


上一篇下一篇

猜你喜欢

热点阅读