作者：李健民
 
在之前的文章中，我们介绍了apriori算法的思想，应用场景，以及先验项、后验项、支持度、置信度、提升度等基本概念。

详见：啤酒与尿布的故事和临床研究有什么关系 ——关联规则的基本概念
 
本文将结合具体数据讲解apriori算法的分析实例和如何在R中实现。

数据说明：近年科学家发现了一种新的流感病毒正在东南亚蔓延，但是人们对这种病毒所知甚少。在显微镜下，人们发现这种病毒可以分成A、B、C、D、E、F这几个亚型，可以导致人体发生symptom1-symptom8这8个常见症状。研究者记录了500个病人的亚型与出现症状，并记录成数据集symptom（请在文章右侧下载数据集及代码）。如果我们进一步探索其中亚型与症状的规律，或许能找到更合适对症治疗的药。数据录入如下图
  
安装和调用相关的R包，本文涉及的包有readxl、ggplot2、arules、arulesViz，可用install.packages()安装

install.packages("readxl")

library(readxl)

#载入数据（可以用绝对路径，也可以用Rstudio中的import Dataset窗口导入）

symptom <- read_excel("C:/ symptom.xlsx")
 
##1.描述统计及可视化
 
#1.1我们想了解C亚型病毒中各种症状的频数表

library(dplyr)#dplyr包是数据清洗常用的包

subC<-filter(symptom, Subtype=="C")#筛选出目标行,共有163个C亚型对象

factor1<-subC[,3:10]#筛选出目标列，第3列至第10列是symptom1-symptom8

factor1.ma<-as.matrix(factor1)#将多列转换成一个矩阵

factor1.ve<-as.vector(t(factor1.ma))#将多矩阵转换成一个向量

dat1 =na.omit(factor1.ve)#删除NA值

dat1=as.data.frame(table(dat1))#转换成table后再转换成dataframe格式

dat1<-arrange(dat1, desc(Freq))#按降序排列

整理好的dat1数据如下图所示：
 
下面我们把dat1的表格传入封装好的图像函数中，以下经过封装的函数可以反复运用在不同的数据表格，提高工作效率

#1.2封装绘图函数及可视化处理

#首先安装和加载需要用的ggplot2包

install.packages("ggplot2")

library(ggplot2)

#编写函数让分析更加自动化fun1 <- function(data, xlab, ylab, xname, yname) {

  ggplot(data, aes(xlab, ylab)) + 

    geom_bar(aes(fill = xlab), stat = 'identity') + 

    labs(x = xname, y = yname) +

    geom_text(aes(label = ylab), hjust = 1.5, colour = 'white') +

    coord_flip() +

    theme_minimal() +

    theme(legend.position = 'none')
}
 
data <- head(dat1，8)#dat1为传入的数据，8是指显示的行数

xname <- '症状'

yname <- '次数'

fun1(data, reorder(data$dat1, data$Freq), data$Freq, xname, yname)

同理，如果你想了解其他亚型的症状频数表，也可以采用以上的方法
 
 
从这张直方图我们可以看到C亚型中，出现symptom3的次数最多，那么C亚型可能就与symptom3最相关了吧！且慢，我们再进一步用apriori算法看看是不是这么回事。
 
##下面是重头戏了，开始运用apriori进行数据探索
 
##2.要运用arules包里的apriori函数需要首先转化成tansaction格式；值得注意的是，apriori函数只能应用在transactions格式，不能直接用在dataframe上，这转换步骤其实是非常麻烦的。为了便于理解。我们可以输入以下例子

#建立具有9行的list

MyList<-list(c("I1", "I2", "I5"),

c("I2", "I4"),

c("I2", "I3"),

c("I1", "I2", "I4"),

c("I1", "I3"),

c("I2", "I3"),

c("I1", "I3"),

c("I1", "I2", "I3", "I5"),

c("I1", "I2", "I3"))

#重命名这个list

names(MyList)<-paste("Tr",c(1:9),sep="")

#把list转换成transactions格式

MyTrans<-as(MyList,"transactions")

#看看transactions什么样子

inspect(MyTrans)
 
#也可以画图理解，数据小的时候可以这样做，大的话就看不清楚了

image(trans)
 
如图，纵坐标对应9行列表，横坐标的每个阴影对于列表中的元素。我们可以把原数据中的dataframe格式先转为list再转为transactions格式，也可以用如下办法：
 
# 安装和加载R包

install.packages("arules")

install.packages("arulesViz")

library(arules)

library(arulesViz)

x<-symptom[,2:10]

x.ma<-as.matrix(x)

x.ve<-as.vector(t(x.ma))

ID<-rep(1:500, each = 9, len = 4500)

new.df<-data.frame(ID,x.ve)

trans <- as(split(new.df$x.ve,new.df$ID),"transactions")#转换成transactions格式
 
##3.探索规则,我们尝试用apriori算法看看C亚型和什么症状最相关

rules1 <- apriori(data = trans, 

                                  parameter = list(support = 0.08, confidence = 0.3,maxlen = 4),

                                  appearance = list(lhs = "C", default = "rhs"))

#support、confidence、maxlen分别指支持度、置信度与最大项，这几个参数都是根据数据的体质而个人定制的，support设置太高的话，可能一些支持度低但置信度高的规则就会遗漏。support和confidence我们可以先从较低数值开始，逐步筛选规则。maxlen指的是先验项和后验项最大总共能出血几个项目，本次例子我们想看看数据中的C亚组有没有和多个症状联系，设置maxlen=4，可惜并没有看到一对多的规则联系。

rules1 <- sort(rules1, decreasing = TRUE, by = "lift")

inspect(rules1)
 
筛选出lift大于1的三条规则

rules_lift<-subset(x=rules1,subset=slot(object=rules1,name="quality")$lift>1)

inspect(rules_lift)
 
总共生成了13条规则，其中有3条规则的提升度是大于1。其中虽然symptom4的支持度和置信度都不是最高，但是提升度却是最高的，这是因为symptom4在C亚型中出现最多，而在其他亚型比较少见。那为什么symptom3的提升度小于1，是个无效关联呢，这是因为symptom3不仅常出现在亚型C中，在其他亚型也很常见，不具有特殊意义。
 
#把规则进行可视化

plot(rules_lift, method = "graph",engine = "interactive") #plot提升度大于1的规则
 

在亚型（C）与症状（symptom4,6,8）之间的圆圈代表关联程度：颜色越深，代表lift越高，体积越大代表support越大，所以C对symptom4的提升度最大，证明两者有较强的关联，也就是说symptom4在亚型C中出现的比较多；C和symptom6的支持度最大，说明symptom6在亚型C中出现的频率比较大，但关联强度不如symptom4，或许symptom6在其他亚型也不少见。
 
小结：个人经验感觉apriori算法对原始数据格式有一定要求，使用前需要转换成transaction格式。对于庞大且错综复杂的数据集，apriori算法能发挥更好的效果，本文用于练习的数据集案例数和维度都较少，不足以体现算法的威力。apriori算法是一种数据探索的方法，对于样本量没有明确要求，并不像RCT和队列研究那样有章可循。arules包中有自带的数据集Groceries，里面有9835行和169列数据，对于这种级别的数据，运用传统统计学处理显然是不太合适的。对于较小的数据集我们可以从较低的support开始，探索可能存在的关联规则，说不定会有令人惊喜的收获。好好理解其中的原理，结合自身的临床背景，能为我们的研究提供新的思路。
 
 
大家可以从文章右侧，下载文章中用到的Code，供下载练习

扫码关注“医咖会”公众号，及时获取最新重磅研究！