R语言之决策树和随机森林

Python057

R语言之决策树和随机森林,第1张

R语言之决策树和随机森林

总结决策树之前先总结一下特征的生成和选择,因为决策树就是一种内嵌型的特征选择过程,它的特征选择和算法是融合在一起的,不需要额外的特征选择。

一、特征生成:

特征生成是指在收集数据之时原始数据就具有的数据特征,这些数据特征由收集的数据决定(其实也就是在产品定型时设定的需要收集的数据特征),当然,在数据预处理时,也可以在此基础上构造一些新的数据特征,这些特征越多越好,表示你考虑问题比较周全,具体那些变量有用或没用,这要交给下一步特征选择来决定。

二、特征选择

特征选择是指在原有数据特征的基础上,去除重要性比较低的特征变量,过滤出有用的特征变量。这里比较困难的是搞清楚什么样的特征比较重要?这需要根据具体的问题具体分析,有些变量的选择可以很直观的看出来,但这种直觉也不一定正确。对于常用特征选择方法主要有:过滤型、包装型、内嵌型。

过滤型:是指你可以根据某个统计量的大小排序来选择特征变量,如相关系数、p值、R值等

包装型:是指在一个特征集合中选取最优的特征子集。具体需要考虑:用什么样的算法来选取?选取的最优的标准是什么?

常用的算法是分步回归包括向前搜索、向后删除、双向搜索

向前搜索:每次选取一个能使模型预测或分类效果最好的特征变量进来,进来后不退出,直到模型改善效果不再明显;

向后删除:是指每次从特征全集中每次删除一个特征变量能使模型预测或分类效果最好,退出后不进来,直到模型改善效果不再明显;

双向搜索:是指每次每次删除一个特征变量或加入一个特征变量能使模型预测或分类效果最好,退出的不进来,进来的不退出,直到模型改善效果不再明显;

这里再提一下特征变量选择的几个标准:p值、R值、AIC(越小效果越好)、BIC(越小效果越好)、熵(越小效果越好)

内嵌型:这里应该主要就是像决策树这样的情况,算法内部完成特征变量的选取。

三、决策树

决策的几个要点:1、如何决策?(也就是如何树如何分叉)------熵和信息增益---这里面包含的就是特征的选择?哪个特征变量包含的信息量大,就排在前面,至于最后树的深度就决定特征变量的个数。

当然不同的算法使用的衡量的标准不同,还有:信息增益比、基尼不纯系数

2、如何剪枝?-----一般是事后剪枝

3、连续性变量如何离散化?-----阈值的选择

熵:是指信息的混合程度(混乱程度),熵【0-1】越大表示该集合中混合的信息越多,也就表明这次的分叉效果不好还是有很多不同类的信息混在一起

信息增益:熵值的减少量,越大越好

决策树模型特点:模型易于解释;存储空间较小,以树的形式存储,决策树是一个弱分类器,不能完全分类,需要把多个弱分类器通过多数投票法组合在一起。

四、R包实现决策树

library(rpart)

library(rpart.plot)

## rpart.control对树进行一些设置

## xval是10折交叉验证

## minsplit是最小分支节点数,这里指大于等于20,那么该节点会继续分划下去,否则停止

## minbucket:叶子节点最小样本数

## maxdepth:树的深度

## cp全称为complexity parameter,指某个点的复杂度,对每一步拆分,模型的拟合优度必须提高的程度

ct <- rpart.control(xval=10, minsplit=20, cp=0.1)

## kyphosis是rpart这个包自带的数据集

## na.action:缺失数据的处理办法,默认为删除因变量缺失的观测而保留自变量缺失的观测。

## method:树的末端数据类型选择相应的变量分割方法:

## 连续性method=“anova”,离散型method=“class”,计数型method=“poisson”,生存分析型method=“exp”

## parms用来设置三个参数:先验概率、损失矩阵、分类纯度的度量方法(gini和information)

## cost是损失矩阵,在剪枝的时候,叶子节点的加权误差与父节点的误差进行比较,考虑损失矩阵的时候,从将“减少-误差”调整为“减少-损失”

data("Kyphosis")

fit <- rpart(Kyphosis~Age + Number + Start,data=kyphosis, method="class",control=ct,parms = list(prior = c(0.65,0.35), split = "information"))

## 作图有2种方法

## 第一种:

par(mfrow=c(1,3))plot(fit)text(fit,use.n=T,all=T,cex=0.9)

## 第二种,这种会更漂亮一些:

rpart.plot(fit, branch=1, branch.type=2, type=1, extra=102,

shadow.col="gray", box.col="green",

border.col="blue", split.col="red",

split.cex=1.2, main="Kyphosis决策树")

## rpart包提供了复杂度损失修剪的修剪方法,printcp会告诉分裂到每一层,cp是多少,平均相对误差是多少

