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

Chapter 11

聚类分析和判别分析都是样本分类的统计方法。
其中主要的区别是，聚类分析是未知分类，无监督学习；判别分析是已知分类，有监督学习。

11.1 距离判别

data(iris)
# cor用于计算相关系数矩阵
cor(iris[, 1:4])
# colMeans是计算列的均值，rowMeans是计算行的均值
m.setosa <- colMeans(iris[1:50, 1:4])
m.setosa
m.versicolor <- colMeans(iris[51:100, 1:4])
m.versicolor
m.virginica <- colMeans(iris[101:150, 1:4])
m.virginica

install.packages("biotools")
library(biotools)
# boxM根据一个分类因子对多元正态数据的协方差矩阵齐性进行Box’s m检验。检验是基于卡方近似。
# boxM(data, grouping)
# iris[, -5]表示数据的前五列，iris[, 5]表示数据的第六列
boxM(iris[, -5], iris[, 5])

# cov用于计算协方差矩阵
v.setosa <- cov(iris[1:50, 1:4])
v.setosa
v.versicolor <- cov(iris[51:100, 1:4])
v.versicolor
v.virginica <- cov(iris[101:150, 1:4])
v.virginica

# 计算马氏距离mahalanobis(x, center, cov, reverse = FALSE，…)
# x表示有p列的数据的向量或矩阵，center表示长度为p的分布或第二数据向量的均值。
# iris[1, 1:4]表示iris数组中第一个数据
mahalanobis(iris[1, 1:4], m.setosa, v.setosa)
mahalanobis(iris[1, 1:4], m.versicolor, v.versicolor)
mahalanobis(iris[1, 1:4], m.virginica, v.virginica)

# iris[, 1:4]第一个参数为空，表示所有行
# data.frame是将所有数据形成有个数据框，便于展示
d.setosa <- mahalanobis(iris[,1:4], m.setosa, v.setosa)
d.versicolor <- mahalanobis(iris[,1:4], m.versicolor, v.versicolor)
d.virginica <- mahalanobis(iris[,1:4], m.virginica, v.virginica)
d <- data.frame(d.setosa, d.versicolor, d.virginica)
head(d)

# apply(d, MARGIN = 1, FUN = which.min)这里MARGIN = 1表示行，MARGIN = 2表示列，
# 后面的FUN = which.min表示函数是对行数据取最小值
# 参考：https://blog.csdn.net/weixin_48364077/article/details/107748853
index <- apply(d, MARGIN = 1, FUN = which.min)
index
# 使用factor函数转化为因子型变量并添加标签
type <- factor(index, labels = c("setosa", "versicolor", "virginica"))
type
# 使用table函数将测量的数据与真实的数据通过列联表进行判定
table(type, iris$Species)
# which函数是获取具体的位置。
# 这里就是真实数据中iris$Species == "versicolor"和
# 我们测量数据中type == "virginica"两个条件都满足的情况，就是我们混肴的数据
which(type == "virginica" & iris$Species == "versicolor")

11.2 K最邻近判别

# 这里需要训练集和测试集
set.seed(1234)
nrow(iris)
# sample(1:150, 100)函数表示从150个数据中随机抽取100个
s <- sample(1:150, 100)
# iris[s, ]表示随机抽取的s行数据是训练集
train <- iris[s, ]
nrow(train)
# iris[-s, ]表示随机抽取的剩下的s行数据是测试集
test <- iris[-s, ]
nrow(test)
# cl就是训练集中的标签
cl <- train[, 5]
cl
install.packages("class")
library(class)
# knn函数实现KNN分类，knn(train, test, cl, k = 1, l = 0, prob = FALSE, use.all = TRUE)
# 前几个参数分别是训练集，测试集，标签，临近点个数
iris.knn <- knn(train[, -5], test[, -5], cl)
iris.knn
# 将训练后进行分类的数据与真实数据的标签进行对比，生成一个混肴矩阵
confusion.matrix <- table(iris.knn, test[, 5])
confusion.matrix
#diag是获取矩阵对角线的值，对这些值求和除以行数就是正确率
accuracy <- sum(diag(confusion.matrix))/nrow(test)
accuracy
# vector表示生成一个长度为20的向量，初值为0
accuracy <- vector(length = 20)
# 写一个循环函数，计算k值在1-20时的精度
for(i in 1:20) {
  iris.knn <- knn(train[, -5], test[, -5], cl, k = i)
  confusion.matrix <- table(iris.knn, test[, 5])
  accuracy[i] <- sum(diag(confusion.matrix))/nrow(test)
}
accuracy
# 使用plot函数将精度绘制出来，选取何时的k值
plot(accuracy, type = "b", xlab = "k")

11.3 Fisher判别

# lda用于线性判别
install.packages("MASS")
library(MASS)
iris.ld <- lda(Species ~ Sepal.Length + Sepal.Width + 
                 Petal.Length + Petal.Width, data = iris)
iris.ld
iris.pred <- predict(iris.ld)
iris.pred$class
table(iris.pred$class, iris$Species)
which(iris.pred$class == "virginica" & iris$Species == "versicolor")
which(iris.pred$class == "versicolor" & iris$Species == "virginica")

LD1 <- iris.pred$x[, 1]
LD2 <- iris.pred$x[, 2]
col <- as.numeric(iris$Species)
pch <- as.numeric(iris$Species)
plot(LD1, LD2, col = col, pch = pch)
# legend函数添加图例，points是画点，cex是缩放比例。
# 参考：https://blog.csdn.net/qq_27390023/article/details/121188434
legend("top", legend = c("setosa", "versicolor", "virginica"), 
       col = 1:3, pch = 1:3)
points(LD1[c(71, 84)], LD2[c(71, 84)], cex = 2)
points(LD1[134], LD2[134], cex = 2)

11.4 Bayes判别

install.packages("klaR")
library(klaR)
iris.bayes <- NaiveBayes(Species ~ Sepal.Length + Sepal.Width + 
                           Petal.Length + Petal.Width, data = iris)
names(iris.bayes)
plot(iris.bayes)
iris.pred <- predict(iris.bayes)
table(iris.pred$class, iris$Species)