> library(pacman)
> p_load(dplyr, readxl, caret)

预测任务：利用建筑物的各种特性，例如表面积和屋顶面积，预测建筑物的能源效率，其中效率以供暖负荷和制冷负荷来表示。

1、读取数据

> enb <- read_xlsx("data_set/ENB2012_data.xlsx")
> 
> names(enb) <- c("relcompactness", "surfarea", "wallarea", "roofarea", "height", 
+                 "orientation", "glazarea", "glazareadist", "heatload", 
+                 "coolload")
> 
> DataExplorer::profile_missing(enb)

## # A tibble: 10 x 3
##    feature        num_missing pct_missing
##    <fct>                <int>       <dbl>
##  1 relcompactness           0           0
##  2 surfarea                 0           0
##  3 wallarea                 0           0
##  4 roofarea                 0           0
##  5 height                   0           0
##  6 orientation              0           0
##  7 glazarea                 0           0
##  8 glazareadist             0           0
##  9 heatload                 0           0
## 10 coolload                 0           0

数据集不存在缺失值。

2、转换为虚拟变量

orientation和glazareadist分别表示建筑朝向和玻璃面积分布情况，应该为因子型变量。

> enb <- mutate(enb, across(c(orientation, glazareadist), as.factor))
> 
> str(enb)

## tibble [768 × 10] (S3: tbl_df/tbl/data.frame)
##  $ relcompactness: num [1:768] 0.98 0.98 0.98 0.98 0.9 0.9 0.9 0.9 0.86 0.86 ...
##  $ surfarea      : num [1:768] 514 514 514 514 564 ...
##  $ wallarea      : num [1:768] 294 294 294 294 318 ...
##  $ roofarea      : num [1:768] 110 110 110 110 122 ...
##  $ height        : num [1:768] 7 7 7 7 7 7 7 7 7 7 ...
##  $ orientation   : Factor w/ 4 levels "2","3","4","5": 1 2 3 4 1 2 3 4 1 2 ...
##  $ glazarea      : num [1:768] 0 0 0 0 0 0 0 0 0 0 ...
##  $ glazareadist  : Factor w/ 6 levels "0","1","2","3",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ heatload      : num [1:768] 15.6 15.6 15.6 15.6 20.8 ...
##  $ coolload      : num [1:768] 21.3 21.3 21.3 21.3 28.3 ...

为了让神经网络能够处理这些因子变量，需要先将它们转换为虚拟变量。

> dum <- dummyVars(heatload + coolload ~ ., data = enb)
> enb.new <- predict(dum, newdata = enb) %>% 
+   as_tibble() %>% 
+   bind_cols(enb[, c(9, 10)])
> str(enb.new)

## tibble [768 × 18] (S3: tbl_df/tbl/data.frame)
##  $ relcompactness: num [1:768] 0.98 0.98 0.98 0.98 0.9 0.9 0.9 0.9 0.86 0.86 ...
##  $ surfarea      : num [1:768] 514 514 514 514 564 ...
##  $ wallarea      : num [1:768] 294 294 294 294 318 ...
##  $ roofarea      : num [1:768] 110 110 110 110 122 ...
##  $ height        : num [1:768] 7 7 7 7 7 7 7 7 7 7 ...
##  $ orientation.2 : num [1:768] 1 0 0 0 1 0 0 0 1 0 ...
##  $ orientation.3 : num [1:768] 0 1 0 0 0 1 0 0 0 1 ...
##  $ orientation.4 : num [1:768] 0 0 1 0 0 0 1 0 0 0 ...
##  $ orientation.5 : num [1:768] 0 0 0 1 0 0 0 1 0 0 ...
##  $ glazarea      : num [1:768] 0 0 0 0 0 0 0 0 0 0 ...
##  $ glazareadist.0: num [1:768] 1 1 1 1 1 1 1 1 1 1 ...
##  $ glazareadist.1: num [1:768] 0 0 0 0 0 0 0 0 0 0 ...
##  $ glazareadist.2: num [1:768] 0 0 0 0 0 0 0 0 0 0 ...
##  $ glazareadist.3: num [1:768] 0 0 0 0 0 0 0 0 0 0 ...
##  $ glazareadist.4: num [1:768] 0 0 0 0 0 0 0 0 0 0 ...
##  $ glazareadist.5: num [1:768] 0 0 0 0 0 0 0 0 0 0 ...
##  $ heatload      : num [1:768] 15.6 15.6 15.6 15.6 20.8 ...
##  $ coolload      : num [1:768] 21.3 21.3 21.3 21.3 28.3 ...

3、标准化

