---
title: "第9章: 决策树和集成学习"
output:
  html_document:
    df_print: default
---

## 统一设置ggplot2的绘图风格
```{r}
library(ggplot2)
theme_set(theme_bw(base_family = "STKaiti"))

```


## 9.1：决策树模型

## 一个决策树的实例

```{r}
library(readr)
library(VIM)
library(caret)
library(rpart)
library(rpart.plot)
library(Metrics)
library(ROCR)
library(readxl)
library(stringr)
library(ggplot2)

## 使用决策树算法，分析泰坦尼克号数据集，
## 读取训练集和测试集

Ttrain <- read_csv("data/chap9/Titanic train.csv")
Ttest <- read_csv("data/chap9/Titanic test.csv")
dim(Ttrain)
colnames(Ttrain)
dim(Ttest)
## 组合数据
Alldata <- rbind.data.frame(Ttrain[,-2],Ttest)
summary(Ttrain)
summary(Ttest)
summary(Alldata)

Survived <- Ttrain$Survived
table(Survived)
## 数据探索与可视化，及特征变换
## 分析数据的缺失值情况，并进行处理
aggr(Alldata)

colnames(Alldata)

# 训练集和测试集要同时作相同的操作
## Cabin缺失值太多，可以直接剔除
Alldata$Cabin <- NULL
## 船票和ID具有识别性所以需要剔除
Alldata$PassengerId <- NULL
Alldata$Ticket <- NULL

## 年龄变量的缺失值可以使用中位数来填补
Alldata$Age[is.na(Alldata$Age)] <- median(Alldata$Age,na.rm = TRUE)
## fare变量的缺失值可以使用均值来填补
Alldata$Fare[is.na(Alldata$Fare)] <- mean(Alldata$Fare,na.rm = TRUE)
## Embarked变量的缺失值，可以使用众数来填补
Embarkedmod <- names(sort(table(Alldata$Embarked),decreasing = T)[1])
Alldata$Embarked[is.na(Alldata$Embarked)] <- Embarkedmod

## 获取新的特征,提取name变量中的特征
newname <- str_split(Alldata$Name," ")
newname <- sapply(newname, function(x) x[2])

sort(table(newname),decreasing = T)
## 名字设置为 Mr.  Miss. Mrs.  Master.，其余的使用other代替
newnamepart <- c("Mr.","Miss.","Mrs.","Master.")
newname[!(newname %in% newnamepart)] <- "other"
Alldata$Name <- as.factor(newname)
Alldata$Sex <- as.factor(Alldata$Sex)
Alldata$Embarked <- as.factor(Alldata$Embarked)
summary(Alldata)
str(Alldata)

```


```{r}

## 有多少兄弟姐妹/配偶同船：'SibSp', 有多少父母/子女同船：'Parch',进行相加生成新的特征
#Alldata$SP <- Alldata$SibSp + Alldata$Parch
#Alldata$SibSp <- NULL
#Alldata$Parch <- NULL
#Alldata$SP <- cut_width(Alldata$SP,3)

#Alldata$Age <- cut_number(Alldata$Age,6)
#Alldata$Fare <- cut_number(Alldata$Fare,5)
#Alldata$Pclass <- as.factor(Alldata$Pclass)

#table(Alldata$Pclass)




```




```{r}

## 与处理好的训练数据和测试数据分开
Ttrainp <- Alldata[1:nrow(Ttrain),]
Ttrainp$Survived <- Survived
Ttestp <- Alldata[(nrow(Ttrain)+1):nrow(Alldata),]
str(Ttrainp)

## 将处理好的训练数据集保存为文件

write.csv(Ttrainp,"data/chap9/Titanic处理后数据.csv",row.names = F)
```



### 数据集切分和模型建立

