使用 R 的非线性回归树

2019-11-18 08:00:00 · 飞浪

介绍回归是一种监督机器学习技术，可以预测连续结果。回归算法有两种类型：线性和非线性。虽然线性模型很有用，但它们依赖于独立变量和因变量之间存在线性关系的假设。在实际业务场景中很难满足这一

介绍

回归是一种监督机器学习技术，可以预测连续结果。回归算法有两种类型：线性和非线性。虽然线性模型很有用，但它们依赖于独立变量和因变量之间存在线性关系的假设。在实际业务场景中很难满足这一假设。这就是非线性回归算法发挥作用的地方，它可以捕捉数据中的非线性。

在本指南中，您将学习如何使用 R 实现非线性回归树。

数据

失业是一个国家的重要社会经济和政治问题，管理失业是任何政府的一项主要任务。在本指南中，我们将构建非线性回归算法来预测经济体中的失业率。

数据来自美国经济时间序列数据（可从https://research.stlouisfed.org/fred2获取）。数据包含 574 行和 5 个变量，如下所述：

psavert：个人储蓄率
pce：个人消费支出，单位：十亿美元
uempmed：失业持续时间中位数，以周为单位
pop：总人口，以百万计
unemploy：失业人数，以百万计。这是因变量。

评估指标

我们将使用两个指标来评估模型的性能：R 平方值和均方根误差 (RMSE)。理想情况下，较低的 RMSE 和较高的 R 平方值表明模型良好。

让我们首先加载所需的库和数据。

      library(plyr)
library(readr)
library(dplyr)
library(caret)
library(ggplot2)
library(rpart)    
library(rpart.plot)
library(randomForest)

dat <- read_csv("reg_data.csv")
glimpse(dat)
    

输出：

      Observations: 574
Variables: 5
$ pce      <dbl> 507.4, 510.5, 516.3, 512.9, 518.1, 525.8, 531.5, 534.2, 544.9...
$ pop      <dbl> 198.7, 198.9, 199.1, 199.3, 199.5, 199.7, 199.8, 199.9, 200.1...
$ psavert  <dbl> 12.5, 12.5, 11.7, 12.5, 12.5, 12.1, 11.7, 12.2, 11.6, 12.2, 1...
$ uempmed  <dbl> 4.5, 4.7, 4.6, 4.9, 4.7, 4.8, 5.1, 4.5, 4.1, 4.6, 4.4, 4.4, 4...
$ unemploy <dbl> 2.9, 2.9, 3.0, 3.1, 3.1, 3.0, 2.9, 3.0, 2.9, 2.7, 2.7, 2.9, 2...
    

输出显示数据集中的所有变量都是数值变量（标记为dbl ）。让我们使用summary()函数查看数据的汇总统计数据。

      summary(dat)

输出：

      pce                 pop              psavert             uempmed         
 Min.   :  507.400   Min.   :198.7000   Min.   : 1.900000   Min.   : 4.000000  
 1st Qu.: 1582.225   1st Qu.:224.8750   1st Qu.: 5.500000   1st Qu.: 6.000000  
 Median : 3953.550   Median :253.0500   Median : 7.700000   Median : 7.500000  
 Mean   : 4843.510   Mean   :257.1913   Mean   : 7.936585   Mean   : 8.610105  
 3rd Qu.: 7667.325   3rd Qu.:290.2500   3rd Qu.:10.500000   3rd Qu.: 9.100000  
 Max.   :12161.500   Max.   :320.9000   Max.   :17.000000   Max.   :25.200000  
    unemploy        
 Min.   : 2.700000  
 1st Qu.: 6.300000  
 Median : 7.500000  
 Mean   : 7.770383  
 3rd Qu.: 8.700000  
 Max.   :15.400000
    

输出显示数据中没有缺失值。它还表明变量的单位和幅度不同，因此需要缩放。

数据分区

我们将在训练集上构建模型，并在测试集上评估其性能。这称为用于评估模型性能的保留验证方法。

下面的第一行代码设置了随机种子，以确保结果的可重复性。第二行创建索引，用于对数据分区进行随机抽样观察。接下来的两行代码创建训练集和测试集，而最后两行打印训练集和测试集的维度。训练集包含 70% 的数据，而测试集包含剩余的 30%。

      set.seed(100) 
index = sample(1:nrow(dat), 0.7*nrow(dat)) 

train = dat[index,] # Create the training data 
test = dat[-index,] # Create the test data

dim(train)
dim(test)
    

输出：

      1] 401   5

[1] 173   5

缩放数值特征

