R语言：商业数据分析实例（2）【k-means, t-test】

适用于初学者。内容包括k-mean和t-test的使用。

使用到的数据：链接：https://pan.baidu.com/s/1yhzQSdquizLayXamM0wygg
提取码：3b7i

前言：k-means实现

k-means算法，是一种最广泛使用的聚类算法。k-means以k作为参数，把数据分为k个组，通过迭代计算过程，将各个分组内的所有数据样本的均值作为该类的中心点，使得组内数据具有较高的相似度，而组间的相似度最低。（引用自：张丹(Conan)）

k-means工作原理：

初始化数据，选择k个对象作为中心点。
遍历整个数据集，计算每个点与每个中心点的距离，将它分配给距离中心最近的组。
重新计算每个组的平均值，作为新的聚类中心。
上面2-3步，过程不断重复，直到函数收敛，不再新的分组情况出现。

案例简介

用到的数据共4张表（cvs格式），为2003-2019年美国纽约市房地产交易数据。

NYC_HISTORICAL包含：交易ID，社区ID，地址，建筑类型，时间，价格，面积等；

BOROUGH包含：BOROUGH_ID和BOROUG名称；

BUILDING_CLASS包含：建筑ID和建筑类型等；

NEIGHBORHOOD 包含：街道ID和BOROUGH_ID等；

本案例需要根据ID整合数据，从中提取目标信息。

本章目标

任务1：使用以下KPI执行k-means聚类。描述所选社区所属的聚类。

自2009年以来，住宅物业的价格中值
自2009年以来，住宅物业的销售数量
自2009年以来，住宅物业的价格标准差
自2009年以来，住宅地产的每平方米价格

任务2：选择另一个社区（如 ID：31），并检验以下假设：从2009年开始，社区（ ID：29）的平均住宅物业价格大于所选的另外一个社区

开始：载入所需的包，并设置工作地址

library(lubridate)    #year()
library(tidyverse)    #csv2
library(factoextra)
library(cluster)
setwd("C:/Users/10098/Desktop/AD571/571，A34")

读取所有csv文件。第四个NYC_HISTORICAL较大（150M），且用分号隔开，故用csv2读取

BOROUGH        <- read.csv("BOROUGH.csv", header=TRUE)
BUILDING_CLASS <- read.csv("BUILDING_CLASS.csv", header=TRUE)
NEIGHBORHOOD   <- read.csv("NEIGHBORHOOD.csv", header=TRUE)
NYC_HISTORICAL <- read_csv2("NYC_HISTORICAL.csv")

将各表数据整合

NYC_HISTORICAL <- mutate(NYC_HISTORICAL, Y = year(SALE_DATE))
NEIGHBORHOOD = left_join(NEIGHBORHOOD,BOROUGH,  by = 'BOROUGH_ID')
df <- NYC_HISTORICAL %>%left_join(BUILDING_CLASS, by = c('BUILDING_CLASS_FINAL_ROLL'='BUILDING_CODE_ID')) %>%left_join(NEIGHBORHOOD, by = 'NEIGHBORHOOD_ID')  %>%select(NEIGHBORHOOD_ID, SALE_DATE,SALE_PRICE,GROSS_SQUARE_FEET,TYPE,SALE_ID,Y) %>%filter(TYPE == 'RESIDENTIAL')  %>%group_by(NEIGHBORHOOD_ID)   %>% subset(Y>=2009)  %>%subset(GROSS_SQUARE_FEET != 0) %>% na.omit()

取出所需KPI ，并重新命名一下列名

median_sale_price  <- summarize(df, median(SALE_PRICE))                        #KPI1
number_of_sales    <- count(df, NEIGHBORHOOD_ID)                               #KPI2
standard_deviation <- summarize(df, sd(SALE_PRICE))                            #KPI3
AVEPRICE_price     <- summarize(df, sum(SALE_PRICE) / sum(GROSS_SQUARE_FEET))  #KPI4df_2 <- median_sale_price   %>%full_join(number_of_sales,    by = 'NEIGHBORHOOD_ID')   %>%full_join(standard_deviation, by = 'NEIGHBORHOOD_ID')   %>%full_join(AVEPRICE_price,     by = 'NEIGHBORHOOD_ID')   colnames(df_2) <- c('ID','median', 'amount', 'sd', 'AVEPRICE_price')

