R语言入门级实例——用igragh包分析社群

引入——

　　本文的主要目的是初步实现R的igraph包的基础功能，包括绘制关系网络图(social relationship)、利用算法进行社群发现(community detecting)。对于R语言零基础的同学非常友好。以下R代码中如有含义不清的，建议尝试先在R编辑器中输入?xxx()进行查询（xxx是函数或语句名）。此外，StackOverflow论坛也帮博主小白看懂了不少报错信息。
　　主要参考资料为《R语言与网站分析》[李明著][机械工业出版社][2014.04] 的9.3节《关系网络分析》。

0.背景

　　现已获得超市中商品的名称、分类以及大量顾客购物篮子中的商品信息，任务是分析哪些商品存在相关性，经常被放在一起购买。题外话，这种分析的一例经典应用就是沃尔玛超市的“啤酒与尿布”，感兴趣者可自行搜索或参见Jocelyn_燕的一篇博客.

1.原始数据及初步处理

　　数据来源是Kaggle竞赛的数据库instacart-market-basket-analysis.下载压缩文件之后，将有用的数据合并到一个Excel文件中，此处需要order_product，order，products，departments的数据.注意，这个文件极大，order_product_prior这个spread sheet里的数据在Excel里已经无法完全显示，博主就截取了前500条信息，形成了mini数据集，以下对数据集的操作都是针对这个mini表进行的.如下：

 　　为了达到参考书上的数据形式，需要先整理这个Excel，形成如下图只有四列数据的形式.这里博主不太熟悉R的操作，就用Python的循环处理了，代码可附在文章最后.

　　这是当初处理数据集的一些文件，由于不会用R完成所有命令，显得很笨拙hhh.

2.数据集导入

　　导入的数据集包含四列，原商品编号过大，不便于处理，p_id、d_id分别是商品、商品分类的新编号，如下图：（这些也是用Python代劳的）

3.建立关系网络与绘图

步骤描述：

引用igraph包，建立空关系网络并设置点数据→
为点数据添加商品号以及商品分类属性→
添加线数据→
plot出来发现是非连通图（存在孤立的点的图），有两个未连通的点（点43,点44）,只用手动对点的个数减2即可
将点的个数修改后，重新跑前面的所有代码即可

这部分代码如下：（完整代码见文末）

 1 #建立空关系网络并设置点数据
 2 library(igraph)
 3 gdata<-graph.empty(directed=F)
 4 #num<-ncol(cart)
 5 num<-ncol(cart)-2 #修改点的个数
 6 gdata<-add.vertices(gdata,num)
 7 
 8 #为点数据添加商品号以及商品分类属性
 9 category<-c();item<-c()
10 for(i in colnames(cart))
11 {
12   if(i!=136&& i!=140)
13   {
14     category<-c(category,data$d_id[which(data$p_id==i)[1]] )
15     item<-c(item,data$p_id[which(data$p_id==i)[1]] ) 
16   }
17 }
18 V(gdata)$category<-category
19 V(gdata)$item<-item
20 
21 #添加线数据
22 #依次遍历每个订单，读取每个订单内的商品ID，并存放于向量item.i
23 for(i in 1:nrow(cart))
24 { item.i<-c()
25 for(j in 1:ncol(cart))
26 {
27   if(cart[i,j]==1)
28   {
29     item.i<-cbind(item.i,colnames(cart)[j])
30   }
31 }
32 #建立向量内不同商品间的关联联系
33 item.i.num<-length(item.i)
34 from<-c();to<-c()
35 for(m in 1:(item.i.num-1))
36 {
37   from<-c(from,item.i[-c((item.i.num-m+1):item.i.num)])
38   to<-c(to,item.i[-c(1:m)])
39 }
40 if(i>1)
41 {
42   edges<-rbind(edges,matrix(c(from,to),nc=2))
43 }
44 else
45 {
46   edges<-matrix(data=c(from,to),nc=2)
47 }
48 }
49 edges0<-edges
50 labels<-union(unique(edges[,1]), unique(edges[,2]))
51 ids<-1: length(labels)#对点的编号重新编码，因为在igraph中边信息的ids必须连续
52 names(ids)<-labels
53 newfrom<-as.character(edges[,1]);newto<-as.character(edges[,2])
54 edges<-matrix (c(ids[newfrom],ids[newto]), nc=2)
55 
56 #添加线信息并设置线权重
57 gdata<-add.edges(gdata,t(edges[-1,]))#t()是矩阵转置函数
58 E(gdata)$weight<-count.multiple(gdata)
59 gdata<-simplify(gdata, remove.multiple=TRUE, remove.loops = TRUE, edge.attr.comb = 'mean')
60 #最后一个参数一定是edge.attr.comb,不是edges.attr.comb
61 dev.off()#关闭图形设备
62 plot(gdata,edge.width=E(gdata)$weight,main="gdata", edge.label=E(gdata)$weight)
63 
64 #发现是非连通图，有两个未连通的点（点43,点44）,只用手动对点的个数减2即可
65 #将点的个数修改后，重新跑前面的所有代码