数值特征需要缩放，否则建模过程可能会受到不利影响。

下面的第一行代码创建一个包含独立数字变量名称的列表。第二行使用caret包中的preProcess函数完成缩放任务。预处理对象仅适用于训练数据，而缩放则应用于训练集和测试集。这在下面的第三行和第四行代码中完成。第五行打印缩放后的训练数据集的摘要。输出显示，现在所有数字特征的平均值均为零，除了未缩放的目标变量unemploy 。

      cols = c('pce', 'pop', 'psavert', 'uempmed')

pre_proc_val <- preProcess(train[,cols], method = c("center", "scale"))

train[,cols] = predict(pre_proc_val, train[,cols])
test[,cols] = predict(pre_proc_val, test[,cols])

summary(train)

输出：

      pce               pop             psavert            uempmed           unemploy     
 Min.   :-1.2135   Min.   :-1.6049   Min.   :-1.89481   Min.   :-1.1200   Min.   : 2.700  
 1st Qu.:-0.8856   1st Qu.:-0.8606   1st Qu.:-0.76518   1st Qu.:-0.6198   1st Qu.: 6.500  
 Median :-0.2501   Median :-0.1416   Median :-0.05513   Median :-0.2447   Median : 7.500  
 Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.00000   Mean   : 0.0000   Mean   : 7.791  
 3rd Qu.: 0.7797   3rd Qu.: 0.9060   3rd Qu.: 0.81629   3rd Qu.: 0.1055   3rd Qu.: 8.600  
 Max.   : 2.1244   Max.   : 1.8047   Max.   : 2.91417   Max.   : 4.1571   Max.   :15.300
    

决策树

决策树，也称为分类和回归树 (CART)，适用于分类和连续输入和输出变量。它们的工作原理是根据独立变量中最重要的分割器将数据分成两个或多个同质集。最佳区分是最小化成本指标的区分。分类树的成本指标通常是熵或基尼系数，而对于回归树，默认指标是均方误差。

决策树的基本工作流程如下：

建模过程从根节点开始，根节点代表整个数据。它被分成两个或多个子节点，也称为分裂。分裂过程持续到满足分裂标准为止，发生分裂的子节点称为决策节点。一旦满足分裂标准，节点就不会再分裂；这样的节点称为叶节点或终端节点。我们还可以通过称为修剪的过程删除子节点。

训练模型

现在我们将使用rpart库创建一个回归树模型。第一步是实例化算法，这在下面的第一行代码中完成。参数method = "anova"指定要构建回归模型。对于分类问题，参数将是class。控制参数提供了控制决策树算法训练方式的选项列表。

在下面的代码中，control=rpart.control(minsplit=20, cp=0.001)参数指定在尝试拆分之前，节点中的最小观察数应为 20。它还指定仅当将错误度量降低 0.001 倍（称为成本复杂度因子）时才会发生拆分。第二行打印模型的摘要。请注意，可以使用?rpart()命令了解所有可能的参数。

      tree_model = rpart(unemploy ~ ., data=train, method = "anova", control=rpart.control(minsplit=20, cp=0.001))

summary(tree_model)

输出：

      #Truncated Output 

Call:
rpart(formula = unemploy ~ ., data = train, method = "anova", 
    control = rpart.control(minsplit = 20, cp = 0.001))
  n= 401 

               CP nsplit     rel error        xerror          xstd
0.508609424691      0 1.00000000000 1.00746991170 0.09222517909
0.240199126781      1 0.49139057531 0.53432047449 0.04138079775
0.062388054346      2 0.25119144853 0.29990084346 0.02550320375
0.027180308031      3 0.18880339418 0.21834245878 0.02083102075
0.024078055174      4 0.16162308615 0.20705841419 0.02021514091
0.019456798303      5 0.13754503098 0.17518743853 0.01711287778
0.015928262852      6 0.11808823267 0.16563982394 0.01633076066
0.014943330509      7 0.10215996982 0.15505479787 0.01516085859
0.012183802564      8 0.08721663931 0.13629795408 0.01491501765
0.007010774289      9 0.07503283675 0.12235179176 0.01517832910
0.006864114984     11 0.06101128817 0.11620014096 0.01353111243
0.004679462622     12 0.05414717319 0.11258856185 0.01342153785
0.002606235055     13 0.04946771057 0.09356044395 0.01191539603
0.002392463727     14 0.04686147551 0.08924901899 0.01192405697
0.002373022016     15 0.04446901178 0.08914954140 0.01192548795
0.001075173879     16 0.04209598977 0.08103175712 0.01180204685
0.001054136373     18 0.03994564201 0.07929979931 0.01179186081
0.001039626829     19 0.03889150564 0.07916127996 0.01179072614
0.001024169682     21 0.03681225198 0.07900442265 0.01178836902
0.001000000000     22 0.03578808230 0.07900442265 0.01178836902

