R语言学习笔记之聚类分析
使用k-means聚类所需的包：
factoextra
cluster #加载包
library(factoextra)
library(cluster)l
#数据准备
使用内置的R数据集USArrests
#load the dataset
data("USArrests")
#remove any missing value (i.e, NA values for not available)
#That might be present in the data
USArrests <- na.omit(USArrests)#view the first 6 rows of the data
head(USArrests, n=6)
在此数据集中，列是变量，行是观测值
在聚类之前我们可以先进行一些必要的数据检查即数据描述性统计，如平均值、标准差等
desc_stats <- data.frame( Min=apply(USArrests, 2, min),#minimum
Med=apply(USArrests, 2, median),#median
Mean=apply(USArrests, 2, mean),#mean
SD=apply(USArrests, 2, sd),#Standard deviation
Max=apply(USArrests, 2, max)#maximum
)
desc_stats <- round(desc_stats, 1)#保留小数点后一位head(desc_stats)
变量有很大的方差及均值时需进行标准化
df <- scale(USArrests)
#数据集群性评估
使用get_clust_tendency()计算Hopkins统计量
res <- get_clust_tendency(df, 40, graph = TRUE)
res$hopkins_stat
## [1] 0.3440875
#Visualize the dissimilarity matrix
res$plot
Hopkins统计量的值<0.5，表明数据是高度可聚合的。另外，从图中也可以看出数据可聚合。
#估计聚合簇数
由于k均值聚类需要指定要生成的聚类数量，因此我们将使用函数clusGap()来计算用于估计最优聚类数。函数fviz_gap_stat()用于可视化。
set.seed(123)
## Compute the gap statistic
gap_stat <- clusGap(df, FUN = kmeans, nstart = 25, K.max = 10, B = 500)
# Plot the result
fviz_gap_stat(gap_stat)
图中显示最佳为聚成四类（k=4）
#进行聚类
set.seed(123)
km.res <- kmeans(df, 4, nstart = 25)
head(km.res$cluster, 20)
# Visualize clusters using factoextra
fviz_cluster(km.res, USArrests)
#检查cluster silhouette图
Recall that the silhouette measures (SiSi) how similar an object ii is to the the other objects in its own cluster versus those in the neighbor cluster. SiSi values range from 1 to - 1:
A value of SiSi close to 1 indicates that the object is well clustered. In the other words, the object ii is similar to the other objects in its group.
A value of SiSi close to -1 indicates that the object is poorly clustered, and that assignment to some other cluster would probably improve the overall results.
sil <- silhouette(km.res$cluster, dist(df))
rownames(sil) <- rownames(USArrests)
head(sil[, 1:3])
#Visualize
fviz_silhouette(sil)
图中可以看出有负值，可以通过函数silhouette()确定是哪个观测值
neg_sil_index <- which(sil[, "sil_width"] <0)
sil[neg_sil_index, , drop = FALSE]
## cluster neighbor sil_width
## Missouri 3 2 -0.07318144
#eclust():增强的聚类分析
与其他聚类分析包相比，eclust()有以下优点：
简化了聚类分析的工作流程
可以用于计算层次聚类和分区聚类
eclust()自动计算最佳聚类簇数。
自动提供Silhouette plot
可以结合ggplot2绘制优美的图形
#使用eclust()的K均值聚类
# Compute k-means
res.km <- eclust(df, "kmeans")
# Gap statistic plot
fviz_gap_stat(res.km$gap_stat)
# Silhouette plotfviz_silhouette(res.km)
## cluster size ave.sil.width
## 1 1 13 0.31
## 2 2 29 0.38
## 3 3 8 0.39
#使用eclust（）的层次聚类
# Enhanced hierarchical clustering
res.hc <- eclust(df, "hclust") # compute hclust
fviz_dend(res.hc, rect = TRUE) # dendrogam
#下面的R代码生成Silhouette plot和分层聚类散点图。
fviz_silhouette(res.hc) # silhouette plot
## cluster size ave.sil.width
## 1 1 19 0.26
## 2 2 19 0.28
## 3 3 12 0.43
fviz_cluster(res.hc) # scatter plot
#Infos
This analysis has been performed using R software (R version 3.3.2)
在r中看函数源代码：在R中，代码可以分为如下几个级别：
首先，是你输入了函数对象名称，你可以直接看到代码的，如要获得函数对象fivenum的代码，就只需要在Console中键入函数对象名称fivenum就可以得到如下结果：
function (x, na.rm = TRUE)
{
xna <- is.na(x)
if (na.rm)
x <- x[!xna]
else if (any(xna))
return(rep.int(NA, 5))
x <- sort(x)
n <- length(x)
if (n == 0)
rep.int(NA, 5)
else {
n4 <- floor((n + 3)/2)/2
d <- c(1, n4, (n + 1)/2, n + 1 - n4, n)
0.5 * (x[floor(d)] + x[ceiling(d)])
}
}
<environment: namespace:stats>
从上面的例子可以看出，这类函数对象的代码是最容易看到的，也是我们学习的最好的材料了，而R中最大多数的函数对象是以这种方式出现的。
其次，我们在输入mean这类函数名次的时候，会出现如下结果：
function (x, ...)
UseMethod("mean")
<environment: namespace:base>
这表示函数作者把函数“封”起来了。这个时候我们可以先试一试methods(mean)，利用methods函数看看mean这个函数都有哪些类型的，我们得到的结果如下：
[1] mean.data.frame mean.Date mean.defaultmean.difftime mean.POSIXct mean.POSIXlt
其实对此可以有一个简单的理解，虽然不够精确。因为在R中，mean函数可以求得属于不同类型对象的平均值，而不同类型对象平均值的求法还是有一些小小差 异的，比如说求一个向量的平均值和求一个数据框的平均值就有所差异，就要编写多个mean函数，然后“封”起来，以一个统一的mean出现，方便我们使 用。这正好也反映了R有一种类似泛型编程语言的性质。
既然我们已经知道mean中还有这么多种类，我们可以输入mean.default试一试就可以得到：
function (x, trim = 0, na.rm = FALSE, ...)
{
if (!is.numeric(x) &&!is.complex(x) &&!is.logical(x)) {
warning("argument is not numeric or logical: returning NA")
return(as.numeric(NA))
}
if (na.rm)
x <- x[!is.na(x)]
trim <- trim[1]
n <- length(x)
if (trim >0 &&n >0) {
if (is.complex(x))
stop("trimmed means are not defined for complex data")
if (trim >= 0.5)
return(stats::median(x, na.rm = FALSE))
lo <- floor(n * trim) + 1
hi <- n + 1 - lo
x <- sort.int(x, partial = unique(c(lo, hi)))[lo:hi]
n <- hi - lo + 1
}
.Internal(mean(x))
}
<environment: namespace:base>
同样就可以得到mean.data.frame、mean.Date、mean.difftime、mean.POSIXct、mean.POSIXlt 的具体内容了。值得注意的是，在R中，出现有多个同样近似功能的函数封装为一个函数的时候（这时候在函数中多半会出现类似UseMethod函数使用的情 况），我们不妨先输入mean.default试一试。这种形式的函数在R中一般作为默认的函数表示。
第三，这是一种特殊的情况，有人认为应该和第二种是一类，但是我还是要提出来单独归类。在这种情况也和第二种的原因有些类似，但并不是完全一致。
也许我们大家都很熟悉plot函数了吧，输入函数名plot的时候，我们会得到如下结果：
function (x, y, ...)
{
if (is.null(attr(x, "class")) &&is.function(x)) {
nms <- names(list(...))
if (missing(y))
y <- {
if (!"from" %in% nms)
0
else if (!"to" %in% nms)
1
else if (!"xlim" %in% nms)
NULL
}
if ("ylab" %in% nms)
plot.function(x, y, ...)
else plot.function(x, y, ylab = paste(deparse(substitute(x)),
"(x)"), ...)
}
else UseMethod("plot")
}
<environment: namespace:graphics>
请注意plot函数中也出现了UseMethod这个函数，但是和mean不同的是，前面有相当多的语句用于处理其他一些事情。这个时候，我们也使用methods(plot)来看看，得到如下结果：
plot.acf* plot.data.frame*plot.Date* plot.decomposed.ts* plot.default
plot.dendrogram*plot.densityplot.ecdf plot.factor*plot.formula*
plot.hclust*plot.histogram* plot.HoltWinters* plot.isoreg*plot.lm
plot.medpolish* plot.mlmplot.POSIXct* plot.POSIXlt* plot.ppr*
plot.prcomp*plot.princomp* plot.profile.nls* plot.spec plot.spec.coherency
plot.spec.phase plot.stepfunplot.stl* plot.table* plot.ts
plot.tskernel* plot.TukeyHSD
不看不知道，一看吓一跳，还以为我们输入plot的输出就是函数本身，结果也许不是如此。可能有人已经理解了，其实最后的UseMethod函数实在默认的调用plot.default函数，赶快再看看plot.default函数吧，发现它再调用plot.xy函数，再看看plot.xy函数，再plot.xy函数中调用了一个.Internal(plot.xy(xy, type, pch, lty, col, bg, cex, lwd, ...))函数，也许这就是真正起作用的函数了吧。思路基本上就是如此了，是否这个时候您可以获得一些阅读查找R函数内容的乐趣。
除了直接输入FUN.default形式外，还可以使用getS3method(FUN,"default")来获得代码。这样就解决了绝大多数函数代码查看的工作了。
在第二种情况种，我们说了一般可以通过FUN.default获得想要的结果。但是只有称为generic的函数才有这种“特权”。而lm等则没有，不过我们也可以尝试使用methods(lm)来看看结果如何，发现：
[1] lm.fit lm.fit.null lm.influence lm.wfit lm.wfit.null
Warning message:
function 'lm' appears not to be generic in: methods(lm)
出现了警告信息，表示说lm不是泛型函数，但是还是给出了结果lm.fit等，大致上可以看成是和lm相关的系列函数吧。这样子就出现了有趣的局面，比如说既有plot.ts，也有ts.plot。
依照第三种情况，我们发现竟然有的函数用星号标识了的，比如plot.stl*等，当我们输入plot.stl，甚至是plot.stl*的时候都会给出 要么找不到这个对象，要么干脆是代码错误的信息。原来凡是用了*标识的函数，都是隐藏起来的函数，估计是怕被人看见（其实这是玩笑话）！我们要看这些函数 的代码，我们该怎么办呢？其实也很容易，使用功能强大的getAnywhere(FUN)，看看这个函数的名称，就可以猜想到它的功能估计是很强大的， Anywhere的内容都可以找到！getAnywhere(plot.stl)的结果如下：
A single object matching 'plot.stl' was found
It was found in the following places
registered S3 method for plot from namespace stats
namespace:stats
with value
function (x, labels = colnames(X), set.pars = list(mar = c(0,
6, 0, 6), oma = c(6, 0, 4, 0), tck = -0.01, mfrow = c(nplot,
1)), main = NULL, range.bars = TRUE, ..., col.range = "light gray")
{
sers <- x$time.series
ncomp <- ncol(sers)
data <- drop(sers %*% rep(1, ncomp))
X <- cbind(data, sers)
colnames(X) <- c("data", colnames(sers))
nplot <- ncomp + 1
if (range.bars)
mx <- min(apply(rx <- apply(X, 2, range), 2, diff))
if (length(set.pars)) {
oldpar <- do.call("par", as.list(names(set.pars)))
on.exit(par(oldpar))
do.call("par", set.pars)
}
for (i in 1:nplot) {
plot(X[, i], type = if (i <nplot)
"l"
else "h", xlab = "", ylab = "", axes = FALSE, ...)
if (range.bars) {
dx <- 1/64 * diff(ux <- par("usr")[1:2])
y <- mean(rx[, i])
rect(ux[2] - dx, y + mx/2, ux[2] - 0.4 * dx, y -
mx/2, col = col.range, xpd = TRUE)
}
if (i == 1 &&!is.null(main))
title(main, line = 2, outer = par("oma")[3] >0)
if (i == nplot)
abline(h = 0)
box()
right <- i%%2 == 0
axis(2, labels = !right)
axis(4, labels = right)
axis(1, labels = i == nplot)
mtext(labels[i], side = 2, 3)
}
mtext("time", side = 1, line = 3)
invisible()
}
<environment: namespace:stats>
注意到前面有一段解释型的语言，描述了我们要找的这个函数放在了什么地方等等。其实对任意我们可以在R中使用的函数，都可以先试一试getAnywhere，看看都有些什么内容。算是一个比较“霸道”的函数。
在上面plot.xy函数中，我们还可以看到.Internal这个函数，类似的也许还可以看到.Primitive、.External、.Call等函数这就和R系统内部工作方式和与外部接口的定义有关了，如果对这些函数有兴趣的话，就要学习组成R系统的源代码了。
最后，如果真的想阅读组成R系统本身的源代码，在各个CRAN中均有下载。你可以得到组成R系统所需要的材料。其中很多C语言（还有就是F）的源代码，均 是精心挑选过的算法，哪怕就是想学从头到尾编写具体的算法，也是学习的好材料。同时，你可以看到R系统内部是如何构成的，理解了这些对于高效使用R有至关 重要的作用。这个范畴的材料就要着重看一看R-Lang和R-inits了。
至此，R中阅读代码的内容就依照我的理解介绍了一下。随后将有一些R代码示例的分析注解、语言本身、R应用的和行业使用的材料翻译和具体例子说明。欢迎大家多多和我交流，一起进步。
一行一行来。basic.stats <- function(x,more=F) { # 建立名叫basic.stats的函数，参数为x和more，more默认是F就是不用输入，但你也可以输入，有额外效果。
stats <- list() #建立名叫stats的列表类型变量
clean.x <- x[!is.na(x)] #把x中的NA全部踢掉，留下有用的数据记为clean.x
stats$mean <- mean(clean.x)# 计算clean.x的均值 赋给列表中的mean单元
stats$std <- sd(clean.x)# 计算clean.x的标准差 赋给列表中的std单元
stats$med <- median(clean.x) # 计算clean.x的中位数 赋给列表中的med单元
if(more) { #如果你在函数中输入2个变量，默认是basic.stats(x)，你可以输入basic.stats(x, y) 有额外效果
stats$skew <- sum(((clean.x-stats$mean)/stats$std)^3)/length(clean.x) #计算偏度 赋给列表中的skew单元
stats$kurt <- sum(((clean.x-stats$mean)/stats$std)^4)/length(clean.x) - 3 #计算峰度 赋给列表中的kurt单元
unlist(stats) #最后拆解列表变量stats 使其变为简单的向量数值变量