View Code

　　画出来的效果如下：

4.社群发现与绘图

　　此处采用自旋玻璃法(spinglass community detecting)进行社群发现。其他社群发现的方法包括中心势、标签传播、随机游走等，这几种方法在算法效率与模拟方式上其实存在不同点。但限于篇幅，此处不再介绍。对这几种方法感兴趣者可自行搜索或参考以下论文（引用格式不够规范，但应该能搜索到）：

[1]J¨org Reichardt & Stefan Bornholdt (2008) Statistical Mechanics of Community Detection <=spinglass相关

[2]M. Girvan & M. E. J. Newman (2001) Community structure in social and biological networks <=中心势betweeness相关

[3]Jierui Xie & Boleslaw K. Szymanski (2013) LabelRank: A Stabilized Label Propagation Algorithm for Community Detection in Networks <=标签传播labelrank相关

[4]Pascal Pons and Matthieu Latapy (2006) Computing Communities in Large Networks Using Random Walks <=随机游走randomwalk相关

总之，在这里spinglass方法适用于购物车商品分析。
　　另外，需要注意：
　　①社群发现必须基于连通图（即，所有点上都在线上，没有孤立的点）；
　　②此处的社群个数对应之后画子图的分组个数。

步骤描述：

对不同商品类别的点配置不同颜色→
建立绘图分组member.list，作为plot函数mark.groups参数的列表对象→
画图并手动添加图例→
可添加点的标签属性vertex.label，呈现原有编号

这部分代码如下：

 1 ##社群发现并绘制关系图（自旋玻璃法）
 2 member<-spinglass.community(gdata, weights= E(gdata)$weight)
 3 V(gdata)$member<-member$membership
 4 member.num<-length(table(V(gdata)$member)); member.num #注意：此处的社群个数对应之后的绘图分组
 5 
 6 #对不同商品类别的点配置不同颜色
 7 mem.col<-rainbow(length(unique(V(gdata)$category)),alpha=0.5)#注意设置alpha值调节对比度
 8 V(gdata)$color<-mem.col[V(gdata)$category]
 9 #建立设置绘图分组(plot函数的mark.groups参数)的列表对象member.list
10 member.list<-list()
11 for(i in 1:member.num)
12 {
13   member.list<-c(member.list, list(which(V(gdata)$member==i)))
14 }
15 #svg(filename=paste(root, "demol.svg",sep=""), width = 14, height = 14)
16 #画图并手动添加图例
17 legend0<-c("dairy eggs","produce","meat seafood","beverages","pantry","bakery","frozen","snacks")
18 #plot(gdata, vertex.size=10, layout=layout.fruchterman.reingold, vertex.color=V(gdata)$color, edge.width=scale(E(gdata)$weight, center=F)+1, mark.groups=member.list)
19 plot(gdata, vertex.label=V(gdata)$item, vertex.size=10, layout=layout.fruchterman.reingold, vertex.color=V(gdata)$color, edge.width=scale(E(gdata)$weight, center=F)+1, mark.groups=member.list)
20 #第二个plot加了label属性
21 legend("topleft",legend=legend0, pch=16, col=mem.col, bty="n", cex=1)

View Code

　　画出图如下（右图为加了lable标签后的效果，所有点恢复了真实编号，而不是左图中临时的连续编号）：

　　OK! 看上去还不错。

　　现在我们得到的图里，每个点的颜色对应左侧图例中的不同商品分类（蛋奶制品、烘焙类、冷冻品、零食等等），点与点之间的连线代表两个曾在同一购物篮子(即订单信息order)中出现过。现在利用算法已经发现了五个可能存在的社群，即，在这个图中关系更密切的点的集合，由浅色“冲积扇”形状色块标出。右图中，点的编号就是原mini数据库中的商品号码。现在就可以研究能不能得出有趣的结论了！

　　对照如下图的数据库，上方右图中编号81,80,31,119的商品位于一个社群中。也许数据量再大些能说明热爱有机蔬果的顾客也偏好矿泉水？

5.绘制子图

　　为了单独研究形成的各个社群，还可以把关系图拆成子图分别绘制。

　　有两种方法画子图：
　　A.设置par，用循环一次性画出；
　　B.依次画每个图，放大后更清晰

 1 #绘制不同社群内的关系图
 2 #svg(filename=paste(root, "demol.svg",sep=""), width = 14, height = 14)
 3 #par(mfcol=c(3,2))
 4 for(i in 1:length(table(member$membership)))
 5 {
 6   tmp.g<-induced.subgraph(gdata,which(V(gdata)$member==i));V(tmp.g)
 7   member.list<-list()
 8   tmp.category<-as.numeric(names(table(V(tmp.g)$category)))
 9   for(j in tmp.category)
10   {
11     member.list<-c(member.list,list(which(V(tmp.g)$category==j)))
12   }
13   plot(tmp.g, vertex.size=10,layout=layout.fruchterman.reingold, edge.width=scale(E(tmp.g)$weight,center=F)+1,mark.groups=member.list,vertex.label=V(tmp.g)$item)
14   #手动添加图例
15   #legend("topleft",legend= ,pch=16,col=mem.col,bty="n",cex=1)
16 }

