用R语言进行关联分析关联是两个或多个变量取值之间存在的一类重要的可被发现的某种规律性。关联分析目的是寻找给定
数据记录集中数据项之间隐藏的关联关系，描述数据之间的密切度。几个基本概念1. 项集这是一个集合的概念，在一篮子商品中的一件消费品即为一项（Item），则若干项的集合为项集，如{啤酒，尿布}构成一个二元项集。2. 关联规则一般记为的形式，X为先决条件，Y为相应的关联结果，用于表示数据内隐含的关联性。如：，表示购买了尿布的消费者往往也会购买啤酒。关联性强度如何，由三个概念——支持度、置信度、提升度来控制和评价。例：有10000个消费者购买了商品，其中购买尿布1000个，购买啤酒2000个，购买面包500个，同时购买尿布和面包800个，同时购买尿布和面包100个。3. 支持度（Support）支持度是指在所有项集中{X, Y}出现的可能性，即项集中同时含有X和Y的概率：该指标作为建立强关联规则的第一个门槛，衡量了所考察关联规则在“量”上的多少。通过设定最小阈值（minsup），剔除“出镜率”较低的无意义规则，保留出现较为
频繁的项集所隐含的规则。设定最小阈值为5%，由于{尿布，啤酒}的支持度为800/10000=8%，满足基本输了要求，成为频繁项集，保留规则；而{尿布，面包}的支持度为100/10000=1%，被剔除。4. 置信度（Confidence）置信度表示在先决条件X发生的条件下，关联结果Y发生的概率：这是生成强关联规则的第二个门槛，衡量了所考察的关联规则在“质”上的可靠性。相似的，我们需要对置信度设定最小阈值（mincon）来实现进一步筛选。具体的，当设定置信度的最小阈值为70%时，置信度为800/1000=80%，而的置信度为800/2000=40%，被剔除。5. 提升度（lift）提升度表示在含有X的条件下同时含有Y的可能性与没有X这个条件下项集中含有Y的可能性之比：该指标与置信度同样衡量规则的可靠性，可以看作是置信度的一种互补指标。R中Apriori算法算法步骤：1. 选出满足支持度最小阈值的所有项集，即频繁项集；2. 从频繁项集中找出满足最小置信度的所有规则。>library(arules) #加载arules包>click_detail =read.transactions("click_detail.txt",format="basket",sep=",",cols=c(1)) #读取txt文档（文档编码为ANSI）>rules <- apriori(click_detail, parameter =list(supp=0.01,conf=0.5,target="rules")) #调用apriori算法>rulesset of419 rules>inspect(rules[1:10]) #查看前十条规则解释1)library(arules)：加载程序包arules，当然如果你前面没有下载过这个包，就要先install.packages(arules)2)click_detail =read.transactions("click_detail.txt",format="basket",sep=",",cols=c(1))：读入数据read.transactions(file, format =c("basket", "single"), sep = NULL,cols = NULL, rm.duplicates =FALSE, encoding = "unknown")file:文件名，对应click_detail中的“click_detail.txt”format:文件格式，可以有两种，分别为“basket”,“single”，click_detail.txt中用的是basket。basket: basket就是篮子，一个顾客买的东西都放到同一个篮子，所有顾客的transactions就是一个个篮子的组合结果。如下形式，每条交易都是独立的。文件形式：item1,item2item1item2,item3读入后：items 1 {item1,item2}2 {item1}3 {item2,item3}single: single的意思，顾名思义，就是单独的交易，简单说，交易记录为：顾客1买了产品1， 顾客1买了产品2，顾客2买了产品3……(产品1，产品2，产品3中可以是单个产品，也可以是多个产品)，如下形式：trans1 item1trans2 item1trans2 item2读入后：items transactionID1 {item1}trans12 {item1, item2}trans2sep:文件中数据是怎么被分隔的，默认为空格，click_detail里面用逗号分隔cols:对basket, col=1,表示第一列是数据的transaction ids(交易号)，如果col=NULL，则表示数据里面没有交易号这一列；对single，col=c(1,2)表示第一列是transaction ids，第二列是item idsrm.duplicates:是否移除重复项，默认为FALSEencoding:写到这里研究了encoding是什么意思，发现前面txt可以不是”ANSI”类型，如果TXT是“UTF-8”，写encoding=”UTF-8”，就OK了.3)rules <- apriori(click_detail,parameter = list(supp=0.01,conf=0.5,target="rules"))：apriori函数apriori(data, parameter = NULL, appearance = NULL, control = NULL)data:数据parameter：设置参数，默认情况下parameter=list(supp=0.1,conf=0.8,maxlen=10,minlen=1,target=”rules”)supp:支持度（support）conf:置信度（confidence）maxlen,minlen:每个项集所含项数的最大最小值target:“rules”或“frequent itemsets”（输出关联规则/频繁项集）apperence:对先决条件X（lhs），关联结果Y（rhs）中具体包含哪些项进行限制，如：设置lhs=beer，将仅输出lhs含有beer这一项的关联规则。默认情况下，所有项都将无限制出现。control：控制函数性能，如可以设定对项集进行升序sort=1或降序sort=-1排序，是否向使用者报告进程（verbose=F/T）补充通过支持度控制：rules.sorted_sup = sort(rules, by=”support”)通过置信度控制：rules.sorted_con = sort(rules, by=”confidence”)通过提升度控制：rules.sorted_lift = sort(rules, by=”lift”)Apriori算法两步法：1. 频繁项集的产生：找出所有满足最小支持度阈值的项集，称为频繁项集；2. 规则的产生：对于每一个频繁项集l，找出其中所有的非空子集；然后，对于每一个这样的子集a，如果support(l)与support(a)的比值大于最小可信度，则存在规则a==>(l-a)。频繁项集产生所需要的计算开销远大于规则产生所需的计算开销频繁项集的产生几个概念：1， 一个包含K个项的数据集，可能产生2^k个候选集 2，先验原理：如果一个项集是频繁的，则它的所有子集也是频繁的（理解了频繁项集的意义，这句话很容易理解的）；相反，如果一个项集是非频繁的，则它所有子集也一定是非频繁的。 3基于支持度（SUPPORT）度量的一个关键性质：一个项集的支持度不会超过它的子集的支持度（很好理解，支持度是共同发生的概率，假设项集{A,B,C}，{A,B}是它的一个自己，A,B,C同时发生的概率肯定不会超过A,B同时发生的概率）。上面这条规则就是Apriori中使用到的，如下图，当寻找频繁项集时，从上往下扫描，当遇到一个项集是非频繁项集（该项集支持度小于Minsup），那么它下面的项集肯定就是非频繁项集，这一部分就剪枝掉了。一个例子（百度到的一个PPT上的）：当我在理解频繁项集的意义时，在R上简单的复现了这个例子，这里采用了eclat算法，跟apriori应该差不多：代码：item <- list(c("bread","milk"),c("bread","diaper","beer","eggs"),c("milk","diaper","beer","coke"),c("bread","milk","diaper","beer"),c("bread","milk","diaper","coke"))names(item) <- paste("tr",c(1:5),sep = "")itemtrans <- as(item,"transactions") #将List转为transactions型rules = eclat(trans,parameter = list(supp = 0.6,target ="frequent itemsets"),control = list(sort=1))inspect(rules) #查看频繁项集运行后结果：>inspect(rules)items support1{beer, diaper}0.62{diaper, milk} 0.63{bread,diaper}0.64{bread,milk} 0.65{beer} 0.66{milk} 0.87{bread} 0.88{diaper} 0.8以上就是该例子的所有频繁项集，然后我发现少了{bread,milk,diaper}这个项集，回到例子一看，这个项集实际上只出现了两次，所以是没有这个项集的。规则的产生每个频繁k项集能产生最多2k-2个关联规则将项集Y划分成两个非空的子集X和Y-X，使得X ->Y-X满足置信度阈值定理：如果规则X->Y-X不满足置信度阈值，则X’->Y-X’的规则一定也不满足置信度阈值，其中X’是X的子集Apriori按下图进行逐层计算，当发现一个不满足置信度的项集后，该项集所有子集的规则都可以剪枝掉了。# 一、R基本操作# 1、将数据文件mydata1.txt按照以下要求整理成标准形式。#(1)读入数据文件mydata.txt命名为insurance。insurance<-read.table("mydata1.txt")head(insurance)dim(insurance)#192个数据#(2)将insurance转换为3列的矩阵。insurance<-matrix(insurance$V1,nrow = 64,ncol = 3)#nrow =192/3=64insurance#(3)将insurance转换为数据框。insurance<-as.data.frame(insurance)class(insurance)#(4)将列名命名为"District", "Holders"和"Claims"。names(insurance)<-c("District", "Holders","Claims")insurance#(5)随机无放回抽取50行数据。sub<-insurance[sample(1:nrow(insurance),50),]#无放回不用设置replacesub#(6)将抽样数据写入result1.txt。write.table(sub,"result1.txt",row.names = FALSE)####################################################################### 2、将数据文件mydata2.txt按照以下要求整理成标准形式。#(1)读入数据文件mydata2.txt命名为iris。iris<-read.table("mydata2.txt")head(iris)dim(iris)#600个数据#(2)将iris转换为4列的矩阵。iris<-matrix(iris$V1,nrow = 150,ncol = 4)#nrow =600/3=150iris#(3)将iris转换为数据框。iris<-as.data.frame(iris)class(iris)#(4)将列名命名为"Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width"。names(iris)<-c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width")iris#(5)随机无放回抽取90行数据。sub<-iris[sample(1:nrow(iris),90),]#无放回不用设置replacesub#(6)将抽样数据写入result2.txt。write.table(sub,"result2.txt",row.names = FALSE)####################################################################### 3.将数据文件data.csv按照以下要求进行数据预处理。#(1)读入数据文件data.csv命名为nhanes2。nhanes2<-read.csv("data.csv")#(2) 载入缺失值处理所需要的包。install.packages("lattice")install.packages("MASS")install.packages("nnet")library(lattice)library(MASS)library(nnet)#(3) 判断nhanes2是否存在缺失值。sum(is.na(nhanes2))#(4) 利用插补法处理chl变量的缺失值。sub=which(is.na(nhanes2[,4])==TRUE)#在数据集中chl变量是第4列，返回nhanes2数据集中第4列为NA的行dataTR<-nhanes2[-sub,]#将第4列不为NA的数存入数据集dataTRdataTE<-nhanes2[sub,]#将第4列为NA的数存入数据集dataTE中dataTE[,4]<-sample(dataTR[,4],length(dataTE[,4]),replace = T)#在非缺失值中简单抽样dataTE #(5) 将插补法处理后的数据写入result3.txt。write.table(dataTE,"result3.txt",row.names = FALSE)#############################################################################################################################################二、函数调用#1、测得某班学术X（身高(cm)）与Y（体重(kg)）的数据如下,试画出散点图,建立线性回归方程，并作进一步分析。# (1) 建立数据集，并画出散点图，考察数据点的分布趋势，看是否呈直线条状分布。x1<-c(171,175,159,155,152,158,154,164,168,166,159,164)#身高y1<-c(57,64,41,38,35,44,41,51,57,49,47,46)#体重#构建数据集model <- data.frame(x1,y1) #探索性分析-做散点图查看数据的分布情况：plot(x1,y1)# (2)进行回归分析,列出回归方程，画拟合线，并对结果进行解读。# 用lm()函数构建
模型lm.reg<-lm(y1~ x1)# 添加回归曲线查看拟合效果 abline(lm.reg) #模型解读summary(lm.reg)# (3)对回归系数进行假设检验。anova(lm.reg) # 回归模型的方差分析summary(lm.reg) #回归系数t检验：提取模型计算结果，其中有t检验的结果# (4)对回归模型进行诊断。#模型检验对方程进行进一步检验，以检查回归方程是否满足模型的先验条件及模型的稳健性。par(mfrow=c(2,2))#画布分面plot(lm.reg)#结果解读：#1.左上图：残差与拟合图，理论上散点应该散乱的分布在横线两侧；#2.右上图：正太Q-Q图，用于检验因变量的正太分布性，若服从正太分布，则散点应分布在一条直线线#3.左下图：齐方差检验，若满足其方差，则散点在水平线周围随机分布#4.右下图：独立性检验，即一个样本是否会影响另一个样本##################################################################2、研究某抗心律失常药对电刺激狗右心室致颤阙的影响,实验测得狗静脉注射不同剂量的抗心律失常药与右心室致颤阙的数据如下,试画出散点图,建立线性回归方程，并作进一步分析。# (1) 建立数据集，并画出散点图，考察数据点的分布趋势，看是否呈直线条状分布。x <- c(1,3,5,7,9)y <- c(8.03, 14.97, 19.23, 27.83, 36.23)#构建数据集model <- data.frame(x,y) #探索性分析-做散点图查看数据的分布情况：plot(model)#画散点图# (2)进行回归分析,列出回归方程，画拟合线，并对结果进行解读。# 用lm()函数构建模型fm <- lm(y ~ x)#建立回归模型fm# 添加回归曲线查看拟合效果abline(fm)# 添加回归曲线至散点图 #模型解读summary(fm)# (3)对回归系数进行假设检验。anova(fm) # 回归模型的方差分析summary(fm) # 提取模型计算结果，其中有t检验的结果# (4)对回归模型进行诊断。#模型检验对方程进行进一步检验，以检查回归方程是否满足模型的先验条件及模型的稳健性。par(mfrow=c(2,2))#画布分面plot(fm)#结果解读：#1.左上图：残差与拟合图，理论上散点应该散乱的分布在横线两侧；#2.右上图：正太Q-Q图，用于检验因变量的正太分布性，若服从正太分布，则散点应分布在一条直线线#3.左下图：齐方差检验，若满足其方差，则散点在水平线周围随机分布#4.右下图：独立性检验，即一个样本是否会影响另一个样本################################################################### 3、countries数据集含有69个国家和地区的出生率与死亡率。# (1) 请使用K-均值聚类将样本点聚为3个类别。countries=read.csv("countries.csv")head(countries)#查看前6行names(countries)=c("country","birth","death")#修改变量名称var=as.character(countries$country)#将变量country转为字符型并赋值给varfor(i in 1:69) row.names(countries)[i]=var[i]#将数据集的行名命名为国家名称km1=kmeans(countries[,-1],center=3)#用kmeans算法对countries数据集进行聚类# (2) 输出聚类结果及各类别的中心点坐标。km1$cluster#获取类别km1$centers#获取中心点坐标# (3) 绘制聚类结果将中心点以星号标识。#画出聚为四类的类别图，标注中心点。plot(countries[,-1],pch=c(1,2,3))#将中心点用星号标示出来points(km1$centers,pch=8,col="red")#对中心点添加标注legend(km1$centers[1,1],km1$centers[1,2],"Center_1",bty="n",xjust=0.5,cex=0.8)legend(km1$centers[2,1],km1$centers[2,2],"Center_2",bty="n",xjust=0.5,cex=0.8)legend(km1$centers[3,1],km1$centers[3,2],"Center_3",bty="n",xjust=0.5,cex=0.8)# (4) 判断与中国大陆同属于一个类别的国家和地区有哪些。cluster_CHINA=km1$cluster[which(countries$country=="CHINA")]which(km1$cluster==cluster_CHINA)###############################################################################################################################三、数据分析# 1、使用arules软件包中的Groceries数据集，该数据集是某一食品杂货店一个月的真实交易数据，使用R完成以下要求：(软件包：arules；数据集：Groceries； 函数：apriori())# (1)利用apriori()函数进行关联分析，支持度为0.01，置信度为0.5。install.packages("arules")library(arules)data("Groceries")rules0<-apriori(Groceries,parameter=list(support=0.01,confidence=0.5))inspect(rules0[1:10])# (2)利用sort()函数按照支持度排序。rules.sorted_sup<-sort(rules0,by="support")inspect(rules.sorted_sup[1:5])# (3)捆绑销售：寻找蛋黄酱（mayonnaise）的捆绑商品。（supp=0.001,conf=0.1,minlen=2, maxlen=6）rules1=apriori(Groceries,parameter=list(minlen=2,maxlen=6,supp=0.001,conf=0.1),appearance=list(rhs="mayonnaise",default="lhs"))inspect(rules1)# (4)查看销量最高的商品。itemsets_apr=apriori(Groceries,parameter=list(supp=0.001,target="frequent itemsets"),control=list(sort=-1))inspect(itemsets_apr[1:5])# (5)适合捆绑销售的商品。（supp=0.001,minlen=2, maxlen=3）itemsets_apr1=eclat(Groceries,parameter=list(supp=0.001,minlen=2,maxlen=3,target="frequent itemsets"),control=list(sort=-1))inspect(itemsets_apr1[1:5])# (6)关联规则的可视化（support=0.001,con=0.5）install.packages("arulesViz")library(arulesViz)rules5=apriori(Groceries,parameter=list(support=0.002,con=0.5))rules5plot(rules5)######################################################################## 2、根据breast-cancer-wisconsin.csv威斯康星州乳腺癌数据集，通过对数据的分析，提取出关键特征来判断乳腺癌患病情况。(软件包：rpart；函数：rpart()。)# (1)属性名依次设置为"编号","肿块厚度","肿块大小","肿块形状","边缘黏附","单个表皮细胞大小","细胞核大小","染色质","细胞核常规","有丝分裂","类别"),并将类别为2的设为"良性",为4的设为"恶性"。install.packages("rpart")library(rpart)install.packages("rpart.plot")library(rpart.plot)#############加载数据breast.cancer<-read.csv('breast-cancer-wisconsin.csv',header=F)head(breast.cancer)#数据整理names(breast.cancer)=c("编号","肿块厚度","肿块大小","肿块形状","边缘黏附","单个表皮细胞大小","细胞核大小","染色质","细胞核常规","有丝分裂","类别")breast.cancer$类别[breast.cancer$类别==2]="良性"breast.cancer$类别[breast.cancer$类别==4]="恶性"head(breast.cancer)# (2)抽取训练数据集为原数据的70%，测试数据集取30%。#数据预处理(分层抽样,划分训练集和测试集)#分别计算良性和恶性组中应抽取测试集样本数，记为a，ba=round(0.3*sum(breast.cancer$类别=="良性"))b=round(0.3*sum(breast.cancer$类别=="恶性"))ab #输出a,b值install.packages("sampling")library(sampling)#使用strata函数对数据集中的“分组油耗”变量进行分层抽样sub=strata(breast.cancer,stratanames="类别",size=c(b,a),method="srswor")sub #所抽出的所有测试集样本信息#生成训练集train1和测试集test1train1=breast.cancer[-sub$ID_unit,]test1=breast.cancer[sub$ID_unit,]nrow(train1)nrow(test1) #显示训练集和测试集的行数，检查两者比例是否为7：3# (3) minsplit=5，建立决策树。#CART建立分类树formula_cla=类别~肿块厚度+肿块大小+肿块形状+边缘黏附+单个表皮细胞大小+细胞核大小+染色质+细胞核常规+有丝分裂cla1=rpart(formula_cla,train1,method="class",minsplit=5)#cla1# (4)选择cp=0.05来剪枝。######修改cp的值cla2=rpart(formula_cla,train1,method="class",minsplit=5,cp=0.05)cla2# (5)画出type为2和4的树图。rpart.plot(cla1,type=2)#修改typerpart.plot(cla1,type=4) # (6)测试数据进行预测，并输出混淆矩阵，给出模型准确率为。#预测pre1=predict(cla1,test1,type="class")pre1table(test1$类别,pre1)#获取混淆矩阵#计算样本错误率error1<-sum(as.numeric(pre1!=test1$类别))/nrow(test1)error1#################################################################### 3、美国科罗拉多州某加油站连续 57 天的OVERSHORTS序列“OVERSHORTS.csv”# (1) 判断该序列的平稳性与纯随机性。# （时序图检验、白噪声检验）install.packages("fUnitRoots")install.packages("TSA")install.packages("forecast")install.packages("zoo")library(fUnitRoots)library(TSA)library(forecast)library(zoo)#读取数据c<-read.csv("OVERSHORTS.csv")#转换为时间序列overshort<-ts(c$overshort,start = 1)#平稳性，纯随机（白噪声检验）## 绘制序列的时间序列图plot.ts(overshort, xlab = "time", ylab = "prop")##对序列做单位根检验unitrootTest(overshort)##对序列做白噪声检验Box.test(overshort, lag = 1, type = "Ljung-Box")# (2) 如果序列平稳且非白噪声，选择适当模型拟合该序列的发展。(10分)# （模型的识别、参数估计（模型显著性、模型参数的显著性））#模型识别##观察自相关，偏自相关图，模型定阶par(mfrow=c(1,2))acf(overshort)###衰减到零是突然的，所以自相关系数1阶截尾pacf(overshort)### 衰减到零不是突然的，所以偏相关系数托尾# 推荐模型为 MA(1)##或者对序列进行模型识别，自动定阶auto.arima(overshort)# 推荐模型为 MA(1)#参数估计###模型检验x.fit<-arima(overshort,order=c(0,0,1),method="ML")x.fit##对残差x.fit$residual进行白噪声检验for(i in 1:2) print(Box.test(x.fit$residual,lag=6*i))##P>0.05,接受原假设，即残差为白噪声，所以拟合模型显著有效####参数检验###模型参数的显著性检验t1<--0.8477/0.1206pt(t1,df=56,lower.tail=T) ###p<0.05参数显著非零t0<--4.7942/1.0253pt(t0,df=56,lower.tail=T) ###p<0.05参数显著非零# (3) 利用拟合模型，预测该加油站未来5天的OVERSHORTS。(10分)# （模型预测、绘制预测图）####模型预测c<-read.csv("OVERSHORTS.csv")x<-ts(c$overshort,start=1)x.fit<-arima(x,order=c(0,0,1))x.fitx.fore<-forecast(x.fit,h=5)#预测x.foreplot(x.fore)###############################################################4、使用是survival软件包中的“pbc”数据集，该数据集记录的是肝硬化数据， 使用R完成一下要求：(软件包：survival；数据集：pbc； 函数：Surv()、survfit()、survdiff()、coxph()、cox.zph(), 将答案保存在“姓名.doc”文件中。)# (1)生成生存分析对象，拟合生存曲线模型。install.packages("survival") #安装survival包library(survival) #加载survival包#使用survival包自带的“pbc”数据集为例（418*20） data("pbc")str(pbc)head(pbc)#生成生存分析对象Sur_Obj<-Surv(pbc$time,pbc$status)Sur_Obj#拟合曲线模型model<-survfit(Sur_Obj~1) summary(model)# (2)两种方法绘制生存曲线。plot(model,ylab = "生存率",xlab="天")#用survminer进行漂亮的展示install.packages("survminer")library(survminer) ggsurvplot(model, data = pbc)# (3)进行单因素比较分析，并进行结果解释。#survdiff(formula)函数进行log-rank检验。survdiff(Sur_Obj~pbc$trt) #trt是分组条件# (4)考虑年龄，性别以及trt是否会影响肝硬化的生存时间,进行多因素分析Cox模型的建立，并进行结果解释。coxmodel<-coxph(Sur_Obj~pbc$age+pbc$sex+pbc$bili)coxmodel# (5)模型诊断——PH检验。zphmodel<-cox.zph(coxmodel)zphmodel############################################################### 5、life.csv为50位急性淋巴细胞白血病病人的数据，包括：入院治疗时取得外辕血中细胞数X1，淋巴结浸润等级X2，出院后有无巩固治疗X3(1表示有巩固治疗，0表示无巩固治疗)；随访后，变量Y=0表示生存期在1年以内，Y=1表示生存时间在1年以上，使用R完成一下要求：(函数：glm()，predict()。)# (1)建立全变量logistic回归，对模型结果进行解释。life<-read.csv("life.csv")#建立全变量logistic回归glm.sol<-glm(Y~X1+X2+X3, family=binomial, data=life)#回归模型解读summary(glm.sol)# (2)预测当X1=5,X2=2,X3=0时，y的概率是多少？pre<-predict(glm.sol, data.frame(X1=5,X2=2,X3=0))p<-exp(pre)/(1+exp(pre))p# (3)预测当X1=5,X2=2,X3=1时，y的概率是多少？(6分)pre<-predict(glm.sol, data.frame(X1=5,X2=2,X3=1))p<-exp(pre)/(1+exp(pre))p# (4)对回归模型参数进行检验，用step()函数做变量筛选。step(glm.sol)glm.new<-glm(Y~X2+X3, family=binomial, data=life)summary(glm.new)# (5)对筛选后的变量进行建模，预测。pre<-predict(glm.new, data.frame(X2=2,X3=0))p<-exp(pre)/(1+exp(pre))ppre<-predict(glm.new, data.frame(X2=2,X3=1))p<-exp(pre)/(1+exp(pre))p
