R语言医学数据分析实战（第十章文中源码解析）-CSDN博客

本文链接：https://blog.csdn.net/karol_peng/article/details/124783122

Chapter 10

聚类分析是在事物的分类面貌尚不清楚的情况下讨论分类问题。
分为两种：Q型聚类（样品聚类），R型聚类（观察指标/变量聚类）
具体分类标准：比较样本/指标的相似程度，将相似度大的归为一类，相似度小的归为一类。

10.1 相似性度量

10.1.1 样品间的距离

# set.seed(1234)用于设定随机数种子，括号里的数字只是一个符号而已。后面如果我们想生成相同的随机数字，可以用1234这个符号
# 使用set.seed只是保证结果的可重复性，可再次复现性。
# 参考：https://blog.csdn.net/vencent_cy/article/details/50350020
# rnorm是随机生成符合正态分布的十个数字，nrow表示矩阵的行数
# dist函数用于计算距离，默认的method是"euclidean"，即欧几里得距离。
# diag = TRUE, upper = TRUE表示显示完整的矩阵，而不是一个下三角矩阵
set.seed(1234)
dat <- matrix(rnorm(10), nrow = 5)
dat
dist(dat)
dist(dat, diag = TRUE, upper = TRUE)

# matrix(c(rnorm(5), rnorm(5, mean = 100, sd = 10)), nrow = 5)表示生成5行数据，前五行是rnorm(5)，
# 后五行是rnorm(5, mean = 100, sd = 10))，即均值为100，标准差是10
set.seed(123)
dat <- matrix(c(rnorm(5), rnorm(5, mean = 100, sd = 10)), nrow = 5)
dat
# colMeans是用于计算所有列的平均值，同理rowMeans计算所有行的平均值
# scale函数将数据进行标准化，我理解的就是将数据的标准差都统一为0
# 参考：https://blog.csdn.net/weixin_39541558/article/details/80054229
colMeans(dat)
dist(dat)
dat.scale <- scale(dat)
dist(dat.scale)

# 下面是一个自己编写的函数，用于计算马氏距离。马氏距离不受量纲的影响，避免变量相关性的影响。
# 记住就行
dist.ma <- function (data) {
    X <- as.matrix(na.omit(data)) 
    V <- cov(X) 
    L <- t(chol(V))  
    stdX <- t(forwardsolve(L, t(X))) 
    d <- dist(stdX)  
    attr(d, "Labels") <- row.names(data)
    d
}
dist.ma(dat)

10.1.2 变量间的相关系数

# 这里set.seed(1234)符号是1234，我们会发现生成的数据跟前面生成的dat数据具有重复性。
# 这就是seed的具体用处
set.seed(1234)
x <- rnorm(20)
y <- rnorm(20)
# cor计算样本的相关系数，可参考具体讲解：https://www.jianshu.com/p/fe15c2c28574
# sum(x * y) / sqrt(sum(x ^ 2) * sum(y ^ 2)) 是计算余弦值，可参考书中的计算公式，很容易获得
cor(x, y)
sum(x * y) / sqrt(sum(x ^ 2) * sum(y ^ 2)) 
# 暂时不理解将数据标准化之后获取的夹角余弦与前面的相似系数等价有什么关系？？？
x.scale <- scale(x)
y.scale <- scale(y)
sum(x.scale * y.scale) / sqrt(sum(x.scale ^ 2) * sum(y.scale ^ 2))

10.2 层次聚类法

10.2.1 类之间相关系数的定义

10.2.2 Q型聚类（样本聚类）