View Code

　　子图如下：

                            

                        

————————————————分割线—————————————————-

6.完整代码

 1 ls()
 2 rm(list = ls())
 3 #初步读取数据
 4 root="C:/Users/asus/Desktop/"
 5 data<-read.csv(file=paste(root,"购物车.csv",sep=""),header=T,encoding="UTF-8");
 6 colname1<-colnames(data)
 7 colname1[1]<-"order_id"
 8 colnames(data)<-colname1
 9 #由于预先对数据集进行了处理，此处不需要书上分离商品名、类别并编号的步骤
10 ##建立关系网络
11 #用cast函数转化格式
12 #重铸函数cast(md,formula,FUN),其中md是已融和的数据，formula描述了想要的结果，
13 #而FUN是数据整合函数，例如mean，也可自定义多值整合函数。默认为统计函数。
14 
15 #install.packages('reshape')
16 library('reshape')
17 data<-cbind(data,value=1)
18 #cast返回数据框，再转换成矩阵
19 cart=as.matrix(cast(data,order_id~p_id,value="value",fill=0))
20 cart[,-1]<-ifelse(cart[,-1]>=1,1,0)#好像有点多余，因为此数据集中每个购物篮子中的某件商品只被记了一次
21 
22 #注：这是最开始的数据准备部分，限于篇幅，后面的部分就是前文各小节代码的拼凑综合，不再重复复制粘贴。

View Code

参考资料：《R语言与网站分析》[李明著][机械工业出版社][2014.04] 的9.3节《关系网络分析》。

R代码部分引用自原书作者，增加了注释，结合R语言语法的变化也有改动。

7.数据预处理部分的Python代码（可以用R的指令代替）

1.对商品重新编号（商品分类的重新编号类似，此处不赘述）

 1 import openpyxl
 2 import re
 3  
 4 def Exceldivide(file_dir):
 5  wb=openpyxl.load_workbook(file_dir)
 6  sheet=wb.get_sheet_by_name('prior_order')
 7  tuple(sheet['A1':'E507'])
 8  t=1
 9  for i in range(2,508):
10      fd=False
11      for j in range(2,i):
12         if sheet.cell(row=i, column=4).value==sheet.cell(row=j, column=4).value:
13             sheet.cell(row=i, column=6).value=sheet.cell(row=j, column=6).value
14             fd=True
15      if fd==False:
16          sheet.cell(row=i, column=6).value=t
17          t+=1
18  return wb
19  
20 g=Exceldivide('C:\Users\asus\Desktop\购物篮子简化版.xlsx')
21 g.save('C:\Users\asus\Desktop\购物篮子简化版.xlsx')

View Code

2.保留被重复购买过的商品（这一步在数据集较大时可省去）

 1 import openpyxl
 2 import re
 3  
 4 def Exceldivide(file_dir):
 5  wb=openpyxl.load_workbook(file_dir)
 6  sheet0=wb.get_sheet_by_name('prior_order')  #
 7  sheet1=wb.get_sheet_by_name('repeat')
 8  sheet2=wb.get_sheet_by_name('order') 
 9  tuple(sheet0['A1':'F507'])
10  tuple(sheet1['A1':'B45'])
11  tuple(sheet2['A1':'D45'])
12  i=1
13  for rows in sheet0['F2':'F507']:
14       for cell0 in rows:
15          for rows2 in sheet1['A2':'A45']:
16              for cell1 in rows2:
17                  if cell0.value==cell1.value:
18                       i+=1
19                       sheet2.cell(row=i, column=1).value=sheet0.cell(row=int(cell0.coordinate[1:]), column=1).value
20                       sheet2.cell(row=i, column=2).value=sheet0.cell(row=int(cell0.coordinate[1:]), column=2).value
21                       sheet2.cell(row=i, column=3).value=sheet0.cell(row=int(cell0.coordinate[1:]), column=3).value
22                       sheet2.cell(row=i, column=4).value=sheet0.cell(row=int(cell0.coordinate[1:]), column=6).value
23  return wb
24  
25 g=Exceldivide('C:\Users\asus\Desktop\购物篮子简化版.xlsx')
26 g.save('C:\Users\asus\Desktop\购物篮子简化版.xlsx')

View Code

小注：写作本文源于博主小白去年一段做RA的经历，当时与队友们共同学习社会网络分析(Social Network Analysis,SNA)，主要参考书是上文提及的《R语言与网站分析》9.3节。博主小白与搭档负责实现书上的两个实例，但由于教材没有提供数据来源、R语言语法近几年的变化，中间费了一番波折，故写作本文，主要内容为博主负责的“购物篮子商品相关性分析”实例，转载请注明来源。如有疏漏，还望指正！

用代码改变世界！就是这样，喵！