在训练神经网络时，为防止饱和现象（因为当优化过程的误差函数的梯度绝对值变得非常小时，非线性神经元激活函数有非常大或非常小的输入，会导致优化过程认为已经达到收敛而终止），需要先对数据进行比例缩放，这样做同时有助于收敛。
将数据维度比例缩放到单位区间[-1, 1]。

> rng <- preProcess(enb.new, method = "range")
> enb.rng <- predict(rng, newdata = enb.new)
> 
> str(enb.rng)

## tibble [768 × 18] (S3: tbl_df/tbl/data.frame)
##  $ relcompactness: num [1:768] 1 1 1 1 0.778 ...
##  $ surfarea      : num [1:768] 0 0 0 0 0.167 ...
##  $ wallarea      : num [1:768] 0.286 0.286 0.286 0.286 0.429 ...
##  $ roofarea      : num [1:768] 0 0 0 0 0.111 ...
##  $ height        : num [1:768] 1 1 1 1 1 1 1 1 1 1 ...
##  $ orientation.2 : num [1:768] 1 0 0 0 1 0 0 0 1 0 ...
##  $ orientation.3 : num [1:768] 0 1 0 0 0 1 0 0 0 1 ...
##  $ orientation.4 : num [1:768] 0 0 1 0 0 0 1 0 0 0 ...
##  $ orientation.5 : num [1:768] 0 0 0 1 0 0 0 1 0 0 ...
##  $ glazarea      : num [1:768] 0 0 0 0 0 0 0 0 0 0 ...
##  $ glazareadist.0: num [1:768] 1 1 1 1 1 1 1 1 1 1 ...
##  $ glazareadist.1: num [1:768] 0 0 0 0 0 0 0 0 0 0 ...
##  $ glazareadist.2: num [1:768] 0 0 0 0 0 0 0 0 0 0 ...
##  $ glazareadist.3: num [1:768] 0 0 0 0 0 0 0 0 0 0 ...
##  $ glazareadist.4: num [1:768] 0 0 0 0 0 0 0 0 0 0 ...
##  $ glazareadist.5: num [1:768] 0 0 0 0 0 0 0 0 0 0 ...
##  $ heatload      : num [1:768] 0.257 0.257 0.257 0.257 0.4 ...
##  $ coolload      : num [1:768] 0.281 0.281 0.281 0.281 0.468 ...

4、拆分训练集和测试集

> set.seed(123)
> ind <- createDataPartition(enb.rng$coolload, p = 0.8, list = F)
> dtrain <- enb.rng[ind, ]
> dtest <- enb.rng[-ind, ]

5、建立神经网络模型

5.1 使用caret包训练

输出变量只选一个，heatload。（暂不清楚caret包能否同时训练两个因变量的模型，以及如何设置参数）

> set.seed(123)
> fit.neur <- train(form = heatload ~ ., method = "neuralnet", data = dtrain[, -18])

5.2 使用neuralnet包训练

输出变量有两个，heatload 和 coolload。10个隐藏层，激活函数选择logistic，误差函数选择sse，它对应的是误差平方和。linear.output = TRUE表示输出层的神经元不应用logistic激活函数，因为这是一个回归任务，需要得到线性的输出，否则输出就会被约束到[0, 1]之间。

> # 使用原始函数可以同时预测两个因变量
> f <- as.formula(paste0("heatload + coolload ~ ", 
+                        paste(names(dtrain)[1:(ncol(dtrain) - 2)],
+                              collapse = " + ")))
> neur <- neuralnet::neuralnet(f, data = dtrain, hidden = 10, 
+                              act.fct = "logistic", linear.output = T,
+                              err.fct = "sse", rep = 1)
> neur$response %>% head()

##    heatload  coolload
## 1 0.2572122 0.2809049
## 2 0.2572122 0.2809049
## 3 0.2572122 0.2809049
## 4 0.4165543 0.3899811
## 5 0.3963332 0.3840560
## 6 0.3685630 0.5036359

6、测试集上的性能

6.1 caret包

> test.pred <- predict(fit.neur, newdata = dtest[, -c(17:18)])
> cor(test.pred, dtest$heatload)

## [1] 0.9966901

预测值与实际值之间的相关性相当高。

6.2 neuralnet包

> test.hat <- predict(neur, newdata = dtest[, -c(17, 18)])
> 
> # 查看预测值与原始值之间的相关度
> cor(test.hat[, 1], dtest$heatload)

## [1] 0.9980793

> cor(test.hat[, 2], dtest$coolload)

## [1] 0.9936744

预测值与实际值之间的相关性非常高，说明模型的性能接近完美。