# head函数用于查看数据的前几行几列的数据
# tolower函数表示将数据中name的大写转化为小写
install.packages("flexclust")
library(flexclust)
data(nutrient)
head(nutrient)
row.names(nutrient) <- tolower(row.names(nutrient))
# 使用str函数查看变量的类型，summary常用于看数据的分布情况。
# 根据这两个，可以分析数据是否属于同一量纲，若差距太大，需要使用scale进行标准化操作
# hclust函数进行层次聚类（第一个参数是距离矩阵，第二个是methods类型），属于flexclust程序包
# rect.hclust是在树状图的分值周围绘制矩形，突出显示相应的簇
# plot(hc1, hang = - 1)可以看出plot函数的强大，它是一个泛型函数，
# 可以理解为根据不同的数据类型去匹配相应的具体函数，这里的hang=-1是将样品标签在同一水平线显示。
str(nutrient)
summary(nutrient)
nutrient.scaled <- scale(nutrient)
## 方法1：使用欧氏距离计算层次分类
d.eu <- dist(nutrient.scaled, method = "euclidean")
hc1 <- hclust(d.eu, method = "average")
plot(hc1, hang = - 1)
rect.hclust(hc1, k = 5)
# NbClust函数确定聚类个数，通常会产生两个图，查找第一幅图下降最快的点和第二幅图上升最快的点，可基本确定聚类个数
install.packages("NbClust")
library(NbClust)
NbClust(nutrient.scaled, distance = "euclidean", method = "average")
# par函数用于设置行列，可参考：https://blog.csdn.net/qq_43448491/article/details/104159079
# cutree函数则是对hclust函数的聚类结果进行剪枝。即选择输出指定类别数的系谱聚类结果。剪枝就是说输出指定指定类别数的系谱聚类结果。
# 这里设置的簇为5类
# sort函数将分类数据进行排序，data.frame将数据生成数据框
par(mfrow = c(1, 1))
cutree(hc1, k = 5)
data.frame(group = sort(cutree(hc1, k = 5)))

## 方法2：使用马式距离计算层次分类.不同的方法可能会得到不同的分类结果
d.ma <- dist.ma(nutrient)
hc2 <- hclust(d.ma, method = "ward.D2")
plot(hc2, hang = -1)
rect.hclust(hc2, k = 5)
data.frame(group = sort(cutree(hc2, k = 5)))


## 处理混合类型的数据
install.packages("cluster")
library(cluster)
data(flower)
str(flower)
# 查看flower中V7和V8数值型变量的分布情况
summary(flower$V7)
summary(flower$V8)
# 数值型变量的取值范围不同，因此需要标准化
flower[, c("V7", "V8")] <- scale(flower[, c("V7", "V8")])
summary(flower)
# daisy是cluster程序包中用于计算距离矩阵的方法，这里不能用dist函数进行计算，因为数据框中包含因子型变量
dmat <- as.matrix(daisy(flower))
dmat[1:6, 1:6]
# 使用agnes进行层次聚类，这里的第一个参数是距离矩阵，默认值methods是average（类平均法）
flower.cluster <- agnes(dmat, diss = TRUE)
# pltree()在当前图形设备上绘制树形图，通常对agnes()或diana()产生的分层聚类数据进行绘制
pltree(flower.cluster, hang = -1, main = "")
rect.hclust(flower.cluster, k = 3)
cutree(flower.cluster, k = 3)

10.2.3 R型聚类（指标）

# matrix中的nr是nrow的简写，即6行.可参考：https://www.jb51.net/article/207456.htm
R <- matrix(c(1, 0.852, 0.671, 0.099, 0.234, 0.376,
              0.852, 1, 0.636, 0.055, 0.174, 0.321,
              0.671, 0.636, 1, 0.153, 0.233, 0.252,
              0.099, 0.055, 0.153, 1, 0.732, 0.627,
              0.234, 0.174, 0.233, 0.732, 1, 0.676,
              0.376, 0.321, 0.252, 0.627, 0.676, 1), 
            nr = 6,
            dimnames = list(c("身高", "下肢长", "手臂长", "腰围", "胸围", "臀围")))
# as.dist表示计算距离
d <- as.dist(1 - R)
hc <- hclust(d)
cutree(hc, k = 2)
plot(hc, hang = -1)
rect.hclust(hc, k = 2)

10.3 k均值聚类法

# 适用于样本量较大的情形。通过不断的迭代进行归类
# kmeans方法通过随机选择k个样品作为中心点计算距离，再进行分类。k值不同，分类不同，因此可以通过层次聚类法决定个数
set.seed(1234)
km <- kmeans(nutrient.scaled, centers = 5)
km
sort(km$cluster)
install.packages("fpc")
library(fpc)
# fpc包还提供了另一个展示聚类结果的函数plotcluster()。需要注意的是，数据将被投影到不同的簇中。
plotcluster(nutrient.scaled, km$cluster)