Variable importance
    pop uempmed     pce psavert 
    32      28       7
    

根据回归树模型，重要的变量是pop、uempmed和pce。

使用 R 的非线性回归树

评估模型

让我们在训练和测试数据集上评估模型的性能。第一步是创建一个函数来计算评估指标 R 平方和 RMSE。第二步是在训练数据上预测和评估模型，而第三步是在测试数据上预测和评估模型。

      # Step 1 - create the evaluation metrics function

eval_results <- function(true, predicted, df) {
  SSE <- sum((predicted - true)^2)
  SST <- sum((true - mean(true))^2)
  R_square <- 1 - SSE / SST
  RMSE = sqrt(SSE/nrow(df))
  
# Model performance metrics
  data.frame(
    RMSE = RMSE,
    Rsquare = R_square
  )
  
}

# Step 2 - predicting and evaluating the model on train data

predictions_train_cart = predict(tree_model, data = train)
eval_results(train$unemploy, predictions_train_cart, train)

# Step 3 - predicting and evaluating the model on test data

predictions_test_cart = predict(tree_model, newdata = test)
eval_results(test$unemploy, predictions_test_cart, test)
    

输出：

      # Training set results

        RMSE   Rsquare
1 0.4644746 0.9642119"

# Test set results

        RMSE   Rsquare
1 0.6106443 0.9591839
    

上述输出显示，决策树模型在训练数据上的 RMSE 和 R 平方值分别为 0.47 百万和 96.4%。对于测试数据，这些指标的结果分别为 0.61 百万和 96%。

随机森林（或引导聚合）

决策树很有用，但问题是它们往往倾向于过度拟合训练数据，导致测试数据方差较大。随机森林算法通过减少决策树的方差来克服这一缺点。它们之所以被称为森林，是因为它们是几棵决策树的集合或整体。决策树和随机森林模型之间的一个主要区别在于分割方式。在随机森林中，不是尝试对所有特征进行分割，而是为每次分割选择一个特征样本，从而减少模型的方差。

训练模型

在 R 中，randomForest包用于训练随机森林算法。下面的第一行代码实例化了随机森林回归模型，第二行打印了模型的摘要。

      rf_model = randomForest(unemploy ~ ., data=train)
summary(rf_model)
    

输出：

      Length Class  Mode     
call              3    -none- call     
type              1    -none- character
predicted       401    -none- numeric  
mse             500    -none- numeric  
rsq             500    -none- numeric  
oob.times       401    -none- numeric  
importance        4    -none- numeric  
importanceSD      0    -none- NULL     
localImport

_{免责声明：本内容来源于第三方作者授权、网友推荐或互联网整理，旨在为广大用户提供学习与参考之用。所有文本和图片版权归原创网站或作者本人所有，其观点并不代表本站立场。如有任何版权侵犯或转载不当之情况，请与我们取得联系，我们将尽快进行相关处理与修改。感谢您的理解与支持！}

_查看原文

技术指南

R语言

阅读全文

使用 R 的非线性回归树

杭州电子商务研究院

5年前 · 面向社会、服务行业、政产学研结合、整合资源、和谐发展

Observations: 574 Variables: 5 $ pce <dbl> 507.4, 510.5, 516.3, 512.9, 518.1, 525.8, 531.5, 534.2, 544.9... $ pop <dbl> 198.7, 198.9, 199.1, 199.3, 199.5, 199.7, 199.8, 199.9, 200.1... $ psavert <dbl> 12.5, 12.5, 11.7, 12.5, 12.5, 12.1, 11.7, 12.2, 11.6, 12.2, 1... $ uempmed <dbl> 4.5, 4.7, 4.6, 4.9, 4.7, 4.8, 5.1, 4.5, 4.1, 4.6, 4.4, 4.4, 4... $ unemploy <dbl> 2.9, 2.9, 3.0, 3.1, 3.1, 3.0, 2.9, 3.0, 2.9, 2.7, 2.7, 2.9, 2...