```{r}

## 切分训练数据集为训练集和测试集，70%训练
set.seed(123)
CDP <- createDataPartition(Ttrainp$Survived,p = 0.8)
train_data <- Ttrainp[CDP$Resample1,]
test_data <- Ttrainp[-CDP$Resample1,]

mod1 <- rpart(Survived~.,data = train_data,method="class",cp = 0.000001)
summary(mod1)
## 看变量重要性
mod1$variable.importance
#cp是每次分割对应的复杂度系数
mod1$cp
## plot cross-validation results
plotcp(mod1)

par(family = "STKaiti")
rpart.plot(mod1, type = 2,extra="auto", under=TRUE, 
           fallen.leaves = FALSE,cex=0.7, main="决策树")

## 查看模型在训练集和测试集上的预测效果
pre_train <- predict(mod1,train_data,type = "prob")
pre_train2<-as.factor(as.vector(ifelse(pre_train[,2]>0.5,1,0)))
pre_test <- predict(mod1,test_data)
pre_test2<-as.factor(as.vector(ifelse(pre_test[,2]>0.5,1,0)))
sprintf("决策树模型在训练集精度为：%f",accuracy(train_data$Survived,pre_train2))
sprintf("决策树模型在测试集精度为：%f",accuracy(test_data$Survived,pre_test2))
## 计算混淆矩阵和模型的精度
cfm <- confusionMatrix(pre_test2,as.factor(test_data$Survived))
cfm$table
cfm



```

    
## 决策树的优化(提高准确率，剪枝等相关操作)

```{r}

bestcp <- mod1$cptable[which.min(mod1$cptable[,"xerror"]),"CP"]
bestcp
# Step3: Prune the tree using the best cp.
mod1.pruned <- prune(mod1, cp = bestcp)

## plot cross-validation results
plotcp(mod1.pruned)
## 可视化剪枝后的决策树
par(family = "STKaiti")
rpart.plot(mod1.pruned, type = 2,extra="auto", under=TRUE, 
           fallen.leaves = FALSE,cex=0.7, main="剪枝后决策树")


## 查看剪枝后模型在训练集和测试集上的预测效果
pre_train_p <- predict(mod1.pruned,train_data)
pre_train_p2<-as.factor(as.vector(ifelse(pre_train_p[,2]>0.5,1,0)))

pre_test_p <- predict(mod1.pruned,test_data)
pre_test_p2<-as.factor(as.vector(ifelse(pre_test_p[,2]>0.5,1,0)))

sprintf("剪枝后决策树模型在训练集精度为：%f",accuracy(train_data$Survived,pre_train_p2))
sprintf("剪枝后决策树模型在测试集精度为：%f",accuracy(test_data$Survived,pre_test_p2))

## 计算混淆矩阵和模型的精度
cfm <- confusionMatrix(pre_test_p2,as.factor(test_data$Survived))
cfm$table
cfm

## 绘制决策树剪枝前后在测试集上的ROC曲线
## 绘制出ROC曲线对比两种模型的效果
## 计算决策树模型的ROC坐标
pr <- prediction(pre_test[,2], test_data$Survived)
prf <- performance(pr, measure = "tpr", x.measure = "fpr")
prfdf <- data.frame(x = prf@x.values[[1]],
                    y = prf@y.values[[1]],
                    model = "rpart")
## 计算剪枝后决策树模型的ROC坐标
pr <- prediction(pre_test_p[,2], test_data$Survived)
prf <- performance(pr, measure = "tpr", x.measure = "fpr")
prfdf2 <- data.frame(x = prf@x.values[[1]],
                    y = prf@y.values[[1]],
                    model = "rpart.prund")

## 合并数据
prfdf <- rbind.data.frame(prfdf,prfdf2)
## plot ROC
ggplot(prfdf,aes(x= x,y = y,colour = model))+
  geom_line(aes(linetype = model),size = 1)+
  theme(aspect.ratio=1)+
  labs(x = "假正例率",y = "真正例率")

## 计算AUC的取值
auc(test_data$Survived,as.vector(pre_test[,2]))
auc(test_data$Survived,as.vector(pre_test_p[,2]))



```


## 9.2：随机森林模型

## 一个实际的数据集作为实例,分类和回归

### 随机森林分类