做kmeans之前先做scale(），它使数据框的列居中和缩放（也就是每个值减去均值除标准差）

为什么要scale(）？

KNN、K-means和SVM 等距离算法受特征范围的影响大。这些算法使用数据点之间的距离来确定相似性。如果不缩放，变量具有不同的尺度，算法可能对具有更高量级的变量赋予更高的权重，进而影响算法的性能

scale_df <- scale(df_2[,-1]) # 此处scale时不用包括第一列‘ID’

开始kmeans，因为会随机取点，可以设置seed保证每次结果不变。

设定不同的center数量做了4次（因为我们不知道分成几个聚类好）

set.seed(123)
k2 <- kmeans(scale_df, center = 2, nstart = 25)
k3 <- kmeans(scale_df, center = 3, nstart = 25)
k4 <- kmeans(scale_df, center = 4, nstart = 25)
k5 <- kmeans(scale_df, center = 5, nstart = 25)

利用fviz_cluster画图看一下（fviz_cluster：基于ggplot2的分区方法的优雅可视化，包括kmeans）

p2 <- fviz_cluster(k2, data = scale_df) + ggtitle("k=2")
p3 <- fviz_cluster(k3, data = scale_df) + ggtitle("k=3")
p4 <- fviz_cluster(k4, data = scale_df) + ggtitle("k=4")
p5 <- fviz_cluster(k5, data = scale_df) + ggtitle("k=5")

单独运行（如p2）即可出图，这个为了好看，将4张图放在一起查看

library(gridExtra)
grid.arrange(p2, p3, p4, p5, nrow = 2)    #把4张图放在一张画布上便于查看

效果如下：

可以看到，这里133社区单独成了一个聚类（由于随机，每个人的结果可能不同）。

这是因为k-means使用的是欧根距离，容易受到异常值的影响。可以用进阶的PAM改进。

（PAM是对k-means的一种改进算法，能降低异常值对于聚类效果的影响）【以后介绍】

如何确定center数量？

上面做了4种情况，但如何确定center数量？不幸的是，这个问题没有明确答案。聚类的最佳数量某种程度上是主观的，取决于用于测量相似性的方法和用于划分的参数。（kassambara）

实际操作中，具体情况还是要根据business goal去设定个数。

这里介绍3种方法：

1.肘法（Elbow method）

即使得集群内总变化 [或集群内总平方和 (WSS )] 尽量小。

f1 <- fviz_nbclust(scale_df, kmeans, method = "wss", nstart = 25)

图如下：

我们看到随着k的变大， WSS越来越小，经验上说4比较好（因为在4已经明显下降，之后的比较平缓，也就是k增大带来的效果越来越不明显）。这个方法之所以叫肘法也源于此。

2.平均轮廓法（Average silhouette method）

计算平均轮廓宽度。较高的平均轮廓宽度表示良好的聚类。

f2 <- fviz_nbclust(scale_df, kmeans, method = "silhouette", nstart = 25)

显然该方法建议我们使用2

3.差距统计法（Gap statistic method）

由R. Tibshirani、G. Walther 和 T. Hastie 发表（斯坦福大学，2001 年）。该方法可以应用于任何聚类方法。其将不同 k 值的集群内变化的总和与其在数据的空参考分布下的预期值进行比较。最佳聚类的估计值将是使间隙统计量最大化的值。

f3 <- fviz_nbclust(scale_df, kmeans, nstart = 25,  method = "gap_stat", nboot = 50)

该方法建议我们使用1。（部分是因为之前提到的异常值导致的）

如果是一般情况图如下：

注：如果想同时显示上述三张表，可使用 grid.arrange(f1, f2, f3, nrow = 3)

分析

按照上述3个方法，选择聚类数（这里选了4）。将聚类属性合并到数据集，计算4个KPI在不用聚类下的均值

df_3  <- mutate(df_2[,-1], cluster=k4$cluster)   #（将所属聚类放入数据）cluster_mean <- df_3 %>%group_by(cluster) %>%summarise_all("mean")

t-test

选择另一个社区（如 ID：31），并检验以下假设：从2009年开始，社区（ ID：29）的平均住宅物业价格大于所选的另外一个社区。

先筛选出两个社区的数据（可以自制函数取代下面冗杂的方法）

df_29 <- NYC_HISTORICAL %>%left_join(BUILDING_CLASS, by = c('BUILDING_CLASS_FINAL_ROLL' = 'BUILDING_CODE_ID'))  %>%left_join(NEIGHBORHOOD, by = 'NEIGHBORHOOD_ID')  %>%select(NEIGHBORHOOD_ID, SALE_DATE,SALE_PRICE,GROSS_SQUARE_FEET,TYPE,SALE_ID,Y) %>%filter(NEIGHBORHOOD_ID =='29') %>%filter(TYPE == 'RESIDENTIAL')  %>%subset(Y>=2009)  %>%subset(GROSS_SQUARE_FEET != 0) %>% group_by(Y)   %>%na.omit() df_31 <- NYC_HISTORICAL %>%left_join(BUILDING_CLASS, by = c('BUILDING_CLASS_FINAL_ROLL' = 'BUILDING_CODE_ID'))  %>%left_join(NEIGHBORHOOD, by = 'NEIGHBORHOOD_ID')  %>%select(NEIGHBORHOOD_ID, SALE_DATE,SALE_PRICE,GROSS_SQUARE_FEET,TYPE,SALE_ID,Y) %>%filter(NEIGHBORHOOD_ID =='31') %>%filter(TYPE == 'RESIDENTIAL')  %>%subset(Y>=2009)  %>%subset(GROSS_SQUARE_FEET != 0) %>% group_by(Y)   %>%na.omit()

进行t-test

X <- df_29$SALE_PRICE   # 29
Y <- df_31$SALE_PRICE   # 31
t.test(X,Y,alternative = "greater", mu = 0, paired = FALSE, conf.level = .95)

p = 0.1105 > 0.05, 否定原假设(X>Y为假)，则X>Y为真，即X>Y。

大功告成！

案例源于波士顿大学，感谢教授Tara Kelly

Reference：
http://blog.fens.me/r-cluster-kmeans/
http://www.sthda.com/english/articles/29-cluster-validation-essentials/96-determiningthe-optimal-number-of-clusters-3-must-know-methods/