pce pop psavert uempmed Min. : 507.400 Min. :198.7000 Min. : 1.900000 Min. : 4.000000 1st Qu.: 1582.225 1st Qu.:224.8750 1st Qu.: 5.500000 1st Qu.: 6.000000 Median : 3953.550 Median :253.0500 Median : 7.700000 Median : 7.500000 Mean : 4843.510 Mean :257.1913 Mean : 7.936585 Mean : 8.610105 3rd Qu.: 7667.325 3rd Qu.:290.2500 3rd Qu.:10.500000 3rd Qu.: 9.100000 Max. :12161.500 Max. :320.9000 Max. :17.000000 Max. :25.200000 unemploy Min. : 2.700000 1st Qu.: 6.300000 Median : 7.500000 Mean : 7.770383 3rd Qu.: 8.700000 Max. :15.400000

cols = c('pce', 'pop', 'psavert', 'uempmed') pre_proc_val <- preProcess(train[,cols], method = c("center", "scale")) train[,cols] = predict(pre_proc_val, train[,cols]) test[,cols] = predict(pre_proc_val, test[,cols]) summary(train)

pce pop psavert uempmed unemploy Min. :-1.2135 Min. :-1.6049 Min. :-1.89481 Min. :-1.1200 Min. : 2.700 1st Qu.:-0.8856 1st Qu.:-0.8606 1st Qu.:-0.76518 1st Qu.:-0.6198 1st Qu.: 6.500 Median :-0.2501 Median :-0.1416 Median :-0.05513 Median :-0.2447 Median : 7.500 Mean : 0.0000 Mean : 0.0000 Mean : 0.00000 Mean : 0.0000 Mean : 7.791 3rd Qu.: 0.7797 3rd Qu.: 0.9060 3rd Qu.: 0.81629 3rd Qu.: 0.1055 3rd Qu.: 8.600 Max. : 2.1244 Max. : 1.8047 Max. : 2.91417 Max. : 4.1571 Max. :15.300

#Truncated Output Call: rpart(formula = unemploy ~ ., data = train, method = "anova", control = rpart.control(minsplit = 20, cp = 0.001)) n= 401 CP nsplit rel error xerror xstd 1 0.508609424691 0 1.00000000000 1.00746991170 0.09222517909 2 0.240199126781 1 0.49139057531 0.53432047449 0.04138079775 3 0.062388054346 2 0.25119144853 0.29990084346 0.02550320375 4 0.027180308031 3 0.18880339418 0.21834245878 0.02083102075 5 0.024078055174 4 0.16162308615 0.20705841419 0.02021514091 6 0.019456798303 5 0.13754503098 0.17518743853 0.01711287778 7 0.015928262852 6 0.11808823267 0.16563982394 0.01633076066 8 0.014943330509 7 0.10215996982 0.15505479787 0.01516085859 9 0.012183802564 8 0.08721663931 0.13629795408 0.01491501765 10 0.007010774289 9 0.07503283675 0.12235179176 0.01517832910 11 0.006864114984 11 0.06101128817 0.11620014096 0.01353111243 12 0.004679462622 12 0.05414717319 0.11258856185 0.01342153785 13 0.002606235055 13 0.04946771057 0.09356044395 0.01191539603 14 0.002392463727 14 0.04686147551 0.08924901899 0.01192405697 15 0.002373022016 15 0.04446901178 0.08914954140 0.01192548795 16 0.001075173879 16 0.04209598977 0.08103175712 0.01180204685 17 0.001054136373 18 0.03994564201 0.07929979931 0.01179186081 18 0.001039626829 19 0.03889150564 0.07916127996 0.01179072614 19 0.001024169682 21 0.03681225198 0.07900442265 0.01178836902 20 0.001000000000 22 0.03578808230 0.07900442265 0.01178836902 Variable importance pop uempmed pce psavert 33 32 28 7

# Step 1 - create the evaluation metrics function eval_results <- function(true, predicted, df) { SSE <- sum((predicted - true)^2) SST <- sum((true - mean(true))^2) R_square <- 1 - SSE / SST RMSE = sqrt(SSE/nrow(df)) # Model performance metrics data.frame( RMSE = RMSE, Rsquare = R_square ) } # Step 2 - predicting and evaluating the model on train data predictions_train_cart = predict(tree_model, data = train) eval_results(train$unemploy, predictions_train_cart, train) # Step 3 - predicting and evaluating the model on test data predictions_test_cart = predict(tree_model, newdata = test) eval_results(test$unemploy, predictions_test_cart, test)

Length Class Mode call 3 -none- call type 1 -none- character predicted 401 -none- numeric mse 500 -none- numeric rsq 500 -none- numeric oob.times 401 -none- numeric importance 4 -none- numeric importanceSD 0 -none- NULL localImport