```{r}
library(tidyr)
library(randomForest)
library(ggRandomForests)
library(caret)

train_data$Survived <- as.factor(train_data$Survived)
rfcla <- randomForest(Survived~.,data = train_data,ntree=200, proximity=TRUE)
summary(rfcla)



## 可视化随机森林的训练过程,随着树的增加训练误差的变化
trainerror <- as.data.frame(plot(rfcla,type = "l"))
colnames(trainerror) <- paste("error",colnames(trainerror),sep = "")
trainerror$ntree <- 1:nrow(trainerror)
trainerror <- gather(trainerror,key = "Type",value = "Error",1:3)
ggplot(trainerror,aes(x = ntree,y = Error))+
  geom_line(aes(linetype = Type,colour = Type))+
  #theme(legend.position = "bottom")+
  ggtitle("随机森林分类模型")+
  theme(plot.title = element_text(hjust = 0.5))
##  模型的精度逐渐趋于平稳

## 或者使用 ggRandomForests包可视化误差
plot(gg_error(rfcla))

## 从randomForest绘制邻近矩阵的标准化坐标。
MDSplot(rfcla,train_data$Survived,k = 2 , palette=c(1, 2),
        pch=20+as.numeric(train_data$Survived))



## 可视化变量的重要性
importance(rfcla)
varImpPlot(rfcla,pch = 20, main = "Importance of Variables")


## 查看模型在测试集上的精度
rfclapre<- predict(rfcla,test_data)
sprintf("随机森林模型测试集精度为：%f",accuracy(test_data$Survived,rfclapre))

```

## 使用所有的训练集训练随机森林模型,预测测试集

```{r}
Ttrainp$Survived <- as.factor(Ttrainp$Survived)
rfclanew <- randomForest(Survived~.,data = Ttrainp,ntree=200, proximity=TRUE)
summary(rfclanew)
rfclanew

## 预测测试集
Ttestpre <- predict(rfclanew,Ttestp)
table(Ttestpre)

```



### 随机森林回归

```{r}

## 使用随机森林的方法，对ENB2012数据进行回归分析
ENB <- read_excel("data/chap9/ENB2012.xlsx")
head(ENB)

summary(ENB)

str(ENB)
## 数据切分为训练集和测试集，训练集70%
set.seed(12)
index <- sample(nrow(ENB),round(nrow(ENB)*0.7))
trainEnb <- ENB[index,]
testENB <- ENB[-index,]
## 建立随机森林回归模型
rfreg <- randomForest(Y1~.,data = trainEnb,ntree=500)
summary(rfreg)

## 可视化模型随着树的增加误差OOB的变化
par(family = "STKaiti")
plot(rfreg,type = "l",col = "red",main = "随机森林回归")

## 使用ggrandomforest包可视化误差
plot(gg_error(rfreg))+labs(title = "随机森林回归")

## 可视化变量的重要性
importance(rfreg)
varImpPlot(rfreg,pch = 20, main = "Importance of Variables")

## 对测试集进行预测，并计算 Mean Squared Error
rfpre <- predict(rfreg,testENB)
sprintf("均方根误差为: %f",mse(testENB$Y1,rfpre))

## 参数搜索,寻找合适的 mtry参数，训练更好的模型
## Tune randomForest for the optimal mtry parameter
set.seed(1234)
rftune <- tuneRF(x = trainEnb[,1:8],y = trainEnb$Y1,
                 stepFactor=1.5,ntreeTry = 500)
print(rftune)
## OOBError误差最小的mtry参数为6

## 建立优化后的随机森林回归模型
rfregbest <- randomForest(Y1~.,data = trainEnb,ntree=500,mtry = 6)

## 可视化两种模型随着树的增加误差OOB的变化
rfregerr <- as.data.frame(plot(rfreg))
colnames(rfregerr) <- "rfregerr"
rfregbesterr <- as.data.frame(plot(rfregbest))
colnames(rfregbesterr) <- "rfregbesterr"
plotrfdata <- cbind.data.frame(rfregerr,rfregbesterr)
plotrfdata$ntree <- 1:nrow(plotrfdata)
plotrfdata <- gather(plotrfdata,key = "Type",value = "Error",1:2)
ggplot(plotrfdata,aes(x = ntree,y = Error))+
  geom_line(aes(linetype = Type,colour = Type),size = 0.9)+
  theme(legend.position = "top")+
  ggtitle("随机森林回归模型")+
  theme(plot.title = element_text(hjust = 0.5))



## 使用优化后的随机森林回归模型，对测试集进行预测，并计算 Mean Squared Error
rfprebest <- predict(rfregbest,testENB)
sprintf("优化后均方根误差为: %f",mse(testENB$Y1,rfprebest))
## 数据准备
index <- order(testENB$Y1)
X <- sort(index)
Y1 <- testENB$Y1[index]
rfpre2 <- rfpre[index]
rfprebest2 <- rfprebest[index]

plotdata <- data.frame(X = X,Y1 = Y1,rfpre =rfpre2,rfprebest = rfprebest2)
plotdata <- gather(plotdata,key="model",value="value",c(-X))

## 可视化模型的预测误差
ggplot(plotdata,aes(x = X,y = value))+
  geom_line(aes(linetype = model,colour = model),size = 0.8)+
  theme(legend.position = c(0.1,0.8),
        plot.title = element_text(hjust = 0.5))+
  ggtitle("随机森林回归模型")

## 随机森林回归的效果非常好

```