## 交叉验证的估计误差(“xerror”列),以及标准误差(“xstd”列),平均相对误差=xerror±xstd

printcp(fit)

## 通过上面的分析来确定cp的值

##调用CP(complexity parameter)与xerror的相关图,一种方法是寻找最小xerror点所对应

#的CP值,并由此CP值决定树的大小,另一种方法是利用1SE方法,寻找xerror+SE的最小点对应的CP值。

plotcp(fit)

##利用以下方法进行修剪:

## prune(fit, cp= fit$cptable[which.min(fit$cptable[,"xerror"]),"CP"])

fit2 <- prune(fit, cp=0.01)

#利用模型预测

ndata=data.frame(...)

predict(fit,newdata=ndata)

#案例

str(iris)

set.seed(1234)#设置随机数种子--使每次运行时产生的一组随机数相同,便于结果的重现

#抽样:从iris数据集中随机抽70%定义为训练数据集,30%为测试数据集(常用)

#这里是对行抽样,ind是一个只含1和2的向量

ind <- sample(2, nrow(iris), replace=TRUE, prob=c(0.7, 0.3))

trainData <- iris[ind==1,]

testData <- iris[ind==2,]

f<-Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width

#训练数据

fit<-rpart(f,trainData)

#预测

re<-predict(fit,testData)

#******************或者用其他包********************

library(party)

#建立决策树模型预测花的种类

myFormula <- Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width

iris_ctree <- ctree(myFormula, data=trainData)

# 查看预测的结果

z<-table(predict(iris_ctree), trainData$Species)

#可以根据以上列联表求出预测的正确率---评估模型

#计算准确度

q<-sum(diag(z))/sum(z)

五、机器集成与随机森林法则

前面说过,决策树的一个特点是:弱分类器,分类不完全,需要利用集成投票的方式来增加精确度和稳健性。

机器集成算法:对于数据集训练多个模型,对于分类问题,可以采用投票的方法,选择票数最多的类别作为最终的类别,而对于回归问题,可以采用取均值的方法,取得的均值作为最终的结果。主要的集成算法有bagging和adaboost算法。

随机森林:随机森林就是利用机器集成多个决策树,主要有两个参数,一个是决策树的个数,一个是每棵树的特征变量个数。

随机森林特点:精确度高、稳健性好,但可解释性差。(可以知道各个变量的重要性)

R包实现机器集成算法:

#adabag包均有函数实现bagging和adaboost的分类建模

#利用全部数据建模

library(adabag)

a<-boosting(Species~.,data=iris)

z0<-table(iris[,5],predict(a,iris)$class)

#计算误差率

E0<-(sum(z0)-sum(diag(z0)))/sum(z0)

barplot(a$importance)

b<-errorevol(a,iris)#计算全体的误差演变

plot(b$error,type="l",main="AdaBoost error vs number of trees") #对误差演变进行画图

a<-bagging(Species~.,data=iris)

z0<-table(iris[,5],predict(a,iris)$class)

#计算误差率

E0<-(sum(z0)-sum(diag(z0)))/sum(z0)

barplot(a$importance)

b<-errorevol(a,iris)#计算全体的误差演变

plot(b$error,type="l",main="AdaBoost error vs number of trees") #对误差演变进行画图

#5折交叉验证

set.seed(1044) #设定随机种子

samp=c(sample(1:50,25),sample(51:100,25),sample(101:150,25)) #进行随机抽样

a=boosting(Species~.,data=iris[samp,]) #利用训练集建立adaboost分类模

z0<-table(iris[samp,5],predict(a,iris[samp,])$class)#训练集结果

z1<-table(iris[-samp,5],predict(a,iris[-samp,])$class)#测试集结果

E0<-(sum(z0)-sum(diag(z0)))/sum(z0)

E1<-(sum(z0)-sum(diag(z0)))/sum(z1)

a=bagging(Species~.,data=iris[samp,]) #利用训练集建立adaboost分类模

z0<-table(iris[samp,5],predict(a,iris[samp,])$class)#训练集结果

z1<-table(iris[-samp,5],predict(a,iris[-samp,])$class)#测试集结果

E0<-(sum(z0)-sum(diag(z0)))/sum(z0)

E1<-(sum(z0)-sum(diag(z0)))/sum(z1)

R包实现随机森林:

#随机森林法则

library(randomForest)

library(foreign)

data("iris")

#抽样数据

ind<-sample(2,nrow(iris),replace = TRUE,prob=c(0.7,0.3))

traning<-iris[ind==1,]

testing<-iris[ind==2,]

#训练数据

rf <- randomForest(Species ~ ., data=traning, ntree=100, proximity=TRUE)

#预测

table(predict(rf),traning$Species)

table(predict(rf,testing),testing$Species)

#查看预测的效果

print(rf)

plot(rf)

#查看重要性

importance(rf)

varImpPlot(rf)

# 一、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