R语言中实现的包：Rtsne
来学习下核心函数 Rtsne() 的主要参数：
下面以常规模型构建过程中产生的risk和riskScore数据为例：
下面结合ggplot2进行绘图：
另一种降维方法：PCA，已在之前的笔记中写过。 【R>>PCA】主成分分析
参考链接：
t-SNE一种高效的降维算法
如何对混合型数据做聚类分析利用聚类分析，我们可以很容易地看清数据集中样本的分布情况。以往介绍聚类分析的文章中通常只介绍如何处理连续型变量，这些文字并没有过多地介绍如何处理混合型数据（如同时包含连续型变量、名义型变量和顺序型变量的数据）。本文将利用 Gower 距离、PAM（partitioning around medoids）算法和轮廓系数来介绍如何对混合型数据做聚类分析。
R语言
本文主要分为三个部分：
距离计算
聚类算法的选择
聚类个数的选择
为了介绍方便，本文直接使用 ISLR 包中的 College 数据集。该数据集包含了自 1995 年以来美国大学的 777 条数据，其中主要有以下几个变量：
连续型变量
录取率
学费
新生数量
分类型变量
公立或私立院校
是否为高水平院校，即所有新生中毕业于排名前 10% 高中的新生数量占比是否大于 50%
本文中涉及到的R包有：
In [3]:
set.seed(1680) # 设置随机种子，使得本文结果具有可重现性
library(dplyr)
library(ISLR)
library(cluster)
library(Rtsne)
library(ggplot2)
Attaching package: ‘dplyr’
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
构建聚类模型之前，我们需要做一些数据清洗工作：
录取率等于录取人数除以总申请人数
判断某个学校是否为高水平院校，需要根据该学校的所有新生中毕业于排名前 10% 高中的新生数量占比是否大于 50% 来决定
In [5]:
college_clean <- College %>%
mutate(name = row.names(.),
accept_rate = Accept/Apps,
isElite = cut(Top10perc,
breaks = c(0, 50, 100),
labels = c("Not Elite", "Elite"),
include.lowest = TRUE)) %>%
mutate(isElite = factor(isElite)) %>%
select(name, accept_rate, Outstate, Enroll,
Grad.Rate, Private, isElite)
glimpse(college_clean)
Observations: 777
Variables: 7
$ name(chr) "Abilene Christian University", "Adelphi University", "...
$ accept_rate (dbl) 0.7421687, 0.8801464, 0.7682073, 0.8369305, 0.7564767, ...
$ Outstate(dbl) 7440, 12280, 11250, 12960, 7560, 13500, 13290, 13868, 1...
$ Enroll (dbl) 721, 512, 336, 137, 55, 158, 103, 489, 227, 172, 472, 4...
$ Grad.Rate (dbl) 60, 56, 54, 59, 15, 55, 63, 73, 80, 52, 73, 76, 74, 68,...
$ Private (fctr) Yes, Yes, Yes, Yes, Yes, Yes, Yes, Yes, Yes, Yes, Yes,...
$ isElite (fctr) Not Elite, Not Elite, Not Elite, Elite, Not Elite, Not...
距离计算
聚类分析的第一步是定义样本之间距离的度量方法，最常用的距离度量方法是欧式距离。然而欧氏距离只适用于连续型变量，所以本文将采用另外一种距离度量方法—— Gower 距离。
Gower 距离
Gower 距离的定义非常简单。首先每个类型的变量都有特殊的距离度量方法，而且该方法会将变量标准化到[0,1]之间。接下来，利用加权线性组合的方法来计算最终的距离矩阵。不同类型变量的计算方法如下所示：
连续型变量：利用归一化的曼哈顿距离
顺序型变量：首先将变量按顺序排列，然后利用经过特殊调整的曼哈顿距离
名义型变量：首先将包含 k 个类别的变量转换成 k 个 0-1 变量，然后利用 Dice 系数做进一步的计算
优点：通俗易懂且计算方便
缺点：非常容易受无标准化的连续型变量异常值影响，所以数据转换过程必不可少；该方法需要耗费较大的内存
利用 daisy 函数，我们只需要一行代码就可以计算出 Gower 距离。需要注意的是，由于新生入学人数是右偏变量，我们需要对其做对数转换。daisy 函数内置了对数转换的功能，你可以调用帮助文档来获取更多的参数说明。
In [6]:
# Remove college name before clustering
gower_dist <- daisy(college_clean[, -1],
metric = "gower",
type = list(logratio = 3))
# Check attributes to ensure the correct methods are being used
# (I = interval, N = nominal)
# Note that despite logratio being called,
# the type remains coded as "I"
summary(gower_dist)
Out[6]:
301476 dissimilarities, summarized :
Min. 1st Qu.Median Mean 3rd Qu. Max.
0.0018601 0.1034400 0.2358700 0.2314500 0.3271400 0.7773500
Metric : mixed Types = I, I, I, I, N, N
Number of objects : 777
此外，我们可以通过观察最相似和最不相似的样本来判断该度量方法的合理性。本案例中，圣托马斯大学和约翰卡罗尔大学最相似，而俄克拉荷马科技和艺术大学和哈佛大学差异最大。
In [7]:
gower_mat <- as.matrix(gower_dist)
# Output most similar pair
college_clean[
which(gower_mat == min(gower_mat[gower_mat != min(gower_mat)]),
arr.ind = TRUE)[1, ], ]
Out[7]:
In [8]:
# Output most dissimilar pair
college_clean[
which(gower_mat == max(gower_mat[gower_mat != max(gower_mat)]),
arr.ind = TRUE)[1, ], ]
Out[8]:
聚类算法的选择
现在我们已经计算好样本间的距离矩阵，接下来需要选择一个合适的聚类算法，本文采用 PAM（partioniong around medoids）算法来构建模型：
PAM 算法的主要步骤：
随机选择 k 个数据点，并将其设为簇中心点
遍历所有样本点，并将样本点归入最近的簇中
对每个簇而言，找出与簇内其他点距离之和最小的点，并将其设为新的簇中心点
重复第2步，直到收敛
该算法和 K-means 算法非常相似。事实上，除了中心点的计算方法不同外，其他步骤都完全一致 。
优点：简单易懂且不易受异常值所影响
缺点：算法时间复杂度为 O(n2)O(n2)
聚类个数的选择
我们将利用轮廓系数来确定最佳的聚类个数，轮廓系数是一个用于衡量聚类离散度的内部指标，该指标的取值范围是[-1,1]，其数值越大越好。通过比较不同聚类个数下轮廓系数的大小，我们可以看出当聚类个数为 3 时，聚类效果最好。
In [9]:
# Calculate silhouette width for many k using PAM
sil_width <- c(NA)
for(i in 2:10){
pam_fit <- pam(gower_dist,
diss = TRUE,
k = i)
sil_width[i] <- pam_fit$silinfo$avg.width
}
# Plot sihouette width (higher is better)
plot(1:10, sil_width,
xlab = "Number of clusters",
ylab = "Silhouette Width")
lines(1:10, sil_width)
聚类结果解释
描述统计量
聚类完毕后，我们可以调用 summary 函数来查看每个簇的汇总信息。从这些汇总信息中我们可以看出：簇1主要是中等学费且学生规模较小的私立非顶尖院校，簇2主要是高收费、低录取率且高毕业率的私立顶尖院校，而簇3则是低学费、低毕业率且学生规模较大的公立非顶尖院校。
In [18]:
pam_fit <- pam(gower_dist, diss = TRUE, k = 3)
pam_results <- college_clean %>%
dplyr::select(-name) %>%
mutate(cluster = pam_fit$clustering) %>%
group_by(cluster) %>%
do(the_summary = summary(.))
print(pam_results$the_summary)
[[1]]
accept_rateOutstate Enroll Grad.Rate Private
Min. :0.3283 Min. : 2340 Min. : 35.0 Min. : 15.00 No : 0
1st Qu.:0.7225 1st Qu.: 8842 1st Qu.: 194.8 1st Qu.: 56.00 Yes:500
Median :0.8004 Median :10905 Median : 308.0 Median : 67.50
Mean :0.7820 Mean :11200 Mean : 418.6 Mean : 66.97
3rd Qu.:0.8581 3rd Qu.:13240 3rd Qu.: 484.8 3rd Qu.: 78.25
Max. :1.0000 Max. :21700 Max. :4615.0 Max. :118.00
isElite cluster
Not Elite:500 Min. :1
Elite: 0 1st Qu.:1
Median :1
Mean :1
3rd Qu.:1
Max. :1
[[2]]
accept_rateOutstate Enroll Grad.Rate Private
Min. :0.1545 Min. : 5224 Min. : 137.0 Min. : 54.00 No : 4
1st Qu.:0.4135 1st Qu.:13850 1st Qu.: 391.0 1st Qu.: 77.00 Yes:65
Median :0.5329 Median :17238 Median : 601.0 Median : 89.00
Mean :0.5392 Mean :16225 Mean : 882.5 Mean : 84.78
3rd Qu.:0.6988 3rd Qu.:18590 3rd Qu.:1191.0 3rd Qu.: 94.00
Max. :0.9605 Max. :20100 Max. :4893.0 Max. :100.00
isElite cluster
Not Elite: 0 Min. :2
Elite:69 1st Qu.:2
Median :2
Mean :2
3rd Qu.:2
Max. :2
[[3]]
accept_rateOutstate Enroll Grad.Rate Private
Min. :0.3746 Min. : 2580 Min. : 153 Min. : 10.00 No :208
1st Qu.:0.6423 1st Qu.: 5295 1st Qu.: 694 1st Qu.: 46.00 Yes: 0
Median :0.7458 Median : 6598 Median :1302 Median : 54.50
Mean :0.7315 Mean : 6698 Mean :1615 Mean : 55.42
3rd Qu.:0.8368 3rd Qu.: 7748 3rd Qu.:2184 3rd Qu.: 65.00
Max. :1.0000 Max. :15516 Max. :6392 Max. :100.00
isElite cluster
Not Elite:199 Min. :3
Elite: 9 1st Qu.:3
Median :3
Mean :3
3rd Qu.:3
Max. :3
PAM 算法的另一个优点是各个簇的中心点是实际的样本点。从聚类结果中我们可以看出，圣弗朗西斯大学是簇1 的中心点，巴朗德学院是簇2 的中心点，而密歇根州州立大学河谷大学是簇3 的中心点。
In [19]:
college_clean[pam_fit$medoids, ]
Out[19]:
可视化方法
t-SNE 是一种降维方法，它可以在保留聚类结构的前提下，将多维信息压缩到二维或三维空间中。借助t-SNE我们可以将 PAM 算法的聚类结果绘制出来，有趣的是私立顶尖院校和公立非顶尖院校这两个簇中间存在一个小聚类簇。
In [22]:
tsne_obj <- Rtsne(gower_dist, is_distance = TRUE)
tsne_data <- tsne_obj$Y %>%
data.frame() %>%
setNames(c("X", "Y")) %>%
mutate(cluster = factor(pam_fit$clustering),
name = college_clean$name)
ggplot(aes(x = X, y = Y), data = tsne_data) +
geom_point(aes(color = cluster))
进一步探究可以发现，这一小簇主要包含一些竞争力较强的公立院校，比如弗吉尼亚大学和加州大学伯克利分校。虽然无法通过轮廓系数指标来证明多分一类是合理的，但是这 13 所院校的确显著不同于其他三个簇的院校。
In [25]:
tsne_data %>%
filter(X >15 &X <25,
Y >-15 &Y <-10) %>%
left_join(college_clean, by = "name") %>%
collect %>%
.[["name"]]
Out[25]:
‘Kansas State University’
‘North Carolina State University at Raleigh’
‘Pennsylvania State Univ. Main Campus’
‘SUNY at Buffalo’
‘Texas A&M Univ. at College Station’
‘University of Georgia’
‘University of Kansas’
‘University of Maryland at College Park’
‘University of Minnesota Twin Cities’
‘University of Missouri at Columbia’
‘University of Tennessee at Knoxville’
‘University of Texas at Austin’
TSNE降维和UMAP一样，应用场景相似，具体的原理感兴趣的或者数据比较好的可以自行去学习。这里就不多说了，我们直接开始分析及作图吧。TSNE分析我们使用Rtsne这个包。
首先加载R包：
读入表达矩阵数据：
TSNE降维分析很简单，就一句话。
得到的Atsne这个文件是一个list，降维的结果在Atsne$Y中，对其行进行命名，并创建画图文件。
画图和之前UMAP一样，采用ggplot即可：
还可以添加置信椭圆：
除了ggplot，这里我们介绍一下另外一种做法，降维图其实就是散点图，有一个画图神包ggscatter可以实现。
看起来效果还是很好的，还可以再加修饰，感兴趣的小伙伴自行探索。
更进一步，还可以做成单细胞那样，看看某一个基因在不同样本中的表达。
挺不错的！！！
想要示例数据的可以打赏截图联系作者获取，记得留下邮箱！
下节预告---limma包分析转录组芯片数据