## 9.3：梯度提升机

梯度提升模型是回归或分类树模型的集合。 两者都是前向学习集合方法，通过逐步改进的估计获得预测结果。 Boosting是一种灵活的非线性回归程序，有助于提高树模型的准确性。 通过将弱分类算法顺序地应用于递增改变的数据，创建一系列决策树，其产生弱预测模型的集合。 

### 使用GBM算法对泰坦尼克号数据集进行建模

```{r}

library(h2o)
library(Metrics)
library(dplyr)
####################################
##说明：下面启动h2o实例时可能需要更新java版本，在Program\Other文件夹里已提供了jdk-8u221-macosx_64.dmg插件，安装后重新运行即可。
###################################
## 启动初始化一个h2o实例; 定义为2核同时计算；
h2o.init(nthreads=2,max_mem_size='4G')  

## 读取数据
train<-h2o.uploadFile("data/chap9/Titanic处理后数据.csv",
                         destination_frame = "train.hex")

str(train)
train$Survived <- as.factor(train$Survived)
colnames(train)

head(train)



## 数据切分为训练集和测试集,70%训练集，30%测试集
splits <- h2o.splitFrame(data = train, ratios = 0.7,seed = 1234)
train_data <- splits[[1]]
test_data <- splits[[2]]

dim(train_data)

head(train_data)
dim(test_data)

## GBM模型
name1 <- colnames(train)
predictors <- name1[1:8]
target <- "Survived"
gbm <- h2o.gbm(x = predictors, y = target,
               training_frame = train_data,
               distribution="bernoulli", ## 二分类模型
               ntrees = 100,  ## 模型使用数的树量  
               learn_rate=0.01, ## 学习率
               sample_rate = 0.8,## 每棵树使用80%的样本
               col_sample_rate = 0.6,## 每次拆分使用80%的特征
               seed = 1234)
summary(gbm)

## 可视化模型中变量的重要性
h2o.varimp_plot(gbm)
## 计算模型在测试集上的预测值和性能
gbmpre <- as.data.frame(h2o.predict(gbm, newdata = test_data))
head(gbmpre)
acc <- accuracy(as.vector(test_data$Survived),gbmpre$predict)
auc <- h2o.auc(h2o.performance(gbm, newdata = test_data)) 
sprintf("GBM model acc: %f",acc)
sprintf("GBM model AUC: %f",auc)

## 使用参数网格搜索，寻找更合适的模型
ntrees_opt <- c(20,50,100,200,500) ## 树的数量
maxdepth_opt <- c(2,4,6,8,10) ##  树的最大深度
balance_opt <- c(TRUE,FALSE) ## 是否对数据进行类别平衡
hyper_par <- list(ntrees=ntrees_opt, max_depth=maxdepth_opt,
                  balance_classes= balance_opt)
## 使用GBM模型进行超参数搜索
grid <- h2o.grid("gbm", hyper_params = hyper_par,grid_id = "gbm_grid_mol.hex",
                 x = predictors, y = target, distribution="bernoulli",
                 training_frame =train_data,learn_rate=0.01)


## 查看模型的输出
sortedGrid <- h2o.getGrid("gbm_grid_mol.hex", sort_by=c("accuracy"),
                          decreasing = TRUE)    
sortedGrid@summary_table%>%head()


## 将搜索的每个模型均作用于测试集，查看在测试集上的精度
grid_models <- lapply(grid@model_ids, 
                      function(model_id) {model = h2o.getModel(model_id) })

acc <- vector()
modelid <- vector()
for (i in 1:length(grid_models)) {
  gbmpre <- as.data.frame(h2o.predict(grid_models[[i]], newdata = test_data))
  acc[i] <- accuracy(as.vector(test_data$Survived),gbmpre$predict)
  modelid[i] <- grid_models[[i]]@model_id
}
data.frame(modelid = modelid,acc = acc) %>%
  inner_join(sortedGrid@summary_table,by = c("modelid"="model_ids"))%>%
  group_by(modelid)%>%
  arrange(desc(acc))%>%
  head()


```


### GBM模型回归


仍然使用建筑能耗数据集

```{r}
library(h2o)
library(Metrics)
library(dplyr)
library(readxl)
library(tidyr)
library(ggplot2)
## 启动初始化一个h2o实例; 定义为2核同时计算；
h2o.init(nthreads=2,max_mem_size='4G')  
## 使用GBM，对ENB2012数据进行回归分析
ENB <- read_excel("data/chap9/ENB2012.xlsx")
head(ENB)

summary(ENB)

str(ENB)

```

```{r}
## 数据切分为训练集和测试集，训练集70%
set.seed(12)
index <- sample(nrow(ENB),round(nrow(ENB)*0.7))
trainEnb <- as.h2o(ENB[index,])
testENB <- as.h2o(ENB[-index,])

## GBM回归模型
name1 <- colnames(trainEnb)
predictors <- name1[1:8]
target <- "Y1"
## 使用训练集训练一个基础GBM回归模型
gbmreg <- h2o.gbm(x = predictors, y = target,
                  training_frame = trainEnb,
                  distribution="AUTO", ## 回归模型
                  ntrees = 100,seed = 1234)
## 查看在测试集上的模型表达能力
h2o.performance(gbmreg,testENB)
## 使用参数网格搜索，寻找更合适的模型
ntrees_opt <- c(50,100,200,500) ## 树的数量
maxdepth_opt <- c(2,4,6,8,10) ##  树的最大深度
hyper_par <- list(ntrees=ntrees_opt, max_depth=maxdepth_opt)
## 使用GBM模型进行超参数搜索
gbm_grid_reg <- h2o.grid(algorithm="gbm", x = predictors,
                         grid_id ="gbm_grid_reg",
                         y = target,distribution="AUTO",
                         training_frame = trainEnb,hyper_params = hyper_par)



## 查看模型的输出
sortedGrid <- h2o.getGrid("gbm_grid_reg", sort_by="mse", 
                          decreasing = FALSE)
sortedGrid@summary_table%>%head()
## 较好的模型为参数ntree=500,maxdepth=10

## 使用新的参数重新训练模型
gbmreg <- h2o.gbm(x = predictors, y = target,
                  training_frame = trainEnb,
                  distribution="AUTO", ## 回归模型
                  ntrees = 500,  ## 模型使用数的树量  
                  max_depth = 10,seed = 1234)

## 查看在测试集上的预测性能
h2o.performance(gbmreg,newdata = testENB)
## 在测试集上的误差为mse = 0.14
## 可视化模型的预测效果
gbmpre <- as.data.frame(h2o.predict(gbmreg,testENB))
testENBdf <- as.data.frame(testENB)
index <- order(testENBdf$Y1)
X <- sort(index)
Y1 <- testENBdf$Y1[index]
gbmprebest2 <- gbmpre$predict[index]

plotdata <- data.frame(X = X,Y1 = Y1,gbmprebest = gbmprebest2)
plotdata <- gather(plotdata,key="model",value="value",c(-X))

## 可视化模型的预测误差
ggplot(plotdata,aes(x = X,y = value))+
  geom_line(aes(linetype = model,colour = model),size = 0.8)+
  theme(legend.position = c(0.1,0.8),
        plot.title = element_text(hjust = 0.5))+
  ggtitle("GBM回归模型")

```


