R 集成学习实用指南（三）
原文：
annas-archive.org/md5/80cdff39151999bd78d34d1c7cef3132译者：飞龙
第九章. 集成回归模型
第三章, 袋装, 到 第八章, 集成诊断，都致力于学习不同类型的集成方法。讨论主要基于分类问题。如果监督学习问题的回归变量/输出是一个数值变量，那么我们就有了一个回归问题，这将在本章中解决。为了演示目的，本章选择了房价问题，数据集来自 Kaggle 竞赛：www.kaggle.com/c/house-prices-advanced-regression-techniques/。数据包括许多变量，包括多达 79 个独立变量，房价作为输出/依赖变量。数据集需要一些预处理，因为一些变量有缺失日期，一些变量有大量级别，其中一些变量只出现得非常少，还有一些变量在超过 20%的观测值中缺失数据。
预处理技术将被变量减少方法所继替，然后我们将拟合重要的回归模型：线性回归、神经网络和回归树。首先将提供回归树的集成扩展，然后我们将应用袋装和随机森林方法。将使用各种提升方法来提高预测。在结论部分将应用堆叠集成方法。
在本章中，我们将涵盖以下内容：
-
数据预处理和可视化
-
变量减少技术
-
回归模型
-
回归数据的袋装和随机森林
-
提升回归模型
-
回归数据的堆叠集成方法
技术要求
我们将需要以下 R 包来完成本章内容：
-
adabag -
caret -
caretEnsemble -
ClustofVar -
FactoMinR -
gbm -
ipred -
missForest -
nnet -
NeuralNetTools -
plyr -
rpart -
RSADBE
预处理房价数据
数据集是从www.kaggle.com选择的，项目的标题是房价：高级回归技术。我们将使用的主要文件是test.csv和train.csv，这些文件可在配套的捆绑包中找到。变量的描述可以在data_description.txt文件中找到。更详细的信息当然可以在www.kaggle.com/c/house-prices-advanced-regression-techniques/获得。训练数据集包含 1460 个观测值，而测试数据集包含 1459 个观测值。房产价格只在训练数据集中已知，在测试数据集中不可用。我们只将使用训练数据集进行模型开发。首先使用read.csv、dim、names和str函数将数据集加载到 R 会话中，并进行初步检查：
> housing_train <- read.csv("../Data/Housing/train.csv",
+ row.names = 1,na.strings = "NA",
+ stringsAsFactors = TRUE)
> housing_test <- read.csv("../Data/Housing/test.csv",
+ row.names = 1,na.strings = "NA",
+ stringsAsFactors = TRUE)
> dim(housing_train)
[1] 1460 80
> dim(housing_test)
[1] 1459 79
> names(housing_train)
[1] "MSSubClass" "MSZoning" "LotFrontage" "LotArea"
[5] "Street" "Alley" "LotShape" "LandContour"
[9] "Utilities" "LotConfig" "LandSlope" "Neighborhood"
[69] "X3SsnPorch" "ScreenPorch" "PoolArea" "PoolQC"
[73] "Fence" "MiscFeature" "MiscVal" "MoSold"
[77] "YrSold" "SaleType" "SaleCondition" "SalePrice"
> str(housing_train)
'data.frame': 1460 obs. of 80 variables:
$ MSSubClass : int 60 20 60 70 60 50 20 60 50 190 ...
$ MSZoning : Factor w/ 5 levels "C (all)","FV",..: 4 4 4 4 4 4 4 4 5 4 ...
$ LotFrontage : int 65 80 68 60 84 85 75 NA 51 50 ...
$ LotArea : int 8450 9600 11250 9550 14260 14115 10084 10382 6120 7420 ...
$ Street : Factor w/ 2 levels "Grvl","Pave": 2 2 2 2 2 2 2 2 2 2 ...
$ Alley : Factor w/ 2 levels "Grvl","Pave": NA NA NA NA NA NA NA NA NA NA ...
$ MiscFeature : Factor w/ 4 levels "Gar2","Othr",..: NA NA NA NA NA 3 NA 3 NA NA ...
$ MiscVal : int 0 0 0 0 0 700 0 350 0 0 ...
$ MoSold : int 2 5 9 2 12 10 8 11 4 1 ...
$ YrSold : int 2008 2007 2008 2006 2008 2009 2007 2009 2008 2008 ...
$ SaleType : Factor w/ 9 levels "COD","Con","ConLD",..: 9 9 9 9 9 9 9 9 9 9 ...
$ SaleCondition: Factor w/ 6 levels "Abnorml","AdjLand",..: 5 5 5 1 5 5 5 5 1 5 ...
$ SalePrice : int 208500 181500 223500 140000 250000 143000 307000 200000 129900 118000 ...
read.csv function enabled importing the data from the comma-separated values file. The size of the imported data frame is evaluated using the dim function, while names gives us the variable names as stored in the original file. The str function gives a quick preview of the variable types and also gives a few of the observations.
数据框的维度表示变量的数量和观测值的数量。所有变量的详细信息可以在data_description.txt文件中找到。可以看出，我们手头上的是一个综合性的数据集。现在，我们在read.csv导入函数中运行了na.strings = "NA"选项，并且很自然地，这暗示我们存在缺失数据。当训练数据和测试数据分区中都有缺失数据时，作者建议将分区的协变量合并，然后进一步检查。首先将协变量合并，然后我们找出每个变量的缺失观测值数量：
> housing <- rbind(housing_train[,1:79],housing_test)
> dim(housing)
[1] 2919 79
> sort(sapply(housing,function(x) sum(is.na(x))),dec=TRUE)
PoolQC MiscFeature Alley Fence FireplaceQu
2909 2814 2721 2348 1420
LotFrontage GarageYrBlt GarageFinish GarageQual GarageCond
486 159 159 159 159
GarageType BsmtCond BsmtExposure BsmtQual BsmtFinType2
157 82 82 81 80
BsmtFinType1 MasVnrType MasVnrArea MSZoning Utilities
79 24 23 4 2
BsmtFullBath BsmtHalfBath Functional Exterior1st Exterior2nd
2 2 2 1 1
BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF Electrical
1 1 1 1 1
KitchenQual GarageCars GarageArea SaleType MSSubClass
1 1 1 1 0
LotArea Street LotShape LandContour LotConfig
0 0 0 0 0
OpenPorchSF EnclosedPorch X3SsnPorch ScreenPorch PoolArea
0 0 0 0 0
MiscVal MoSold YrSold SaleCondition
0 0 0 0
rbind函数将训练和测试数据集中的数据合并。is.na(x)代码检查x中每个元素的值是否存在，应用sum函数后告诉我们变量的缺失观测值数量。然后使用sapply函数对housing的每个变量应用此函数。使用带有dec=TRUE参数的sort函数按降序对变量的缺失观测值计数进行排序，因此它使我们能够在开始时找到缺失值最多的变量。
读者可能会想知道关于观察结果整理背后的理由。整理观察结果的直观推理是，虽然某些变量可能有缺失数据，在训练数据中比在测试数据中多，或者相反，重要的是整体缺失百分比不要超过观测值的某个特定阈值。尽管我们有缺失数据插补技术，但在缺失数据百分比过高时使用它们可能会使我们错过特征的重要模式。因此，我们任意选择限制变量，如果超过 10%的值缺失。如果任何变量的缺失百分比超过 10%，我们将避免进一步分析该变量。首先，我们识别出超过 10%的变量，然后从主数据框中移除它们。下面的 R 代码块给出了我们想要的结果：
> miss_variables <- names(which(sapply(housing,
+ function(x) sum(is.na(x)))>0.1*nrow(housing_train)))
> miss_variables
[1] "LotFrontage" "Alley" "FireplaceQu" "GarageType"
[5] "GarageYrBlt" "GarageFinish" "GarageQual" "GarageCond"
[9] "PoolQC" "Fence" "MiscFeature"
> length(miss_variables)
[1] 11
> housing[,miss_variables] <- NULL
> dim(housing)
[1] 2919 68
首先识别出超过 10%缺失观测值的变量，然后存储在miss_variables字符向量中，我们有 11 个变量符合这一标准。这些变量通过NULL赋值被消除。
接下来，我们找到因子变量的水平数（不同的）。我们定义了一个函数，find_df，它将找到因子变量的水平数。对于数值和整数变量，它将返回1。这个练习的目的很快就会变得清楚。find_df函数将在下一个代码块中创建：
> find_df <- function(x){
+ if(class(x)=="numeric") mdf <- 1
+ if(class(x)=="integer") mdf <- 1
+ if(class(x) =="factor") mdf <- length(levels(x))
+ if(class(x) =="character") mdf <- length(unique(x))
+ return(mdf)
+ }
> sapply(housing,find_df)
MSSubClass MSZoning LotArea Street LotShape
1 4 1 2 3
LandContour Utilities LotConfig LandSlope Neighborhood
2 3 4 2 25
Condition1 Condition2 BldgType HouseStyle OverallQual
3 2 3 4 1
X3SsnPorch ScreenPorch PoolArea MiscVal MoSold
1 1 1 1 1
YrSold SaleType SaleCondition
1 4 4
> dim(housing)
[1] 2919 68
我们需要检查67个变量，在消除11个具有超过 10%缺失观测值的变量之后。其中一些可能不是因子变量。find_df函数显示，对于因子变量，水平数从 2 到 25 不等。现在对于Condition2和Exterior1st变量出现了一个快速问题：
> round(table(housing$Condition2)/nrow(housing),2)
Artery Feedr Norm PosA PosN RRAe RRAn RRNn
0.00 0.00 0.99 0.00 0.00 0.00 0.00 0.00
> round(table(housing$Exterior1st)/nrow(housing),2)
AsbShng AsphShn BrkComm BrkFace CBlock CemntBd HdBoard ImStucc MetalSd
0.02 0.00 0.00 0.03 0.00 0.04 0.15 0.00 0.15
Plywood Stone Stucco VinylSd Wd Sdng WdShing
0.08 0.00 0.01 0.35 0.14 0.02
在许多实际问题中，似乎存在一些因子变量，它们的一些水平出现频率非常低。现在，如果我们测试/验证分区中有新的水平，我们就无法进行预测。从统计学的角度来看，我们遇到了一个技术问题：失去了太多的自由度。这里采取了一种基本的方法，我们只是简单地将所有观察结果汇总到Others这个总类别中。创建了一个Truncate_Factor函数，它有两个参数：x和alpha。x对象是要传递给函数的变量，而alpha是任何变量频率会被汇总到Others中的指定比例。
注意
如果因子变量在测试数据集中有新的水平，就没有分析方法是能够包含其影响的。因此，在我们有太多不常见水平的情况下，某些水平没有被包括在训练数据集中的可能性很高，预测结果将不会为测试观察结果提供输出。
现在，创建Truncate_Factor函数：
> Truncate_Factor <- function(x,alpha){
+ xc <- as.character(x); n <- length(x)
+ if(length(unique(x))<=20) {
+ critical <- n*alpha
+ xc[xc %in% names(which((prop.table(table(xc)))<alpha))] <- "Others"
+ }
+ xc <- as.factor(xc)
+ return(xc)
+ }
> for(i in 1:ncol(housing)){
+ if(any(class(housing[,i]) == c('character','factor')))
+ housing[,i] = Truncate_Factor(housing[,i],0.05)
+ }
> table(housing$Condition2)/nrow(housing)
Norm Others
0.99 0.01
> table(housing$Exterior1st)/nrow(housing)
HdBoard MetalSd Others Plywood VinylSd Wd Sdng
0.151 0.154 0.126 0.076 0.351 0.141
我们现在可以看到，“其他”级别出现的频率更高，如果我们随机创建分区，未知级别的出现问题很可能不会发生。
你可能还记得，我们之前已经消除了具有过多缺失观测值的变量。这并不意味着我们没有缺失数据，这一点可以很快地发现：
> sum(is.na(housing))
[1] 474
> prod(dim(housing))
[1] 198492
474个值不能被忽视。缺失数据插补是填充缺失值的重要方法。尽管 EM 算法是实现这一目标的一种流行方法，但我们将应用随机森林技术来模拟缺失观测值。missForest包在第四章中介绍，即随机森林，并使用一个示例来模拟缺失值。我们将应用此函数到房价数据框上。由于此函数中默认选择的变量数是mtry=5，而我们房价中有 68 个变量，因此用于分割节点的变量数将更改为大约 p/3，因此在下一个 R 块中可以看到mtry=20的选项。在 8GB RAM 的机器上，下一行代码需要运行几个小时。接下来，我们将应用missForest函数，保存插补对象以供将来参考，并使用插补值创建测试和训练数据集：
> housing_impute <- missForest(housing,maxiter = 10,ntree=500,mtry=20)
missForest iteration 1 in progress...done!
missForest iteration 2 in progress...done!
missForest iteration 3 in progress...done!
missForest iteration 4 in progress...done!
missForest iteration 5 in progress...done!
missForest iteration 6 in progress...done!
missForest iteration 7 in progress...done!
There were 14 warnings (use warnings() to see them)
> save(housing_impute,file=
+ '../Data/Housing/housing_covariates_impute.Rdata')
> ht_imp <- cbind(housing_impute$ximp[1:nrow(housing_train),],
+ housing_train$SalePrice)
> save(ht_imp,file='../Data/Housing/ht_imp.Rdata')
> htest_imp <- housing_impute$ximp[(nrow(housing_train)+1):nrow(
+ housing),]
> save(htest_imp,file='../Data/Housing/htest_imp.Rdata')
读者当然应该在他们的本地机器上运行missForest代码行。然而，为了节省时间，读者也可以跳过这一行，然后从代码包中加载ht_imp和htest_imp对象。下一节将展示一种可视化大数据集和两种数据降维方法的方式。
可视化和变量降维
在上一节中，房价数据经历了大量的分析预处理，我们现在准备进一步分析它。首先，我们从可视化开始。由于我们有大量的变量，R 可视化设备上的可视化稍微有些困难。正如前几章所见，为了可视化随机森林和其他大型、复杂结构，我们将初始化一个 PDF 设备并将图表存储在其中。在房价数据集中，主要变量是房价，因此我们将首先将输出变量命名为SalePrice。我们需要以方便展示众多变量与SalePrice之间的关系的方式来可视化数据。自变量可以是数值型或分类型。如果变量是数值型，散点图将表明变量与SalePrice回归量之间的关系。如果自变量是分类型/因子型，我们将对因子的每个级别可视化箱线图。pdf、plot和boxplot函数将有助于生成所需的图表：
> load("../Data/Housing/ht_imp_author.Rdata")
> names(ht_imp)[69] <- "SalePrice"
> SP <- ht_imp$SalePrice
> pdf("../Output/Visualizing_Housing_Data.pdf")
> for(i in 1:68){
+ if(class(ht_imp[,i])=="numeric") {
+ plot(ht_imp[,i],SP,xlab=names(ht_imp)[i],ylab="Sales Price")
+ title(paste("Scatter plot of Sales Price against ",
+ names(ht_imp)[i]))
+ }
+ if(class(ht_imp[,i])=="factor") {
+ boxplot(SP~ht_imp[,i],xlab=names(ht_imp)[i],ylab=
+ "Sales Price",notch=TRUE)
+ title(paste("Boxplot of Salesprice by ",names(ht_imp)[i]))
+ }
+ }
> dev.off()
null device
1
ht_imp对象是从ht_imp_author.Rdata文件中加载的。注意，如果你在自己的机器上运行missForest函数并在此文件上工作，那么结果将与ht_imp_author.Rdata不同。pdf函数已知可以启动一个同名文件，如之前多次所见。对于数值变量，检查if条件，并显示散点图，其中xlab使用变量的实际名称作为x轴上的标签名称。title函数将paste函数的输出贴上，而paste函数确保我们有一个适合生成的图的标题。对于因子变量，也测试了类似的条件。我们现在将查看一些有趣的图表。SalePrice与MSSubClass的第一个图表（见Visualizing_Housing_Data.pdf文件）如下：
图 1：销售价格与 MSSubClass 的散点图
注意
注意，尽管我们指定了MSSubClass变量为数值变量，但散点图并没有给出同样的印象。在这里，MSSubClass变量的值围绕一个特定的点杂乱无章，然后尺度跳到下一个值。
简而言之，它似乎不是一个连续变量，这可以通过以下方式轻松验证：
> table(ht_imp$MSSubClass)
20 30 40 45 50 60 70 75 80 85 90 120 160 180 190
536 69 4 12 144 299 60 16 58 20 52 87 63 10 30
练习：读者应将MSSubClass变量转换为因子，然后应用Truncate_Factor以减少噪声。在Visualizing_Housing_Data.pdf文件中识别其他表现出此特性的数值变量。
现在让我们看看MSZoning因子变量的箱线图：
图 2：MSZoning 三个级别的销售价格箱线图
超出胡须的点表示存在异常值。然而，对于复杂问题，解释也可能非常错误。缺口是箱线图显示中的有用技巧。如果两个变量级别的缺口不重叠，则表示级别是显著的，因此信息是有用的，如SalePrice与MSZoning级别箱线图的显示所示。
接下来展示的是SalePrice与LotArea的散点图：
图 3：销售价格与 Lot Area 的散点图
明显地，散点图显示这两个变量SalePrice和LotArea之间没有有意义的关系。在下面的图中，SalePrice与TotalBsmtSF之间可以看到不同类型的显示：
图 4：销售价格与 TotalBsmtSF 的散点图
我们可以清楚地看到图的最右侧TotalBsmtSF值存在一个异常值。同时，TotalBsmtSF在 0 值处也存在值的杂乱，这可能是其他某个变量的控制结果。或者，可能会发现该变量存在零膨胀，因此它可能是一个混合变量。同样，所有其他图表也可以进行解释。接下来获得SalePrice与其他数值变量之间的相关性：
> cor(ht_imp[sapply(ht_imp,is.numeric)])[,1]
MSSubClass LotArea OverallQual OverallCond YearBuilt
1.0000 -0.1398 0.0326 -0.0593 0.0279
YearRemodAdd MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF
0.0406 0.0206 -0.0698 -0.0656 -0.1408
TotalBsmtSF X1stFlrSF X2ndFlrSF LowQualFinSF GrLivArea
-0.2385 -0.2518 0.3079 0.0465 0.0749
BsmtFullBath BsmtHalfBath FullBath HalfBath BedroomAbvGr
0.0035 -0.0023 0.1316 0.1774 -0.0234
KitchenAbvGr TotRmsAbvGrd Fireplaces GarageCars GarageArea
0.2817 0.0404 -0.0456 -0.0401 -0.0987
WoodDeckSF OpenPorchSF EnclosedPorch X3SsnPorch ScreenPorch
-0.0126 -0.0061 -0.0120 -0.0438 -0.0260
PoolArea MiscVal MoSold YrSold SalePrice
0.0083 -0.0077 -0.0136 -0.0214 -0.0843
练习：解释Visualizing_Housing_Data.pdf文件中的所有关系，并按照其绝对值在前面 R 代码中对相关性进行排序。
我们使用了感兴趣的变量进行可视化，这进而导致了有价值的见解。正如之前所述，p = 68意味着有很多协变量/独立变量。在大数据中，复杂性将在北方方向增加，而且众所周知，对于许多实际应用，我们拥有数千个独立变量。虽然大多数可视化技术都有洞察力，但一个缺点是我们很少能深入了解高阶关系。例如，当涉及到三个或更多变量时，在图形显示中很少能充分展现它们之间的关系。因此，部署那些在不牺牲信息的情况下减少变量数量的方法是重要的。这里将要讨论的两种数据降维方法是主成分分析和变量聚类。
主成分分析（PCA）是从多元统计的更广泛领域抽取的一种方法。在数据降维方面，这种方法很有用，因为它在给定原始变量数量的情况下，试图用尽可能少的新变量来覆盖原始数据的大部分方差。这里简要介绍了 PCA。
假设我们有一个观测值的随机向量 。给定随机向量
，主成分分析（PCA）找到一个新的主成分向量
，使得每个 Yi 都是
的线性组合。此外，主成分满足 Y1 的方差高于 Y2 的方差，且两者不相关；Y2 的方差高于 Y3 和 Y1 的方差，Y2 和 Y3 不相关，以此类推。这与
有关，它们之间都不相关。主成分被设置成使得
的大部分方差累积在前几个主成分中（关于此信息的更多信息，请参阅 Tattar 等人（2016）的第十五章）。因此，我们可以实现大量的数据降维。然而，PCA 的基本前提是
是一个连续随机变量的向量。在我们的数据集中，我们也有因子变量。因此，我们不能为了我们的目的使用 PCA。一种粗略的方法是忽略因子变量，并在连续变量上简单地执行数据降维。相反，我们会使用
混合数据的因子分析，并且 FactoMineR 软件包中提供了执行此操作的软件功能。
由于数据降维只需要在协变量上执行，而我们没有纵向数据，因此数据降维应用于所有可用的观测值集，而不仅仅是训练数据集。在整份数据集上执行数据降维的理由与截断因子变量级别的数量相同。housing_impute 数据框可在 housing_covariates_impute.Rdata 中找到。我们首先将其加载，然后应用 FAMD 函数来执行混合数据的因子分析：
> load("../Data/Housing/housing_covariates_impute.Rdata")
> housing_covariates <- housing_impute$ximp
> housing_cov_famd <- FAMD(housing_covariates,ncp=68,graph=FALSE)
> colnames(housing_cov_famd$eig) <- c("Component","Variance",
+ "Cumulative")
> housing_cov_famd$eig
Component Variance Cumulative
comp 1 12.2267562274 9.3334017003 9.33340170
comp 2 5.4502085801 4.1604645650 13.49386627
comp 3 4.5547218487 3.4768869074 16.97075317
comp 4 4.0710151565 3.1076451576 20.07839833
comp 5 3.1669428163 2.4175136002 22.49591193
comp 6 2.8331129142 2.1626816139 24.65859354
comp 7 2.6471571767 2.0207306692 26.67932421
comp 8 2.1871762983 1.6696002277 28.34892444
comp 9 2.1563067109 1.6460356572 29.99496010
comp 10 2.0083000432 1.5330534681 31.52801357
comp 66 0.7691341212 0.5871252834 80.58667899
comp 67 0.7648033308 0.5838193365 81.17049833
comp 68 0.7559712365 0.5770772798 81.74757561
> windows(height=100,width=200)
> pareto.chart(housing_cov_famd$eig[,2])
在FAMD函数中，ncp选项被设置为 68，因为这是我们拥有的变量数量。我们还想看看主成分如何响应数据集。如果选择graph=TRUE选项，函数将显示相关图表。housing_cov_famd$eig的colnames被更改为默认名称，因为这些名称并不能公正地反映它生成的输出。我们可以从特征值分析中看到，总共 68 个成分并没有完全覆盖数据中可用的全部变异。此外，即使是解释了 50%变异的成分，我们也需要从中选择 26 个。因此，这里的变量减少似乎并不富有成效。然而，这并不意味着在下一组分析中性能会差。当使用质量控制包qcc中的pareto.chart函数对频率数据进行处理时，会得到一个帕累托图。正如百分比所展示的，很明显，如果我们需要原始变量中 90%的变异由主成分解释，那么我们需要近 60 个主成分。因此，减少的变量数量仅为 8，解释也增加了一个复杂性。这并不是好消息。然而，我们仍然会保存主成分的数据：
> save(housing_cov_famd,file='../Data/Housing/Housing_FAMD.Rdata')
> Housing_FAMD_Data <- housing_cov_famd$ind$coord
> save(Housing_FAMD_Data,file='../Data/Housing/
+ Housing_FAMD_Data.Rdata')
图 5：主成分贡献的帕累托图
练习：探索使用来自 R 包PCAmix的PCAmix函数，通过主成分分析来减少变量的数量。
变量聚类
变量可以像我们对观测值那样进行分组。为了实现这一点，我们将使用来自ClustOfVar包的kmeansvar函数。变量聚类包需要分别指定定量（数值）变量，定性（因子）变量也需要分别指定。此外，我们还需要指定需要多少个变量簇。init选项有助于在这里进行指定。is.numeric和is.factor函数用于识别数值和因子变量，并设置变量簇：
> Housing_VarClust <- kmeansvar(
+ X.quanti = housing_covariates[sapply(housing_covariates,
+ is.numeric)],
+ X.quali = housing_covariates[sapply(housing_covariates,
+ is.factor)],init=4)
Error: Some categorical variables have same names of categories,
rename categories or use the option rename.level=TRUE to rename it automatically
哎呀！这是一个错误。重要的是要记住，所有不常见的因子变量级别都被标记为“其他”。可能存在其他级别在多个变量中具有相同的名称，这在调查数据中是一个非常常见的标签选择，包括如“非常不满意 < 不满意 < 一般 < 好 < 优秀”这样的选项。这种变量级别的选择可以在多个问题中相同。然而，我们需要所有变量的级别名称在所有变量中都是独特的。手动重命名标签将是徒劳的，并且是时间上的极大浪费。因此，我们将使用一组将在变量中唯一的名称来解决这个问题，即变量名称本身。变量名称将与变量级别连接，因此我们将在整个变量中拥有独特的因子级别。使用paste0函数和plyr包中的mapvalues，我们首先执行级别重命名操作，然后再次应用kmeansvar：
> hc2 <- housing_covariates
> for(i in 1:ncol(hc2)){
+ if(class(hc2[,i])=="factor") {
+ hc2[,i] <- mapvalues(hc2[,i],from=levels(hc2[,i]),
+ to=paste0(names(hc2)[i],"_",levels(hc2[,i])))
+ }
+ }
> Housing_VarClust <- kmeansvar(
+ X.quanti = hc2[sapply(hc2,is.numeric)],
+ X.quali = hc2[sapply(hc2,is.factor)], init=4)
> Housing_VarClust$cluster
MSSubClass LotArea OverallQual OverallCond YearBuilt
2 1 1 4 4
YearRemodAdd MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF
4 3 1 2 4
BsmtCond BsmtExposure BsmtFinType1 BsmtFinType2 Heating
3 1 4 2 3
HeatingQC CentralAir Electrical KitchenQual Functional
4 1 4 4 4
PavedDrive SaleType SaleCondition
1 4 4
> summary(Housing_VarClust)
Call:
kmeansvar(X.quanti = hc2[sapply(hc2, is.numeric)], X.quali = hc2[sapply(hc2, is.factor)], init = 4)
number of iterations: 2
Data:
number of observations: 2919
number of variables: 68
number of numerical variables: 34
number of categorical variables: 34
number of clusters: 4
Cluster 1 :
squared loading correlation
X1stFlrSF 0.6059 0.778
TotalBsmtSF 0.5913 0.769
OverallQual 0.5676 0.753
PoolArea 0.0166 0.129
MiscVal 0.0059 0.077
MoSold 0.0024 0.049
Cluster 2 :
squared loading correlation
X2ndFlrSF 0.8584 -0.927
HouseStyle 0.7734 NA
TotRmsAbvGrd 0.5185 -0.720
BsmtFinType2 0.0490 NA
BsmtFinSF2 0.0408 0.202
X3SsnPorch 0.0039 0.063
Cluster 3 :
squared loading correlation
MasVnrType 0.83189 NA
MasVnrArea 0.82585 -0.909
Heating 0.03532 NA
BsmtCond 0.02681 NA
Utilities 0.00763 NA
YrSold 0.00084 0.029
Cluster 4 :
squared loading correlation
Neighborhood 0.7955 NA
YearBuilt 0.7314 -0.855
BsmtQual 0.6792 NA
BsmtHalfBath 0.0087 0.093
Street 0.0041 NA
Condition2 0.0015 NA
Gain in cohesion (in %): 11.56
重要的问题是，尽管我们已经将变量分组标记，但我们如何使用它们？答案是每个组内变量的系数。要显示系数，请在clustvar对象Housing_VarClust旁边运行$coef：
> Housing_VarClust$coef
$cluster1
[,1]
const -7.1e+00
LotArea 2.1e-05
OverallQual 2.2e-01
CentralAir_N -5.3e-01
CentralAir_Y 3.8e-02
PavedDrive_N -5.2e-01
PavedDrive_Others -2.8e-01
PavedDrive_Y 5.0e-02
$cluster2
[,1]
const 3.79789
MSSubClass -0.00472
BsmtFinSF2 0.00066
HouseStyle_1.5Fin -0.11967
HouseStyle_1Story 0.41892
HouseStyle_2Story -0.69610
HouseStyle_Others 0.10816
BsmtFinType2_Others 0.33286
BsmtFinType2_Unf -0.04491
$cluster3
[,1]
const -33.1748
MasVnrArea -0.0039
YrSold 0.0167
BsmtCond_TA -0.0365
Heating_GasA -0.0179
Heating_Others 1.1425
$cluster4
[,1]
const 45.30644
OverallCond 0.09221
YearBuilt -0.01009
SaleCondition_Normal 0.03647
SaleCondition_Others 0.20598
SaleCondition_Partial -0.58877
现在，对于数据中的观测值，相应的变量将与变量聚类的系数相乘，以获得该变量聚类的单个向量。因此，我们将 68 个变量减少到 4 个变量。
练习：使用之前显示的系数，获取housing_covariates数据框的聚类变量。
房地产问题的数据预处理现在已完成。在下一节中，我们将为回归数据构建基础学习器。
回归模型
弗朗西斯·高尔顿爵士在十九世纪末发明了简单的线性回归模型。所用的例子是研究父母的身高如何影响孩子的身高。这项研究使用了数据，并奠定了回归分析的基础。父母和孩子的身高之间的相关性是众所周知的，高尔顿利用 928 对身高测量数据，开发了线性回归模型。然而，在伽尔顿正式发明之前，这种方法可能已经在非正式地使用了。简单的线性回归模型由一个单一输入（自变量）和一个单一输出组成。
在这种监督学习方法中，目标变量/输出/因变量是一个连续变量，它也可以取区间值，包括非负数和实数。输入/自变量没有限制，因此它可以取数值、分类或其他我们之前用于分类问题的任何形式。有趣的是，线性回归模型的出现时间比分类回归模型（如逻辑回归模型）要早得多。机器学习问题更常基于分类问题进行概念化，而集成方法，特别是提升方法，是通过将分类作为动机来开发的。主要原因在于误差改进提供了良好的直觉，次要原因可能是因为著名的机器学习示例，如数字识别、垃圾邮件分类等。
简单线性回归的扩展是多元线性回归，其中我们允许有多个独立变量。我们将完全放弃简单和多元回归的惯例，并简单地遵循回归。作为基学习器，线性回归模型首先被介绍。有趣的数据库将被用来启动线性回归模型。
线性回归模型
更正式地说，设为包含p个独立变量的集合，而Y是感兴趣的变量。我们需要从回归变量
的角度来理解回归变量Y。线性回归模型由以下公式给出：
Y 与回归变量之间的关系是线性的；是截距项；
是回归系数；而
是误差项。需要指出的是，这里的线性是指回归系数的线性。同时，需要注意的是，回归变量可以采取任何形式，有时也可以被看作是其他形式，包括对数、指数和二次形式。误差项
通常假设服从具有未知方差和零均值的正态分布。关于线性回归模型的更多细节，可以参考 Draper 和 Smith (1999)、Chatterjee 和 Hadi (2012)以及 Montgomery 等人(2005)的著作。关于使用 R 软件实现此技术的信息，请参阅 Tattar 等人(2016)的第十二章或 Tattar(2017)的第六章。
首先，我们将使用 Galton 数据集来解释线性回归模型的核心概念。数据是从RSADBE包中加载的，使用lm函数，我们可以构建模型：
> data(galton)
> cor(galton)
child parent
child 1.00 0.46
parent 0.46 1.00
> plot(galton)
> head(galton)
child parent
1 62 70
2 62 68
3 62 66
4 62 64
5 62 64
6 62 68
> cp_lm <- lm(child~parent,data=galton)
> summary(cp_lm)
Call:
lm(formula = child ~ parent, data = galton)
Residuals:
Min 1Q Median 3Q Max
-7.805 -1.366 0.049 1.634 5.926
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 23.9415 2.8109 8.52 <2e-16 ***
parent 0.6463 0.0411 15.71 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.2 on 926 degrees of freedom
Multiple R-squared: 0.21, Adjusted R-squared: 0.21
F-statistic: 247 on 1 and 926 DF, p-value: <2e-16
图 6：儿童身高与父母身高对比 - 散点图
这段代码块告诉我们什么？首先，我们将从 RSADBE 包中加载 galton 数据，然后查看父母和孩子的cor相关系数。相关系数是0.46，这似乎是一个强烈的正相关。散点图也表明了这种正相关，因此我们继续构建以父母高度为函数的孩子的线性回归模型。建议首先查看与模型相关的 p 值，在这种情况下，它在summary(cp_lm)的最后行给出，为<2e-16。较小的 p 值意味着我们拒绝模型无意义的零假设，因此当前的拟合模型是有用的。与截距项和变量项相关的 p 值都是<2e-16，这也意味着这些项是显著的。回归系数0.6463意味着如果父母身高增加一英寸，孩子的身高将增加回归系数的量级。
Multiple R-squared（技术上简单的 R 平方）和Adjusted R-squared的值都是0.21，这是预期的，因为我们模型中只有一个变量。R 平方的解释是，如果我们将其乘以 100（在这种情况下是 21%），那么得到的数字就是数据（孩子的高度）变化中由拟合值解释的百分比。这个指标的值越高，模型就越好。在这个例子中，这意味着父母的高度只能解释孩子身高变化的约 21%。这意味着我们需要考虑其他变量。在这种情况下，一个起点可能是考虑父母双方的高度。如果你添加更多的变量，多重 R 平方值将会持续增加，因此更倾向于使用更稳健的调整 R 平方值。是否有可能获得完美的 R 平方，例如 1 或 100%？
一个名为Mantel的数据集在捆绑包中在线可用，我们将构建一个线性回归模型来检查其 R 平方。为此，我们导入数据集并对其运行lm函数：
> Mantel <- read.csv("../Data/Mantel.csv")
> Mantel
Y X1 X2 X3
1 5 1 1004 6.0
2 6 200 806 7.3
3 8 -50 1058 11.0
4 9 909 100 13.0
5 11 506 505 13.1
> Mantel_lm <- lm(Y~.,data=Mantel)
> summary(Mantel_lm)
Call:
lm(formula = Y ~ ., data = Mantel)
Residuals:
1 2 3 4 5
-2.49e-13 2.92e-13 3.73e-14 -3.89e-14 -4.14e-14
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.00e+03 2.73e-10 -3.67e+12 1.7e-13 ***
X1 1.00e+00 2.73e-13 3.67e+12 1.7e-13 ***
X2 1.00e+00 2.73e-13 3.67e+12 1.7e-13 ***
X3 1.33e-14 2.16e-13 6.00e-02 0.96
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.9e-13 on 1 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 4.99e+25 on 3 and 1 DF, p-value: 1.04e-13
在这里，我们可以看到 R 平方是完美的。在我们开始分析房价数据的严肃任务之前，让我们先玩一会儿吧。
对于galton数据集，我们将添加一个新的变量，称为frankenstein，这个变量将是拟合模型cp_lm的残差。将创建一个新的数据集，它将使用残差来增强galton数据集；然后，将使用lm函数拟合线性模型，并检查其 R 平方：
> d2 <- cbind(galton,residuals(cp_lm))
> names(d2)
[1] "child" "parent" "residuals(cp_lm)"
> names(d2) <- c("child","parent","frankenstein")
> cpf_lm <- lm(child~.,d2)
> summary(cpf_lm)
Call:
lm(formula = child ~ ., data = d2)
Residuals:
Min 1Q Median 3Q Max
-2.60e-15 -7.40e-16 -3.00e-16 2.10e-16 1.02e-13
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.39e+01 5.74e-15 4.17e+15 <2e-16 ***
parent 6.46e-01 8.40e-17 7.69e+15 <2e-16 ***
frankenstein 1.00e+00 6.71e-17 1.49e+16 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.6e-15 on 925 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 1.41e+32 on 2 and 925 DF, p-value: <2e-16
Warning message:
In summary.lm(cpf_lm) : essentially perfect fit: summary may be unreliable
不要忽略警告函数。你可能还记得，对于Mantel数据集没有显示这样的警告函数。这是因为通过向frankenstein变量添加一点噪声，可以消除这个警告，从而使它更加可怕：
> d2$frankenstein <- jitter(d2$frankenstein)
> summary(lm(child~.,d2))
Call:
lm(formula = child ~ ., data = d2)
Residuals:
Min 1Q Median 3Q Max
-0.004072 -0.002052 0.000009 0.001962 0.004121
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.39e+01 2.92e-03 8210 <2e-16 ***
parent 6.46e-01 4.27e-05 15143 <2e-16 ***
frankenstein 1.00e+00 3.41e-05 29331 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.0023 on 925 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 5.45e+08 on 2 and 925 DF, p-value: <2e-16
因此，我们已经掌握了获得完美 R 平方的艺术。玩耍的时间已经结束，让我们继续到住房数据集。我们之前已经将住房数据集保存为训练和测试块，分别命名为ht_imp.Rdata和htest_imp.Rdata文件。作者将文件名版本修改为_author，以便更清晰地表示。然后，我们将训练块分为训练和测试两部分。然后，我们使用load函数导入数据，使用sample函数对其进行分区，然后使用lm函数构建回归模型：
> load("../Data/Housing/ht_imp_author.Rdata")
> load("../Data/Housing/htest_imp_author.Rdata")
> ls()
[1] "ht_imp" "htest_imp"
> Y <- "SalePrice"
> X <- names(ht_imp)[-69]
> set.seed(12345)
> BV <- sample(c("Build","Validate"),nrow(ht_imp),replace = TRUE,
+ prob=c(0.7,0.3))
> HT_Build <- ht_imp[BV=="Build",]
> HT_Validate <- ht_imp[BV=="Validate",]
> HT_Formula <- as.formula("SalePrice~.")
> HT_LM_01 <- lm(HT_Formula,data=HT_Build)
> summary(HT_LM_01)
Call:
lm(formula = HT_Formula, data = HT_Build)
Residuals:
Min 1Q Median 3Q Max
-268498 -12222 -409 11351 240990
Coefficients: (2 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.87e+03 1.53e+06 0.00 0.99850
MSSubClass -1.52e+02 7.95e+01 -1.91 0.05583 .
MSZoningRL 8.55e+03 6.27e+03 1.36 0.17317
MSZoningRM 1.20e+04 7.50e+03 1.60 0.11011
LotArea 4.90e-01 1.21e-01 4.04 5.8e-05 ***
StreetPave 2.81e+04 1.70e+04 1.65 0.09979 .
LotShapeOthers -3.59e+03 6.12e+03 -0.59 0.55733
LotShapeReg 1.25e+03 2.40e+03 0.52 0.60111
LandContourOthers -1.22e+04 3.99e+03 -3.05 0.00236 **
UtilitiesOthers -5.76e+04 3.25e+04 -1.77 0.07637 .
LotConfigCulDSac 1.21e+04 4.96e+03 2.44 0.01477 *
LotConfigInside -1.62e+03 2.58e+03 -0.63 0.52972
LotConfigOthers -1.28e+04 5.57e+03 -2.30 0.02144 *
EnclosedPorch 6.95e+00 1.91e+01 0.36 0.71628
X3SsnPorch 3.81e+01 3.87e+01 0.98 0.32497
ScreenPorch 3.78e+01 2.01e+01 1.88 0.05988 .
PoolArea 5.13e+01 2.60e+01 1.98 0.04842 *
MiscVal 5.13e-02 6.57e+00 0.01 0.99377
MoSold -4.38e+02 3.67e+02 -1.19 0.23313
YrSold -1.01e+02 7.53e+02 -0.13 0.89376
SaleTypeOthers -4.88e+04 2.19e+04 -2.23 0.02598 *
SaleTypeWD -5.10e+04 2.20e+04 -2.32 0.02061 *
SaleConditionNormal 1.93e+03 4.31e+03 0.45 0.65421
SaleConditionOthers 1.87e+03 7.42e+03 0.25 0.80168
SaleConditionPartial -3.21e+04 2.21e+04 -1.45 0.14641
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 28400 on 861 degrees of freedom
Multiple R-squared: 0.884, Adjusted R-squared: 0.867
F-statistic: 51.1 on 129 and 861 DF, p-value: <2e-16
在拟合三个更多基础学习器之后，将进行拟合线性模型的准确性评估。调整后的 R 平方值约为 87%。然而，我们有 68 个变量，我们可以从之前的总结中的 p 值看到，许多变量没有 p 值小于 0.05 或 0.1。因此，我们需要消除不显著的变量。step函数可以应用于许多拟合的回归模型，以消除不显著的变量，同时保留模型的大部分特征。
在 R 会话中运行step函数会导致控制台中出现大量的输出。初始输出因空间限制而丢失。因此，作者使用 RStudio 中的从 R 脚本编译报告选项运行脚本，选择 MS Word 作为报告输出格式，并保存该文件。以下是该文件结果的简略版：
## Start: AIC=20446.87
## SalePrice ~ MSSubClass + MSZoning + LotArea + Street + LotShape +
## LandContour + Utilities + LotConfig + LandSlope + Neighborhood +
## Condition1 + Condition2 + BldgType + HouseStyle + OverallQual +
## OverallCond + YearBuilt + YearRemodAdd + RoofStyle + RoofMatl +
## Exterior1st + Exterior2nd + MasVnrType + MasVnrArea + ExterQual +
## ExterCond + Foundation + BsmtQual + BsmtCond + BsmtExposure +
## BsmtFinType1 + BsmtFinSF1 + BsmtFinType2 + BsmtFinSF2 + BsmtUnfSF +
## TotalBsmtSF + Heating + HeatingQC + CentralAir + Electrical +
## X1stFlrSF + X2ndFlrSF + LowQualFinSF + GrLivArea + BsmtFullBath +
## BsmtHalfBath + FullBath + HalfBath + BedroomAbvGr + KitchenAbvGr +
## KitchenQual + TotRmsAbvGrd + Functional + Fireplaces + GarageCars +
## GarageArea + PavedDrive + WoodDeckSF + OpenPorchSF + EnclosedPorch +
## X3SsnPorch + ScreenPorch + PoolArea + MiscVal + MoSold +
## YrSold + SaleType + SaleCondition
##
##
## Step: AIC=20446.87
## SalePrice ~ MSSubClass + MSZoning + LotArea + Street + LotShape +
## LandContour + Utilities + LotConfig + LandSlope + Neighborhood +
## Condition1 + Condition2 + BldgType + HouseStyle + OverallQual +
## OverallCond + YearBuilt + YearRemodAdd + RoofStyle + RoofMatl +
## Exterior1st + Exterior2nd + MasVnrType + MasVnrArea + ExterQual +
## ExterCond + Foundation + BsmtQual + BsmtCond + BsmtExposure +
## BsmtFinType1 + BsmtFinSF1 + BsmtFinType2 + BsmtFinSF2 + BsmtUnfSF +
## TotalBsmtSF + Heating + HeatingQC + CentralAir + Electrical +
## X1stFlrSF + X2ndFlrSF + LowQualFinSF + BsmtFullBath + BsmtHalfBath +
## FullBath + HalfBath + BedroomAbvGr + KitchenAbvGr + KitchenQual +
## TotRmsAbvGrd + Functional + Fireplaces + GarageCars + GarageArea +
## PavedDrive + WoodDeckSF + OpenPorchSF + EnclosedPorch + X3SsnPorch +
## ScreenPorch + PoolArea + MiscVal + MoSold + YrSold + SaleType +
## SaleCondition
##
##
## Step: AIC=20446.87
## SalePrice ~ MSSubClass + MSZoning + LotArea + Street + LotShape +
## LandContour + Utilities + LotConfig + LandSlope + Neighborhood +
## Condition1 + Condition2 + BldgType + HouseStyle + OverallQual +
## OverallCond + YearBuilt + YearRemodAdd + RoofStyle + RoofMatl +
## Exterior1st + Exterior2nd + MasVnrType + MasVnrArea + ExterQual +
## ExterCond + Foundation + BsmtQual + BsmtCond + BsmtExposure +
## BsmtFinType1 + BsmtFinSF1 + BsmtFinType2 + BsmtFinSF2 + BsmtUnfSF +
## Heating + HeatingQC + CentralAir + Electrical + X1stFlrSF +
## X2ndFlrSF + LowQualFinSF + BsmtFullBath + BsmtHalfBath +
## FullBath + HalfBath + BedroomAbvGr + KitchenAbvGr + KitchenQual +
## TotRmsAbvGrd + Functional + Fireplaces + GarageCars + GarageArea +
## PavedDrive + WoodDeckSF + OpenPorchSF + EnclosedPorch + X3SsnPorch +
## ScreenPorch + PoolArea + MiscVal + MoSold + YrSold + SaleType +
## SaleCondition
##
## Df Sum of Sq RSS AIC
## - Exterior2nd 5 2.6926e+09 6.9890e+11 20441
## - HeatingQC 3 8.4960e+08 6.9706e+11 20442
## - MasVnrType 3 9.3578e+08 6.9714e+11 20442
## - OverallQual 1 3.2987e+10 7.2919e+11 20491
## - X2ndFlrSF 1 3.9790e+10 7.3600e+11 20500
## - Neighborhood 24 1.6770e+11 8.6391e+11 20613
##
## Step: AIC=20440.69
## SalePrice ~ MSSubClass + MSZoning + LotArea + Street + LotShape +
## LandContour + Utilities + LotConfig + LandSlope + Neighborhood +
## Condition1 + Condition2 + BldgType + HouseStyle + OverallQual +
## OverallCond + YearBuilt + YearRemodAdd + RoofStyle + RoofMatl +
## Exterior1st + MasVnrType + MasVnrArea + ExterQual + ExterCond +
## Foundation + BsmtQual + BsmtCond + BsmtExposure + BsmtFinType1 +
## BsmtFinSF1 + BsmtFinType2 + BsmtFinSF2 + BsmtUnfSF + Heating +
## HeatingQC + CentralAir + Electrical + X1stFlrSF + X2ndFlrSF +
## LowQualFinSF + BsmtFullBath + BsmtHalfBath + FullBath + HalfBath +
## BedroomAbvGr + KitchenAbvGr + KitchenQual + TotRmsAbvGrd +
## Functional + Fireplaces + GarageCars + GarageArea + PavedDrive +
## WoodDeckSF + OpenPorchSF + EnclosedPorch + X3SsnPorch + ScreenPorch +
## PoolArea + MiscVal + MoSold + YrSold + SaleType + SaleCondition
## Step: AIC=20386.81
## SalePrice ~ MSSubClass + LotArea + Street + LandContour + Utilities +
## LotConfig + LandSlope + Neighborhood + Condition1 + Condition2 +
## BldgType + HouseStyle + OverallQual + OverallCond + YearBuilt +
## RoofStyle + RoofMatl + Exterior1st + BsmtQual + BsmtCond +
## BsmtExposure + BsmtFinType1 + BsmtFinSF1 + BsmtFinType2 +
## X1stFlrSF + X2ndFlrSF + LowQualFinSF + BsmtFullBath + FullBath +
## HalfBath + KitchenAbvGr + KitchenQual + TotRmsAbvGrd + Functional +
## Fireplaces + GarageCars + WoodDeckSF + ScreenPorch + PoolArea +
## MoSold + SaleType
##
## Df Sum of Sq RSS AIC
## <none> 7.1467e+11 20387
## - KitchenAbvGr 1 1.4477e+09 7.1612e+11 20387
## - MoSold 1 1.6301e+09 7.1630e+11 20387
## - BldgType 2 3.1228e+09 7.1779e+11 20387
## - Utilities 1 1.7130e+09 7.1639e+11 20387
## - BsmtCond 1 1.7554e+09 7.1643e+11 20387
## - BsmtFinType2 1 1.8708e+09 7.1654e+11 20387
## - YearBuilt 1 2.0543e+09 7.1673e+11 20388
## - Street 1 2.1163e+09 7.1679e+11 20388
## - LowQualFinSF 1 2.1785e+09 7.1685e+11 20388
## - ScreenPorch 1 2.2387e+09 7.1691e+11 20388
## - MSSubClass 1 2.2823e+09 7.1695e+11 20388
## - LandSlope 1 2.5566e+09 7.1723e+11 20388
## - PoolArea 1 2.6036e+09 7.1728e+11 20388
## - Exterior1st 5 9.1221e+09 7.2379e+11 20389
## - Functional 1 3.4117e+09 7.1808e+11 20390
## - Condition1 2 4.9604e+09 7.1963e+11 20390
## - BsmtFinSF1 1 3.9442e+09 7.1862e+11 20390
## - Condition2 1 4.0659e+09 7.1874e+11 20390
## - RoofStyle 2 6.1817e+09 7.2085e+11 20391
## - HalfBath 1 5.3010e+09 7.1997e+11 20392
## - FullBath 1 5.4987e+09 7.2017e+11 20392
## - Fireplaces 1 6.0438e+09 7.2072e+11 20393
## - TotRmsAbvGrd 1 7.0166e+09 7.2169e+11 20395
## - LandContour 1 7.7036e+09 7.2238e+11 20395
## - WoodDeckSF 1 8.8947e+09 7.2357e+11 20397
## - LotConfig 3 1.2015e+10 7.2669e+11 20397
## - RoofMatl 1 9.0967e+09 7.2377e+11 20397
## - BsmtFullBath 1 9.4178e+09 7.2409e+11 20398
## - HouseStyle 3 1.2940e+10 7.2761e+11 20399
## - BsmtFinType1 5 1.7704e+10 7.3238e+11 20401
## - SaleType 2 1.5305e+10 7.2998e+11 20404
## - LotArea 1 1.4293e+10 7.2897e+11 20404
## - OverallCond 1 1.8131e+10 7.3280e+11 20410
## - BsmtQual 3 2.3916e+10 7.3859e+11 20413
## - X1stFlrSF 1 2.1106e+10 7.3578e+11 20414
## - BsmtExposure 3 2.8182e+10 7.4285e+11 20419
## - GarageCars 1 2.6886e+10 7.4156e+11 20421
## - KitchenQual 3 3.1267e+10 7.4594e+11 20423
## - OverallQual 1 3.7361e+10 7.5203e+11 20435
## - X2ndFlrSF 1 4.3546e+10 7.5822e+11 20443
## - Neighborhood 24 1.8921e+11 9.0389e+11 20572
model的总结如下：
> summary(HT_LM_Final)
Call:
lm(formula = SalePrice ~ MSSubClass + LotArea + Street + LandContour +
Utilities + LotConfig + LandSlope + Neighborhood + Condition1 +
Condition2 + BldgType + HouseStyle + OverallQual + OverallCond +
YearBuilt + RoofStyle + RoofMatl + Exterior1st + BsmtQual +
BsmtCond + BsmtExposure + BsmtFinType1 + BsmtFinSF1 + BsmtFinType2 +
X1stFlrSF + X2ndFlrSF + LowQualFinSF + BsmtFullBath + FullBath +
HalfBath + KitchenAbvGr + KitchenQual + TotRmsAbvGrd + Functional +
Fireplaces + GarageCars + WoodDeckSF + ScreenPorch + PoolArea +
MoSold + SaleType, data = HT_Build)
Residuals:
Min 1Q Median 3Q Max
-272899 -11717 -42 11228 235349
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.64e+05 1.78e+05 -1.48 0.13894
MSSubClass -1.27e+02 7.46e+01 -1.70 0.08965 .
LotArea 4.75e-01 1.12e-01 4.25 2.3e-05 ***
MoSold -4.99e+02 3.48e+02 -1.44 0.15136
SaleTypeOthers -1.69e+04 5.85e+03 -2.89 0.00396 **
SaleTypeWD -1.76e+04 4.00e+03 -4.40 1.2e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 28100 on 904 degrees of freedom
Multiple R-squared: 0.881, Adjusted R-squared: 0.87
F-statistic: 78.1 on 86 and 904 DF, p-value: <2e-16
覆盖step函数的小模块可以在Housing_Step_LM.R文件中找到，使用 R Markdown 生成的输出保存在名为Housing_Step_LM.docx的文件中。step函数的输出超过四十三页，但我们不必检查每一步中省略的变量。只需说，已经消除了许多不显著的变量，而没有丢失模型的特征。验证分区的准确性评估将在稍后看到。接下来，我们将扩展线性回归模型到非线性模型，并解决神经网络问题。
练习：使用主成分和变量聚类变量构建线性回归模型。使用相关变量集的准确性（R 平方）是否提高了线性回归模型？
神经网络
神经网络架构在第一章的统计/机器学习模型部分中介绍，集成技术介绍。神经网络能够处理非线性关系，隐藏神经元数量的选择、传递函数的选择以及学习率（或衰减率）的选择在构建有用的回归模型时提供了很大的灵活性。Haykin (2009)和 Ripley (1996)提供了神经网络理论的两个详细解释。
我们已经探讨了神经网络在分类问题中的应用，并看到了堆叠集成模型的实际操作。对于回归模型，我们需要通过linout=TRUE选项告诉nnet函数输出/因变量是连续变量。在这里，我们将构建一个具有五个隐藏神经元的神经网络，size=5，并运行函数最多 100 次迭代，maxit=100：
> HT_NN <- nnet(HT_Formula,data=HT_Build,linout=TRUE,maxit=100,size=5)
# weights: 666
initial value 38535430702344.617187
final value 5951814083616.587891
converged
> summary(HT_NN)
a 131-5-1 network with 666 weights
options were - linear output units
b->h1 i1->h1 i2->h1 i3->h1 i4->h1 i5->h1 i6->h1 i7->h1
-1.0e-02 6.5e-01 -8.0e-02 4.6e-01 5.0e-02 -4.0e-02 3.9e-01 1.3e-01
i8->h1 i9->h1 i10->h1 i11->h1 i12->h1 i13->h1 i14->h1 i15->h1
2.1e-01 4.6e-01 1.9e-01 5.2e-01 -6.6e-01 3.2e-01 -3.0e-02 2.2e-01
i16->h1 i17->h1 i18->h1 i19->h1 i20->h1 i21->h1 i22->h1 i23->h1
-2.5e-01 -1.2e-01 3.3e-01 -2.8e-01 -4.6e-01 -3.8e-01 -4.1e-01 -3.2e-01
-4.0e-01 -2.9e-01 -5.1e-01 -2.6e-01 2.5e-01 -6.0e-01 1.0e-02 1.5e-01
i120->h5 i121->h5 i122->h5 i123->h5 i124->h5 i125->h5 i126->h5 i127->h5
3.7e-01 -2.0e-01 2.0e-01 1.0e-02 -3.3e-01 -2.4e-01 -1.9e-01 7.0e-01
i128->h5 i129->h5 i130->h5 i131->h5
-1.3e-01 -3.4e-01 -6.9e-01 -6.6e-01
b->o h1->o h2->o h3->o h4->o h5->o
6.3e+04 6.3e+04 6.3e+04 -9.1e+04 4.7e-01 -8.4e+03
注意，神经网络架构并不非常实用。然而，有时我们被要求展示我们所构建的内容。因此，我们将使用来自NeuralNetTools包的plotnet函数来生成网络。由于变量太多（本例中为 68 个），我们将绘图保存到Housing_NN.pdf PDF 文件中，读者可以打开并放大查看：
> pdf("../Output/Housing_NN.pdf",height = 25, width=60)
> plotnet(HT_NN) # very chaotic network
> dev.off()
RStudioGD
2
神经网络的预测将很快进行。
练习 1：使用不同的衰减选项构建神经网络；默认值为 0。在 0-0.2 的范围内变化衰减值，增量分别为 0.01、0.05 等。
练习 2：使用reltol值、衰减值以及这些变量的组合来改进神经网络拟合。
回归树
回归树构成了住房数据集的第三个基学习器，并为回归问题提供了决策树结构。决策树的自然优势也自然地转移到回归树上。如第三章中所述，Bagging，许多决策树的选项也适用于回归树。
我们将使用rpart库中的rpart函数和默认设置来构建回归树。使用绘图和文本函数，我们设置了回归树：
> HT_rtree <- rpart(HT_Formula,data=HT_Build)
> windows(height=100,width=100)
> plot(HT_rtree,uniform = TRUE); text(HT_rtree)
> HT_rtree$variable.importance
OverallQual Neighborhood YearBuilt ExterQual KitchenQual
3.2e+12 2.0e+12 1.7e+12 1.7e+12 1.4e+12
Foundation GarageCars GrLivArea GarageArea X1stFlrSF
1.3e+12 8.0e+11 6.9e+11 6.1e+11 3.8e+11
X2ndFlrSF TotalBsmtSF TotRmsAbvGrd BsmtQual MasVnrArea
3.8e+11 3.2e+11 2.7e+11 2.7e+11 1.8e+11
FullBath HalfBath HouseStyle BsmtFinSF1 YearRemodAdd
1.7e+11 1.3e+11 1.2e+11 1.1e+11 5.3e+10
MSZoning BsmtFinType1 BedroomAbvGr Exterior1st BsmtFullBath
4.6e+10 4.4e+10 4.0e+10 2.4e+10 1.1e+10
LotArea
5.7e+09
图 7：房屋销售价格的回归树
哪些变量在这里很重要？这个答案由变量重要性指标提供。我们从HT_rtree中提取变量重要性，条形长度最高的变量是所有变量中最重要的。现在我们将使用barplot函数对HT_rtree进行绘图：
> barplot(HT_rtree$variable.importance,las=2,yaxt="n")
图 8：住房模型回归树的变量重要性
练习：探索回归树的剪枝选项。
接下来，我们将查看三个基础学习器在验证数据集上的性能。
回归模型的预测
我们将房价训练数据集分为两个部分：训练和验证。现在我们将使用构建的模型并检查它们的性能如何。我们将通过查看 MAPE 指标（|实际-预测|/实际）来完成这项工作。使用带有newdata选项的predict函数，首先获得预测值，然后计算数据验证部分观测值的 MAPE：
> HT_LM_01_val_hat <- predict(HT_LM_01,newdata = HT_Validate[,-69])
Warning message:
In predict.lm(HT_LM_01, newdata = HT_Validate[, -69]) :
prediction from a rank-deficient fit may be misleading
> mean(abs(HT_LM_01_val_hat - HT_Validate$SalePrice)/HT_Validate$SalePrice)
[1] 0.11
> HT_LM_Final_val_hat <- predict(HT_LM_Final,newdata = HT_Validate[,-69])
> mean(abs(HT_LM_Final_val_hat - HT_Validate$SalePrice)/HT_Validate$SalePrice)
[1] 0.11
> HT_NN_val_hat <- predict(HT_NN,newdata = HT_Validate[,-69])
> mean(abs(HT_NN_val_hat - HT_Validate$SalePrice)/HT_Validate$SalePrice)
[1] 0.37
> HT_rtree_val_hat <- predict(HT_rtree,newdata = HT_Validate[,-69])
> mean(abs(HT_rtree_val_hat - HT_Validate$SalePrice)/HT_Validate$SalePrice)
[1] 0.17
线性回归模型HT_LM_01和最有效的线性模型（根据 AIC）HT_LM_Final都给出了相同的准确度（精确到两位数），这两个模型的 MAPE 为0.11。具有五个隐藏神经元的神经网络模型HT_NN的 MAPE 为0.37，这是一个不好的结果。这再次证实了众所周知的事实：复杂性并不一定意味着准确性。回归树HT_rtree的准确度为0.17。
预测的价格在以下程序中进行了可视化：
> windows(height = 100,width = 100)
> plot(HT_Validate$SalePrice,HT_LM_01_val_hat,col="blue",
+ xlab="Sales Price",ylab="Predicted Value")
> points(HT_Validate$SalePrice,HT_LM_Final_val_hat,col="green")
> points(HT_Validate$SalePrice,HT_NN_val_hat,col="red")
> points(HT_Validate$SalePrice,HT_rtree_val_hat,col="yellow")
> legend(x=6e+05,y=4e+05,lty=3,
+ legend=c("Linear","Best Linear","Neural Network","Regression Tree"),
+ col=c("blue","green","red","yellow"))
图 9：预测房价
现在我们已经设置了基础学习器，是时候构建它们的集成模型了。我们将基于决策树的同质基础学习器构建集成模型。
Bagging 和随机森林
第三章，Bagging，和第四章，随机森林，展示了如何提高基本决策树的稳定性和准确性。在本节中，我们将主要使用决策树作为基础学习器，并以与第三章和第四章中相同的方式创建树集成。
split函数是 Bagging 和随机森林算法在分类和回归树中的主要区别。因此，不出所料，我们可以继续使用与分类问题中使用的相同函数和包来处理回归问题。我们将首先使用ipred包中的bagging函数来为房价数据设置 Bagging 算法：
> housing_bagging <- bagging(formula = HT_Formula,data=ht_imp,nbagg=500,
+ coob=TRUE,keepX=TRUE)
> housing_bagging$err
[1] 35820
与第三章中的相同方式一样，Bagging 对象中的树可以被保存到 PDF 文件中：Bagging：
> pdf("../Output/Housing_Bagging.pdf")
> for(i in 1:500){
+ temp <- housing_bagging$mtrees[[i]]
+ plot(temp$btree)
+ text(temp$btree,use.n=TRUE)
+ }
> dev.off()
RStudioGD
2
由于 ipred 包没有直接给出变量重要性，而了解哪些变量重要始终是一个重要的衡量标准，因此我们运行了一个与 第三章 中使用的类似循环和程序，即 Bagging，以获取变量重要性图：
> VI <- data.frame(matrix(0,nrow=500,ncol=ncol(ht_imp)-1))
> vnames <- names(ht_imp)[-69]
> names(VI) <- vnames
> for(i in 1:500){
+ VI[i,] <- as.numeric(housing_bagging$mtrees[[i]]$btree$variable.importance[vnames])
+ }
> Bagging_VI <- colMeans(VI,na.rm = TRUE)
> Bagging_VI <- sort(Bagging_VI,dec=TRUE)
> barplot(Bagging_VI,las=2,yaxt="n")
> title("Variable Importance of Bagging")
图 10：住房数据 bagging 算法的变量重要性图
练习： 比较 图 10 和 图 8，以决定回归树中是否存在过拟合问题。
bagging 是否提高了预测性能？这是我们需要评估的重要标准。使用带有 newdata 选项的 predict 函数，我们再次计算 MAPE，如下所示：
> HT_bagging_val_hat <- predict(housing_bagging,newdata = HT_Validate[,-69])
> mean(abs(HT_bagging_val_hat - HT_Validate$SalePrice)/HT_Validate$SalePrice)
[1] 0.13
简单回归树的 MAPE 为 17%，现在降至 13%。这引出了下一个练习。
练习： 使用 rpart.control 中的某些剪枝选项来提高 bagging 的性能。
在 bagging 之后，下一步是随机森林。我们将使用同名的 randomForest 函数。在这里，我们探索了 500 棵树来构建这个森林。对于回归数据，节点分裂时随机采样的协变量数量的默认设置是 mtry = p/3，其中 p 是协变量的数量。我们将使用默认选择。randomForest 函数用于设置树集成，然后使用在 第四章 中定义的 plot_rf 函数，即 随机森林，将森林的树保存到 PDF 文件中：
> housing_RF <- randomForest(formula=HT_Formula,data=ht_imp,ntree=500,
+ replace=TRUE,importance=TRUE)
> pdf("../Output/Housing_RF.pdf",height=100,width=500)
Error in pdf("../Output/Housing_RF.pdf", height = 100, width = 500) :
cannot open file '../Output/Housing_RF.pdf'
> plot_RF(housing_RF)
[1] 1
[1] 2
[1] 3
[1] 498
[1] 499
[1] 500
> dev.off()
null device
1
> windows(height=100,width=200)
> varImpPlot(housing_RF2)
下图给出了随机森林的变量重要性图：
图 11：住房数据的随机森林变量重要性
练习： 找出两个变量重要性图 %lncMSE 和 IncNodePurity 之间的差异。同时，比较随机森林的变量重要性图和 bagging 图，并对此进行评论。
我们的森林有多准确？使用 predict 函数，我们将得到答案：
> HT_RF_val_hat <- predict(housing_RF,newdata = HT_Validate[,-69])
> mean(abs(HT_RF_val_hat - HT_Validate$SalePrice)/HT_Validate$SalePrice)
[1] 0.038
这真是太棒了，随机森林通过大幅降低 MAPE 从 0.17 到 0.038 显著提高了精度。这是迄今为止构建的所有模型中的绝对赢家。
练习： 尽管精度有所提高，但尝试基于剪枝树构建森林，并计算精度。
让我们看看提升如何改变树的性能。
提升回归模型
第五章《提升》介绍了当有感兴趣的分类变量时，提升法对树的提升方法。将提升法应用于回归问题需要大量的计算改变。更多信息，请参阅 Zemel 和 Pitassi（2001）的论文，papers.nips.cc/paper/1797-a-gradient-based-boosting-algorithm-for-regression-problems.pdf，或 Ridgeway 等人（1999）的论文，dimacs.rutgers.edu/Research/MMS/PAPERS/BNBR.pdf。
将使用gbm库中的gbm函数来提升由随机森林生成的弱学习器。我们生成了一千棵树，n.trees=1e3，并使用shrinkage因子为0.05，然后使用梯度提升算法对回归数据进行提升：
> housing_gbm <- gbm(formula=HT_Formula,data=HT_Build,distribution = "gaussian",
+ n.trees=1e3,shrinkage = 0.05,keep.data=TRUE,
+ interaction.depth=1,
+ cv.folds=3,n.cores = 1)
> summary(housing_gbm)
var rel.inf
OverallQual OverallQual 29.22608012
GrLivArea GrLivArea 18.85043432
Neighborhood Neighborhood 13.79949556
PoolArea PoolArea 0.00000000
MiscVal MiscVal 0.00000000
YrSold YrSold 0.00000000
本总结按降序给出了变量的重要性。可以通过gbm.perf函数查看提升法的性能，由于我们的目标始终是生成对新数据表现良好的技术，因此还叠加了以下出袋曲线：
> windows(height=100,width=200)
> par(mfrow=c(1,2))
> gbm.perf(housing_gbm,method="OOB",plot.it=TRUE,
+ oobag.curve = TRUE,overlay=TRUE)
[1] 135
图 12：住房数据的提升收敛
提升法在第137次迭代时收敛。接下来，我们查看提升法在验证数据上的性能：
> HT_gbm_val_hat <- predict(housing_gbm,newdata = HT_Validate[,-69])
Using 475 trees...
> mean(abs(HT_gbm_val_hat - HT_Validate$SalePrice)/HT_Validate$SalePrice)
[1] 0.11
MAPE 已从 17%降至 11%。然而，随机森林仍然是迄今为止最准确的模型。
回归模型的堆叠方法
线性回归模型、神经网络和回归树是这里将要堆叠的三个方法。我们将需要caret和caretEnsemble包来完成这项任务。堆叠集成方法已在第七章《一般集成技术》中详细介绍。首先，我们指定训练任务的控制参数，指定算法列表，并创建堆叠集成：
> control <- trainControl(method="repeatedcv", number=10, repeats=3,
+ savePredictions=TRUE, classProbs=TRUE)
> algorithmList <- c('lm', 'rpart')
> set.seed(12345)
> Emodels <- caretList(HT_Formula, data=HT_Build, trControl=control,
+ methodList=algorithmList,
+ tuneList=list(
+ nnet=caretModelSpec(method='nnet', trace=FALSE,
+ linout=TRUE)
+
+ )
+ )
There were 37 warnings (use warnings() to see them)
神经网络通过caretModelSpec指定。Emodels需要重新采样以进行进一步分析：
> Enresults <- resamples(Emodels)
> summary(Enresults)
Call:
summary.resamples(object = Enresults)
Models: nnet, lm, rpart
Number of resamples: 30
MAE
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
nnet 30462 43466 47098 47879 53335 58286 0
lm 16153 18878 20348 20138 21337 23865 0
rpart 30369 33946 35688 35921 37354 42437 0
RMSE
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
nnet 42598 66632 70197 69272 73089 85971 0
lm 22508 26137 29192 34347 39803 66875 0
rpart 38721 46508 50528 50980 55705 65337 0
Rsquared
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
nnet 0.0064 0.16 0.32 0.31 0.44 0.74 4
lm 0.4628 0.77 0.85 0.81 0.88 0.92 0
rpart 0.4805 0.55 0.57 0.58 0.61 0.69 0
> dotplot(Enresults)
接下来显示dotplot：
图 13：住房数据的 R-square、MAE 和 RMSE
我们可以从图 13中看到，三个模型的 R-square 相似，尽管三个模型的 MAE 和 RMSE 有显著差异。可以使用modelCor函数找到模型的相关性：
> modelCor(Enresults)
nnet lm rpart
nnet 1.000 0.033 -0.44
lm 0.033 1.000 0.29
rpart -0.441 0.288 1.00
我们现在将集成方法应用于验证数据：
> HT_Validate_Predictions <- rowMeans(predict(Emodels,newdata = HT_Validate))
Warning message:
In predict.lm(modelFit, newdata) :
prediction from a rank-deficient fit may be misleading
> mean(abs(HT_Validate_Predictions - HT_Validate$SalePrice)/HT_Validate$SalePrice)
[1] 0.16
注意
注意，神经网络的默认结果，我们没有指定隐藏层的大小。16%的 MAPE 并不理想，我们最好使用随机森林集成。
练习：对主成分和变量聚类数据进行堆叠集成方法的操作。
摘要
在本章中，我们扩展了书中早期学习的大多数模型和方法。本章从住房数据的详细示例开始，我们进行了可视化和预处理。主成分方法有助于减少数据，变量聚类方法也有助于完成同样的任务。随后介绍了线性回归模型、神经网络和回归树作为基学习器的方法。袋装、提升和随机森林算法是一些有助于改进模型的方法。这些方法基于同质集成方法。本章最后以堆叠集成方法为三个异质基学习器进行了总结。
下一章将讨论不同结构的有删失观测数据，这种数据被称为生存数据，它通常出现在临床试验的研究中。
第十章：集成生存模型
原始的胆汁性肝硬化数据在前两章中使用了 Jackknife 方法进行介绍。临床试验中的观察结果通常会受到截尾的影响，而 Jackknife 方法通过伪值的概念帮助完成不完整的观察。由于伪值很可能会相互依赖，广义估计方程框架使得在感兴趣的时间点估计十二个协变量的影响成为可能。伪值的概念和广义估计方程框架使得从业者更容易解释结果。然而，如果截尾观察的数量异常高，这种方法可能并不适用。此外，也 preferable to have statistical methods that preserve the incompleteness of the observations and yet make good use of them. 在常规回归框架中，可以将以时间为因变量，误差项遵循适当的生存分布的（线性）回归框架建立起来。然而，事实证明这种方法并不可靠，在许多情况下，它被认为是不可靠的，或者收敛根本不发生。在《回归模型与生命表》（www.stat.cmu.edu/~ryantibs/journalclub/cox_1972.pdf）中，Cox（1972）在提出比例风险模型时，在生存数据的回归建模方面取得了突破。那么，什么是风险模型？
本章将首先介绍核心生存分析概念，如风险率、累积风险函数和生存函数。还将讨论和通过 R 程序可视化一些参数生存模型。对于给定的数据，我们将研究如何通过非参数方法进行寿命分布的推断。然后，对于pbc数据集中的感兴趣事件时间，将展示生存函数和累积风险函数的估计。通过 logrank 检验对不同段落的pbc数据进行假设检验。回归模型将从参数回归模型的一个简单示例开始，以指数分布为例。众所周知，参数模型对于临床试验数据并不非常有效。这导致 Cox 比例风险回归模型的一个重要变体，即半参数模型，在这种模型中，基线风险率被完全未指定，协变量的影响通过风险率上的指数线性项进行建模。
生存树是决策树在生存数据上的一个重要变体，其分裂标准基于logrank检验。自然地，我们将对生存数据的集成方法感兴趣，因此我们在最后一节中开发了生存随机森林。
本章我们将涵盖以下主题：
-
生存分析的基本概念，如风险率、累积风险函数和生存函数
-
Nelson-Aalen 和 Kaplan-Meier 估计量分别作为累积风险函数和生存函数的估计量
-
Logrank 测试用于比较生存曲线
-
参数和半参数方法分析独立协变量对风险率的影响
-
基于 Logrank 测试的生存树
-
随机森林作为生存数据的集成方法
生存分析的核心概念
生存分析处理截尾数据，在临床试验中观察到寿命时，参数模型通常不适用。
用 T 表示生存时间，或感兴趣的事件发生的时间，我们自然有，这是一个连续随机变量。假设寿命累积分布为 F，相关的密度函数为 f。我们将定义进一步分析所需的重要概念。我们将探讨生存函数的概念。
假设 T 是寿命的连续随机变量，相关的累积分布函数为 F。在时间 t 的生存函数是观察到的时间点仍然存活的可能性，它由以下定义：
生存函数可以采取不同的形式。让我们通过每个分布的例子来了解生存函数之间的差异，以获得更清晰的图像。
指数分布：假设一个电子元件的寿命分布遵循指数分布，其速率。那么，其密度函数如下：
累积分布函数如下：
指数分布的均值和方差分别为和
。指数分布的生存函数如下：
指数分布的均值。指数分布由单个参数驱动，并且它还享有一种优雅的性质，称为无记忆性质（参见第六章，Tattar 等人，2016 年）。
伽马分布：如果寿命随机变量的概率密度函数 f 具有以下形式，我们说其遵循速率和形状
的伽马分布：
均值和方差分别为 和
。累积分布函数的闭式形式，因此生存函数也不存在。
威布尔分布：如果一个寿命随机变量遵循具有率 和形状
的威布尔分布，那么其概率密度函数
f 的形式如下：
威布尔分布的累积分布函数如下所示：
生存函数如下：
接下来，我们将定义风险率的概念，它也被称为瞬时故障率。
用 T 表示寿命随机变量，用 F 表示相关的累积分布函数，那么在时间 t 的风险率定义为以下：
估计风险率的问题与密度函数的问题一样困难，因此累积函数的概念将是有用的。
用 T 表示寿命随机变量， 表示相关风险率，那么累积风险函数的定义如下：
这三个量之间存在以下关系：
期望值与生存函数的关系如下：
在下一个 R 程序中，我们将可视化三个概率分布的三个生存量。首先，我们将使用par和mfrow函数设置九个图的图形设备。程序以指数分布为例进行解释。考虑时间区间 0-100，并在程序中创建一个名为Time的数值对象。我们将从使用dexp函数计算Time对象的密度函数值开始。这意味着dexp(Time)将为 0-100 之间的每个点计算密度函数f(t)的值。由于生存函数与累积分布函数通过相关联，而
pexp给出了时间点t处的F值，因此指数分布的生存函数计算为1-pexp()。危险率、密度函数和生存函数通过相关联，并且可以轻松获得。通过使用生存函数的值和以下关系，可以获得累积危险函数：
程序随后将重复用于伽马分布和威布尔分布，并对适当的参数进行更改，如下面的代码所示：
> par(mfrow=c(3,3))
> Time <- seq(0,100,1)
> lambda <- 1/20
> expdens <- dexp(Time,rate=lambda)
> expsurv <- 1-pexp(Time,rate=lambda)
> exphaz <- expdens/expsurv
> expcumhaz <- -log(expsurv)
> plot(Time,exphaz,"l",xlab="Time",ylab="Hazard Rate",ylim=c(0,0.1))
> plot(Time,expcumhaz,"l",xlab="Time",ylab="Cumulative Hazard Function")
> mtext("Exponential Distribution")
> plot(Time,expsurv,"l",xlab="Time",ylab="Survival Function")
>
> # Gamma Distribution
> lambda <- 1/10; k <- 2
> gammadens <- dgamma(Time,rate=lambda,shape=k)
> gammasurv <- 1-pgamma(Time,rate=lambda,shape=k)
> gammahaz <- gammadens/gammasurv
> gammacumhaz <- -log(gammasurv)
> plot(Time,gammahaz,"l",xlab="Time",ylab="Hazard Rate")
> plot(Time,gammacumhaz,"l",xlab="Time",ylab="Cumulative Hazard Function")
> mtext("Gamma Distribution")
> plot(Time,gammasurv,"l",xlab="Time",ylab="Survival Function")
>
> # Weibull Distribution
> lambda <- 25; k <- 2
> Weibulldens <- dweibull(Time,scale=lambda,shape=k)
> Weibullsurv <- 1-pweibull(Time,scale=lambda,shape=k)
> Weibullhaz <- Weibulldens/Weibullsurv
> Weibullcumhaz <- -log(Weibullsurv)
> plot(Time,Weibullhaz,"l",xlab="Time",ylab="Hazard Rate")
> plot(Time,Weibullcumhaz,"l",xlab="Time",ylab="Cumulative Hazard Function")
> mtext("Weibull Distribution")
> plot(Time,Weibullsurv,"l",xlab="Time",ylab="Survival Function")
图 1：指数、伽马和威布尔分布的危脸率、累积危脸函数和生存函数
重复前面的程序，使用不同的参数值，并为危险函数、累积危险函数和生存函数的变化准备观察总结。分别总结指数、伽马和威布尔分布的观察结果。
现在，我们需要查看模型与pbc数据集的拟合程度。在这里，我们将指数、伽马和威布尔分布拟合到pbc数据集中感兴趣的生命周期。请注意，由于我们有截尾数据，不能简单地丢弃不完整的观察值，因为 418 个中有 257 个是不完整的观察值。虽然我们不能深入探讨生存数据的最大似然估计的数学原理，但在此处重要的是要注意，完整观察对似然函数的贡献是f(t)，如果它是不完整/截尾的，则是S(t)。因此，对于软件来说，知道哪个观察值是完整的，哪个是不完整的非常重要。在这里，我们将使用survival包中的Surv函数来指定这一点，然后使用flexsurv包中的flexsurvreg函数来拟合适当的生命周期分布。dist选项有助于设置适当的分布，如下面的程序所示：
> pbc <- survival::pbc
> Surv(pbc$time,pbc$status==2)
[1] 400 4500+ 1012 1925 1504+ 2503 1832+ 2466 2400 51 3762
[12] 304 3577+ 1217 3584 3672+ 769 131 4232+ 1356 3445+ 673
[397] 1328+ 1375+ 1260+ 1223+ 935 943+ 1141+ 1092+ 1150+ 703 1129+
[408] 1086+ 1067+ 1072+ 1119+ 1097+ 989+ 681 1103+ 1055+ 691+ 976+
> pbc_exp <- flexsurvreg(Surv(time,status==2)~1,data=pbc,dist="exponential")
> pbc_exp
Call:
flexsurvreg(formula = Surv(time, status == 2) ~ 1, data = pbc,
dist = "exponential")
Estimates:
est L95% U95% se
rate 2.01e-04 1.72e-04 2.34e-04 1.58e-05
N = 418, Events: 161, Censored: 257
Total time at risk: 801633
Log-likelihood = -1531.593, df = 1
AIC = 3065.187
> windows(height=100,width=100)
> plot(pbc_exp,ylim=c(0,1),col="black")
> pbc_gamma <- flexsurvreg(Surv(time,status==2)~1,data=pbc,dist="gamma")
> pbc_gamma
Call:
flexsurvreg(formula = Surv(time, status == 2) ~ 1, data = pbc,
dist = "gamma")
Estimates:
est L95% U95% se
shape 1.10e+00 9.21e-01 1.30e+00 9.68e-02
rate 2.33e-04 1.70e-04 3.21e-04 3.78e-05
N = 418, Events: 161, Censored: 257
Total time at risk: 801633
Log-likelihood = -1531.074, df = 2
AIC = 3066.147
> plot(pbc_gamma,col="blue",add=TRUE)
> pbc_Weibull <- flexsurvreg(Surv(time,status==2)~1,data=pbc,dist="weibull")
> pbc_Weibull
Call:
flexsurvreg(formula = Surv(time, status == 2) ~ 1, data = pbc,
dist = "weibull")
Estimates:
est L95% U95% se
shape 1.08e+00 9.42e-01 1.24e+00 7.48e-02
scale 4.71e+03 3.96e+03 5.59e+03 4.13e+02
N = 418, Events: 161, Censored: 257
Total time at risk: 801633
Log-likelihood = -1531.017, df = 2
AIC = 3066.035
> plot(pbc_Weibull,col="orange",add=TRUE)
> legend(3000,1,c("Exponential","Gamma","Weibull"),
+ col=c("black","blue","orange"),merge=TRUE,lty=2)
结果图如下所示：
图 2：对截尾数据进行指数、伽马和威布尔分布的拟合
拟合的指数模型的 AIC 值是 3065.187，拟合的伽马模型的 AIC 值是 3066.147，而韦伯尔分布的 AIC 值是 3066.035。标准越低越好。因此，根据 AIC 标准，指数模型是最好的拟合。然后是韦伯尔分布和伽马分布。在这里，单参数的指数分布比更复杂的伽马和韦伯尔模型拟合得更好。
现在需要对 R 程序进行一些解释。pbc数据集通过survival::pbc加载，因为 R 包randomForestSRC也有一个同名数据集，它是稍有不同的版本。因此，survival::pbc代码确保我们继续加载pbc数据集，就像我们在早期实例中所做的那样。对我们来说，感兴趣的事件由status==2和Surv(pbc$time,pbc$status==2)指示，它创建了一个包含在数值对象中的完整观察值的生存对象。如果status不是2，则观察值是删失的，这由跟在+号后面的数字表示。Surv(time,status==2)~1代码创建必要的公式，这对于应用生存函数很有用。dist="exponential"选项确保在生存数据上拟合指数分布。当在控制台上运行拟合模型pbc_exp时，我们得到拟合模型的摘要，并返回模型参数的估计值、95%置信区间和参数估计的标准误差。我们还得到完整和删失观察值的数量、所有患者的总风险时间、似然函数值和 AIC。注意三个拟合分布的自由度是如何变化的。
这里详细描述的参数模型给出了生存概念的一个概念。当我们没有足够的证据来构建一个参数模型时，我们求助于非参数和半参数模型来进行统计推断。在下一节中，我们将继续分析pbc数据。
非参数推断
生存数据会受到删失的影响，我们需要引入一个新的量来捕捉这个信息。假设我们有一个来自的独立同分布的寿命随机变量的n样本，并且我们知道感兴趣的事件可能已经发生，或者它将在未来的某个时间发生。通过克罗内克指示变量
捕捉了额外的信息：
因此，我们在 Ts 和 s，
中有 n 对随机观测值。为了获得累积风险函数和生存函数的估计值，我们需要额外的符号。让
表示 Ts 中事件发生时的唯一时间点。接下来，我们用
表示在时间
和
之前处于风险中的观测值数量，而
表示在该时间发生的事件数量。使用这些量，我们现在提出使用以下方法来估计累积风险函数：
估计量 是著名的 Nelson-Aalen 估计量。Nelson-Aalen 估计量具有统计性质，包括（i）它是累积风险函数的非参数最大似然估计量，以及（ii）它遵循渐近正态分布。生存函数的估计量如下：
估计量 是著名的 Kaplan-Meier 估计量。通过应用函数-Δ定理，Nelson-Aalen 估计量的性质被传递到 Kaplan-Meier 估计量。需要注意的是，Kaplan-Meier 估计量再次是非参数最大似然估计量，并且渐近地遵循正态分布。我们现在将探讨如何使用 R 软件获取给定数据集的估计值。
我们已经使用 Surv(pbc$time, pbc$status==2) 代码创建了生存对象。现在，在生存对象上应用 survfit 函数，我们在 pbc_sf survfit 对象中设置了 Kaplan-Meier 估计器：
> pbc_sf <- survfit(Surv(time,status==2)~1,pbc)
> pbc_sf
Call: survfit(formula = Surv(time, status == 2) ~ 1, data = pbc)
n events median 0.95LCL 0.95UCL
418 161 3395 3090 3853
输出结果显示我们共有 418 个观测值。其中，161 个经历了感兴趣的事件。我们希望获取不同时间点的生存函数。中位生存时间为 3395，该点估计值的置信区间为 3090 至 3853。然而，如果你找到整体时间的平均值，完整观测值的平均时间，以及截尾观测值的平均时间，这些值都不会接近显示的 3395 值。一段快速代码显示了以下结果：
> median(pbc$time)
[1] 1730
> median(pbc$time[pbc$status==2])
[1] 1083
> median(pbc$time[pbc$status!=2])
[1] 2157
你可能会问自己，为什么估计的中位生存时间与这些中位数之间有如此大的差异？答案很快就会变得清晰。
我们将使用 summary 函数来获取这个答案。对于观察到的十个时间分位数，包括截尾时间，我们将获得 Kaplan-Meier 估计值及其相关的 95% 置信区间，该置信区间基于方差估计，而方差估计又基于 Greenwood 公式：
> summary(pbc_sf,times=as.numeric(quantile(pbc$time,seq(0,1,0.1))))
Call: survfit(formula = Surv(time, status == 2) ~ 1, data = pbc)
time n.risk n.event survival std.err lower 95% CI upper 95% CI
41 418 2 0.995 0.00338 0.989 1.000
607 376 39 0.902 0.01455 0.874 0.931
975 334 31 0.827 0.01860 0.791 0.864
1218 292 19 0.778 0.02061 0.738 0.819
1435 251 8 0.755 0.02155 0.714 0.798
1730 209 13 0.713 0.02323 0.669 0.760
2107 167 12 0.668 0.02514 0.621 0.719
2465 126 9 0.628 0.02702 0.577 0.683
2852 84 10 0.569 0.03032 0.512 0.632
3524 42 10 0.478 0.03680 0.411 0.556
4795 1 8 0.353 0.04876 0.270 0.463
我们现在已经得到了每个十分位时间点的 Kaplan-Meier 估计值、每个点的标准误差和置信区间。使用plot函数，我们现在将可视化pbc数据集的拟合 Kaplan-Meier 估计值：
> plot(pbc_sf,xlab="Time",ylab="Survival Function Confidence Bands")
图 3：PBC 数据集的 Kaplan-Meier 估计值
现在，如果你观察生存时间接近 0.5 的时间点，中位生存时间 3395 的早期答案就足够清晰了。接下来，我们来看累积风险函数。
要获得累积风险函数，我们将对生存对象应用coxph函数，并使用basehaz函数获取基线累积风险函数，如下面的代码所示：
> pbc_na <- basehaz(coxph(Surv(time,status==2)~1,pbc))
> pbc_na
hazard time
1 0.004790426 41
2 0.007194272 43
3 0.009603911 51
4 0.012019370 71
396 1.030767970 4509
397 1.030767970 4523
398 1.030767970 4556
399 1.030767970 4795
我们将使用以下代码创建 Nelson-Aalen 估计值的可视化表示：
> plot(pbc_na$time,pbc_na$hazard,"l",xlab="Time",ylab="Cumulative Hazard Function")
下面的图表说明了 Nelson-Aalen 估计值：
图 4：累积风险函数的 Nelson-Aalen 估计值
注意，有人可能会倾向于使用关系来获得 Kaplan-Meier 估计值，或者反之亦然。让我们使用以下代码来检查一下：
> str(exp(-pbc_na$hazard))
num [1:399] 0.995 0.993 0.99 0.988 0.986 ...
> str(summary(pbc_sf,times=pbc_na$time)$surv)
num [1:399] 0.995 0.993 0.99 0.988 0.986 ...
这里似乎一切正常，所以让我们全面检查一下：
> round(exp(-pbc_na$hazard),4)==round(summary(pbc_sf,
+ times=pbc_na$time)$surv,4)
[1] TRUE TRUE TRUE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE
[12] TRUE TRUE FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE FALSE
[23] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
[34] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[375] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[386] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[397] FALSE FALSE FALSE
经过一段时间后，估计值会有很大差异，因此我们分别计算这两个量。接下来，我们将探讨如何进行统计检验以比较生存曲线的相等性。
如前所述，我们将公式输入到survfit函数中。它出现在Surv公式中的额外指定'~1'。由于公式对于生存数据的进一步分析至关重要，我们现在可以很好地利用这个框架。如果我们用 1 替换为一个分类变量，如性别，那么我们将为分类变量的每个水平获得生存曲线。对于pbc数据，我们将绘制生存曲线。分别对男性和女性的 Kaplan-Meier 估计值进行绘图。
>plot(survfit(Surv(time,status==2)~sex,pbc),conf.int=TRUE,xlab="Time",+ ylab="Survival Probability", col=c("red","blue"))
图 5：PBC 数据的性别生存曲线比较
生存曲线（用蓝色和红色连续线表示）清楚地显示了差异，我们需要评估观察到的差异是否具有统计学意义。为此，我们应用survdiff函数并检查差异是否显著，如下面的代码所示：
> survdiff(Surv(time,status==2)~sex,pbc)
Call:
survdiff(formula = Surv(time, status == 2) ~ sex, data = pbc)
N Observed Expected (O-E)²/E (O-E)²/V
sex=m 44 24 17.3 2.640 2.98
sex=f 374 137 143.7 0.317 2.98
Chisq= 3 on 1 degrees of freedom, p= 0.0845
p 值为0.0845，因此如果选择的显著性水平是 95%，我们得出结论，差异不显著。
读者注意：显著性水平是预先确定的。如果你在进行分析之前将其固定在 95%，然后查看 p 值发现它在 0.05 和 0.10 之间，不要改变水平。坚持之前达成的协议。
在到目前为止的分析中，我们考虑了参数和非参数方法，现在我们需要开发一个更大的框架。需要明确评估协变量的影响，我们将在下一节探讨这个主题。
回归模型 – 参数和 Cox 比例风险模型
你可能还记得，生存数据包括完整的以及截尾的观测值，我们看到了对于 pbc 数据集，寿命看起来像是 400，4500+，1012，1925，1504+，……。尽管寿命是连续的随机变量，但形式为 的回归模型在这里并不合适。事实上，在 20 世纪 70 年代，有许多尝试修正和改进这种形式的模型，而且结果往往是有害的。我们将一个通用的风险回归模型定义为以下形式：
在这里，t 是寿命，
是协变量向量，
是回归系数向量，
是基线风险率。一个特别感兴趣的相对风险模型如下：
我们将专注于这一类模型。首先，考虑参数风险回归。这意味着我们将通过参数模型指定风险率 ，例如，通过指数分布。但这意味着什么呢？这意味着基线风险函数具有以下形式：
因此，风险回归模型如下：
然后估计问题就是找到 和
。
flexsurv 包中的 R 函数 survreg 将有助于拟合参数风险回归模型。它将在 pbc 数据集上展示连续性。survival 公式将被扩展以包括模型中的所有协变量，如下面的代码所示：
> pbc_Exp <- survreg(Surv(time,status==2)~trt + age + sex + ascites
+ hepato + spiders + edema + bili + chol + albumin
+ copper + alk.phos + ast + trig + platelet +
+ protime + stage,
+ dist="exponential", pbc)
> pbc_Exp_summary <- summary(pbc_Exp)
> pbc_Exp_summary
Call:
survreg(formula = Surv(time, status == 2) ~ trt + age + sex +
ascites + hepato + spiders + edema + bili + chol + albumin +
copper + alk.phos + ast + trig + platelet + protime + stage,
data = pbc, dist = "exponential")
Value Std. Error z p
(Intercept) 1.33e+01 1.90e+00 7.006 2.46e-12
trt 6.36e-02 2.08e-01 0.305 7.60e-01
age -2.81e-02 1.13e-02 -2.486 1.29e-02
sexf 3.42e-01 3.00e-01 1.140 2.54e-01
ascites -1.01e-01 3.70e-01 -0.273 7.85e-01
hepato -7.76e-02 2.43e-01 -0.320 7.49e-01
spiders -1.21e-01 2.36e-01 -0.513 6.08e-01
edema -8.06e-01 3.78e-01 -2.130 3.32e-02
bili -4.80e-02 2.49e-02 -1.929 5.37e-02
chol -4.64e-04 4.35e-04 -1.067 2.86e-01
albumin 4.09e-01 2.87e-01 1.427 1.54e-01
copper -1.63e-03 1.11e-03 -1.466 1.43e-01
alk.phos -4.51e-05 3.81e-05 -1.182 2.37e-01
ast -3.68e-03 1.92e-03 -1.917 5.52e-02
trig 1.70e-04 1.37e-03 0.124 9.01e-01
platelet -2.02e-04 1.16e-03 -0.174 8.62e-01
protime -2.53e-01 9.78e-02 -2.589 9.61e-03
stage -3.66e-01 1.66e-01 -2.204 2.75e-02
Scale fixed at 1
Exponential distribution
Loglik(model)= -984.1 Loglik(intercept only)= -1054.6
Chisq= 141.17 on 17 degrees of freedom, p= 0
Number of Newton-Raphson Iterations: 5
n=276 (142 observations deleted due to missingness)
在这里，迭代五次后收敛。p 值几乎等于零，这意味着拟合的模型是显著的。然而，并非所有与协变量相关的 p 值都表明显著性。我们将使用以下代码来查找：
> round(pbc_Exp_summary$table[,4],4)
(Intercept) trt age sexf ascites hepato
0.0000 0.7601 0.0129 0.2541 0.7846 0.7491
spiders edema bili chol albumin copper
0.6083 0.0332 0.0537 0.2859 0.1536 0.1426
alk.phos ast trig platelet protime stage
0.2372 0.0552 0.9010 0.8621 0.0096 0.0275
> AIC(pbc_exp)
[1] 3065.187
AIC 值也非常高，我们试图看看是否可以改进。因此，我们将step函数应用于拟合的指数风险回归模型，并消除不显著的协变量，如下面的代码所示：
> pbc_Exp_eff <- step(pbc_Exp)
Start: AIC=2004.12
Surv(time, status == 2) ~ trt + age + sex + ascites + hepato +
spiders + edema + bili + chol + albumin + copper + alk.phos +
ast + trig + platelet + protime + stage
Df AIC
- ascites 1 2002.2
- trt 1 2002.2
- hepato 1 2002.2
- spiders 1 2002.4
- chol 1 2003.2
- sex 1 2003.4
- alk.phos 1 2003.4
- albumin 1 2004.1
<none> 2004.1
- ast 1 2005.5
- bili 1 2005.5
- edema 1 2006.5
- stage 1 2007.2
- protime 1 2008.0
- age 1 2008.3
- trig 1 2020.4
- copper 1 2020.7
- platelet 1 2021.8
Step: AIC=2002.19
Surv(time, status == 2) ~ trt + age + sex + hepato + spiders +
edema + bili + chol + albumin + copper + alk.phos + ast +
trig + platelet + protime + stage
Step: AIC=1994.61
Surv(time, status == 2) ~ age + edema + bili + albumin + copper +
ast + trig + platelet + protime + stage
Df AIC
<none> 1994.6
- albumin 1 1995.5
- edema 1 1996.3
- ast 1 1996.9
- bili 1 1998.8
- protime 1 1999.2
- stage 1 1999.7
- age 1 2000.9
- platelet 1 2012.1
- copper 1 2014.9
- trig 1 2198.7
There were 50 or more warnings (use warnings() to see the first 50)
> pbc_Exp_eff_summary <- summary(pbc_Exp_eff)
> round(pbc_Exp_eff_summary$table[,4],4)
(Intercept) age edema bili albumin copper
0.0000 0.0037 0.0507 0.0077 0.0849 0.0170
ast trig platelet protime stage
0.0281 0.9489 0.7521 0.0055 0.0100
> AIC(pbc_Exp_eff)
[1] 1994.607
我们在这里看到，当前模型中的所有协变量都是显著的，除了trig和platelet变量。AIC 值也急剧下降。
在生命科学中，参数模型通常不可接受。的灵活框架是著名的 Cox 比例风险模型。它是一个半参数回归模型，因为基线风险函数
是完全未指定的。Cox（1972）提出了这个模型，这是统计学中最重要的模型之一。基线风险函数
的唯一要求是它必须是非负的，并且与之相关的概率分布必须是合适的概率分布。在这个模型中，回归系数向量
通过将基线风险函数视为一个干扰因素来估计。其推断基于部分似然函数的重要概念；有关完整细节，请参阅 Cox（1975）。在这里，我们只指定 Cox 比例风险模型的形式，并将感兴趣的读者指引到 Kalbfleisch 和 Prentice（2002）：
我们将使用survival包中的coxph函数来拟合比例风险回归模型。由于某些技术原因，我们必须从pbc数据集中省略所有包含缺失值的行，并且剩余步骤与指数风险回归模型的拟合平行：
> pbc2 <- na.omit(pbc)
> pbc_coxph <- coxph(Surv(time,status==2)~trt + age + sex + ascites
+ hepato + spiders + edema + bili + chol + albumin
+ copper + alk.phos + ast + trig + platelet +
+ protime + stage, pbc2)
> pbc_coxph_summary <- summary(pbc_coxph)
> pbc_coxph_summary
Call:
coxph(formula = Surv(time, status == 2) ~ trt + age + sex + ascites +
hepato + spiders + edema + bili + chol + albumin + copper +
alk.phos + ast + trig + platelet + protime + stage, data = pbc2)
n= 276, number of events= 111
coef exp(coef) se(coef) z Pr(>|z|)
trt -1.242e-01 8.832e-01 2.147e-01 -0.579 0.56290
age 2.890e-02 1.029e+00 1.164e-02 2.482 0.01305 *
sexf -3.656e-01 6.938e-01 3.113e-01 -1.174 0.24022
ascites 8.833e-02 1.092e+00 3.872e-01 0.228 0.81955
hepato 2.552e-02 1.026e+00 2.510e-01 0.102 0.91900
spiders 1.012e-01 1.107e+00 2.435e-01 0.416 0.67760
edema 1.011e+00 2.749e+00 3.941e-01 2.566 0.01029 *
bili 8.001e-02 1.083e+00 2.550e-02 3.138 0.00170 **
chol 4.918e-04 1.000e+00 4.442e-04 1.107 0.26829
albumin -7.408e-01 4.767e-01 3.078e-01 -2.407 0.01608 *
copper 2.490e-03 1.002e+00 1.170e-03 2.128 0.03337 *
alk.phos 1.048e-06 1.000e+00 3.969e-05 0.026 0.97893
ast 4.070e-03 1.004e+00 1.958e-03 2.078 0.03767 *
trig -9.758e-04 9.990e-01 1.333e-03 -0.732 0.46414
platelet 9.019e-04 1.001e+00 1.184e-03 0.762 0.44629
protime 2.324e-01 1.262e+00 1.061e-01 2.190 0.02850 *
stage 4.545e-01 1.575e+00 1.754e-01 2.591 0.00958 ** ---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
trt 0.8832 1.1323 0.5798 1.3453
age 1.0293 0.9715 1.0061 1.0531
sexf 0.6938 1.4414 0.3769 1.2771
ascites 1.0924 0.9155 0.5114 2.3332
hepato 1.0259 0.9748 0.6273 1.6777
spiders 1.1066 0.9037 0.6865 1.7835
edema 2.7487 0.3638 1.2697 5.9505
bili 1.0833 0.9231 1.0305 1.1388
chol 1.0005 0.9995 0.9996 1.0014
albumin 0.4767 2.0977 0.2608 0.8714
copper 1.0025 0.9975 1.0002 1.0048
alk.phos 1.0000 1.0000 0.9999 1.0001
ast 1.0041 0.9959 1.0002 1.0079
trig 0.9990 1.0010 0.9964 1.0016
platelet 1.0009 0.9991 0.9986 1.0032
protime 1.2617 0.7926 1.0247 1.5534
stage 1.5754 0.6348 1.1170 2.2219
Concordance= 0.849 (se = 0.031 )
Rsquare= 0.455 (max possible= 0.981 )
Likelihood ratio test= 167.7 on 17 df, p=0
Wald test = 174.1 on 17 df, p=0
Score (logrank) test = 283.7 on 17 df, p=0
> round(pbc_coxph_summary$coefficients[,5],4)
trt age sexf ascites hepato spiders edema bili
0.5629 0.0131 0.2402 0.8195 0.9190 0.6776 0.0103 0.0017
chol albumin copper alk.phos ast trig platelet protime
0.2683 0.0161 0.0334 0.9789 0.0377 0.4641 0.4463 0.0285
stage
0.0096
> AIC(pbc_coxph)
[1] 966.6642
由于我们发现许多变量不显著，我们将尝试通过使用step函数并计算改进的 AIC 值来改进它，如下面的代码所示：
> pbc_coxph_eff <- step(pbc_coxph)
Start: AIC=966.66
Surv(time, status == 2) ~ trt + age + sex + ascites + hepato +
spiders + edema + bili + chol + albumin + copper + alk.phos +
ast + trig + platelet + protime + stage
Df AIC
- alk.phos 1 964.66
- hepato 1 964.67
- ascites 1 964.72
- spiders 1 964.84
- trt 1 965.00
- trig 1 965.22
- platelet 1 965.24
- chol 1 965.82
- sex 1 965.99
<none> 966.66
- ast 1 968.69
- copper 1 968.85
- protime 1 968.99
- albumin 1 970.35
- age 1 970.84
- edema 1 971.00
- stage 1 971.83
- bili 1 973.34
Step: AIC=952.58
Surv(time, status == 2) ~ age + edema + bili + albumin + copper +
ast + protime + stage
Df AIC
<none> 952.58
- protime 1 955.06
- ast 1 955.79
- edema 1 955.95
- albumin 1 957.27
- copper 1 958.18
- age 1 959.97
- stage 1 960.11
- bili 1 966.57
> pbc_coxph_eff_summary <- summary(pbc_coxph_eff)
> pbc_coxph_eff_summary
Call:
coxph(formula = Surv(time, status == 2) ~ age + edema + bili +
albumin + copper + ast + protime + stage, data = pbc2)
n= 276, number of events= 111
coef exp(coef) se(coef) z Pr(>|z|)
age 0.0313836 1.0318812 0.0102036 3.076 0.00210 **
edema 0.8217952 2.2745795 0.3471465 2.367 0.01792 *
bili 0.0851214 1.0888492 0.0193352 4.402 1.07e-05 ***
albumin -0.7185954 0.4874364 0.2724486 -2.638 0.00835 **
copper 0.0028535 1.0028576 0.0009832 2.902 0.00370 **
ast 0.0043769 1.0043865 0.0018067 2.423 0.01541 *
protime 0.2275175 1.2554794 0.1013729 2.244 0.02481 *
stage 0.4327939 1.5415584 0.1456307 2.972 0.00296 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
age 1.0319 0.9691 1.0114 1.0527
edema 2.2746 0.4396 1.1519 4.4915
bili 1.0888 0.9184 1.0484 1.1309
albumin 0.4874 2.0515 0.2858 0.8314
copper 1.0029 0.9972 1.0009 1.0048
ast 1.0044 0.9956 1.0008 1.0079
protime 1.2555 0.7965 1.0292 1.5314
stage 1.5416 0.6487 1.1588 2.0508
Concordance= 0.845 (se = 0.031 )
Rsquare= 0.448 (max possible= 0.981 )
Likelihood ratio test= 163.8 on 8 df, p=0
Wald test = 176.1 on 8 df, p=0
Score (logrank) test = 257.5 on 8 df, p=0
> round(pbc_coxph_eff_summary$coefficients[,5],4)
age edema bili albumin copper ast protime stage
0.0021 0.0179 0.0000 0.0084 0.0037 0.0154 0.0248 0.0030
> AIC(pbc_coxph_eff)
[1] 952.5814
我们现在可以看到，pbc_coxph_eff模型中的几乎所有变量都是显著的。AIC 值也与其早期值相比有所下降。
在大多数生存数据分析中，目的是找到有效的协变量及其对风险率的影响。在生存数据的情况下，不存在像 AUC 这样的分类问题准确度度量。在类似（尽管重要）的方面，对于给定协变量选择的寿命预测可能不再重要。在第二章中，我们探讨了伪值方法和 Jackknife 方法的应用，这些方法提供了易于解释的便利。在下一节中，我们将探讨不同的方法。
生存树
参数风险回归模型有时被从业者视为一种限制性模型类别，Cox 比例风险回归有时比其参数对应物更受欢迎。与参数模型相比，解释有时会丢失，常规从业者发现很难与风险回归模型联系起来。当然，一个替代方案是在伪观察值上构建生存树。这种尝试可以在 Tattar（2016）未发表的论文中看到。Gordon 和 Olshen（1985）首次尝试构建生存树，许多科学家继续构建它。LeBlanc 和 Crowley（1992）是建立生存树的重要贡献者之一。Zhang 和 Singer（2010）也对相关方法进行了系统发展，他们书中第 7-10 章讨论了生存树。基本前提保持不变，我们需要好的分割标准来创建生存树。
LeBlanc 和 Crowley 提出了基于饱和模型对数似然与最大对数似然之间的节点偏差度量的分割标准，然后通过用 Nelson-Aalen 估计器估计的基线累积风险函数替换，来近似未知的完整似然。正如 Hamad 等人（2011）所建议的，请注意，这种方法的优势在于任何能够执行泊松树实现的软件都能够创建生存树。这种方法在 Therneau 和 Atkinson（2017）创建的rpartR 包中被利用。生存树的终端节点可以通过属于该节点的观察值的 Kaplan-Meier 生存函数来总结。
我们现在将使用rpart库为pbc数据集设置生存树：
> pbc_stree <- rpart(Surv(time,status==2)~trt + age + sex + ascites +
+ hepato + spiders + edema + bili + chol + albumin +
+ copper + alk.phos + ast + trig + platelet +
+ protime + stage,
+ pbc)
> pbc_stree
n= 418
node), split, n, deviance, yval
* denotes terminal node
1) root 418 555.680700 1.0000000
2) bili< 2.25 269 232.626500 0.4872828
4) age< 51.24162 133 76.376990 0.2395736
8) alk.phos< 1776 103 30.340440 0.1036750 *
9) alk.phos>=1776 30 33.092480 0.6135695 *
5) age>=51.24162 136 136.800300 0.7744285
10) protime< 10.85 98 80.889260 0.5121123 *
11) protime>=10.85 38 43.381670 1.4335430
22) age< 65.38125 26 24.045820 0.9480269
44) bili< 0.75 8 5.188547 0.3149747 *
45) bili>=0.75 18 12.549130 1.3803650 *
23) age>=65.38125 12 8.392462 3.2681510 *
3) bili>=2.25 149 206.521900 2.6972690
6) protime< 11.25 94 98.798830 1.7717830
12) stage< 3.5 57 56.734150 1.2620350
24) age< 43.5332 25 16.656000 0.6044998 *
25) age>=43.5332 32 32.986760 1.7985470 *
13) stage>=3.5 37 32.946240 2.8313470 *
7) protime>=11.25 55 76.597760 5.1836880
14) ascites< 0.5 41 52.276360 4.1601690
28) age< 42.68172 7 6.829564 1.4344660 *
29) age>=42.68172 34 37.566600 5.1138380 *
15) ascites>=0.5 14 17.013010 7.9062910 *
通过在控制台运行pbc_stree，显示生存树的文本表示。节点yval末尾的星号（*）表示删失。
我们现在将查看与生存树相关的cptable和变量重要性，如下程序所示：
> pbc_stree$cptable
CP nsplit rel error xerror xstd
1 0.20971086 0 1.0000000 1.0026398 0.04761402
2 0.05601298 1 0.7902891 0.8179390 0.04469359
3 0.03500063 2 0.7342762 0.8198147 0.04731189
4 0.02329408 3 0.6992755 0.8247148 0.04876085
5 0.02254786 4 0.6759815 0.8325230 0.05132829
6 0.01969366 5 0.6534336 0.8363125 0.05160507
7 0.01640950 6 0.6337399 0.8373375 0.05268262
8 0.01366665 7 0.6173304 0.8336743 0.05406099
9 0.01276163 9 0.5899971 0.8209714 0.05587843
10 0.01135209 10 0.5772355 0.8248827 0.05612403
11 0.01000000 11 0.5658834 0.8255873 0.05599763
> pbc_stree$variable.importance
bili protime age albumin edema stage
138.007841 70.867308 54.548224 32.239919 25.576170 15.231256
ascites alk.phos platelet sex copper ast
14.094208 13.440866 10.017966 2.452776 2.114888 1.691910
变量重要性显示，bili变量是最重要的，其次是protime和age。我们用以下代码的视觉表示来结束本节：
> windows(height=100,width=60)
> plot(pbc_stree,uniform = TRUE)
> text(pbc_stree,use.n=TRUE)
图 6：PBC 数据集的生存树
就像单个树模型中存在的过拟合问题一样，我们需要考虑生存树集成的重要替代方法。
集成生存模型
随机森林包randomForestSRC将继续用于创建与生存数据相关的随机森林。实际上，包名中的SRC的S代表生存！rfsrc函数的使用方法与前面章节相同，现在我们将提供一个Surv对象，如下面的代码所示：
> pbc_rf <- rfsrc(Surv(time,status==2)~trt + age + sex + ascites +
+ hepato + spiders + edema + bili + chol + albumin+
+ copper + alk.phos + ast +trig + platelet + protime+
+ stage, ntree=500, tree.err = TRUE, pbc)
我们将找到设置此随机森林的一些基本设置：
> pbc_rf$splitrule
[1] "logrankCR"
> pbc_rf$nodesize
[1] 6
> pbc_rf$mtry
[1] 5
因此，分裂标准基于 log-rank 检验，终端节点中的最小观测数是六个，每个分裂考虑的变量数是五个。
接下来，我们将找到两个事件的变量重要性，并应用我们在前面章节中使用的var.select函数。然后，我们将展示迭代运行的一部分，如下面的代码所示：
> vimp(pbc_rf)$importance
event.1 event.2
trt 0.0013257796 0.0002109143
age 0.0348848966 0.0166352983
sex 0.0007755091 0.0004011303
ascites 0.0008513276 0.0107212361
hepato 0.0050666763 0.0015445001
spiders -0.0001136547 0.0003552415
edema 0.0006227470 0.0147982184
bili 0.0696654202 0.0709713627
chol 0.0002483833 0.0107024051
albumin -0.0106392917 0.0115188264
copper 0.0185417386 0.0255099568
alk.phos 0.0041407841 0.0022297323
ast 0.0029317937 0.0063469825
trig -0.0040190463 0.0014371745
platelet 0.0021021396 0.0002388797
protime 0.0057968358 0.0133710339
stage -0.0023944666 0.0042808034
> var.select(pbc_rf, method = "vh.vimp", nrep = 50)
--------------------- Iteration: 1 ---------------------
selecting variables using Variable Hunting (VIMP) ...
PE: 38.8889 dim: 2
--------------------- Iteration: 2 ---------------------
selecting variables using Variable Hunting (VIMP) ...
PE: 23.6842 dim: 2
--------------------- Iteration: 3 ---------------------
selecting variables using Variable Hunting (VIMP) ...
PE: 50.6667 dim: 1
--------------------- Iteration: 48 ---------------------
selecting variables using Variable Hunting (VIMP) ...
PE: 38.3562 dim: 2
--------------------- Iteration: 49 ---------------------
selecting variables using Variable Hunting (VIMP) ...
PE: 19.4737 dim: 2
--------------------- Iteration: 50 ---------------------
selecting variables using Variable Hunting (VIMP) ...
PE: 55.914 dim: 2
fitting forests to final selected variables ...
-----------------------------------------------------------
family : surv-CR
var. selection : Variable Hunting (VIMP)
conservativeness : medium
dimension : 17
sample size : 276
K-fold : 5
no. reps : 50
nstep : 1
ntree : 500
nsplit : 10
mvars : 4
nodesize : 2
refitted forest : TRUE
model size : 2.16 +/- 0.6809
PE (K-fold) : 43.3402 +/- 17.7472
Top variables:
rel.freq
bili 28
hepato 24
ast 24
-----------------------------------------------------------
接下来，随着树的数量增加，错误率图将使用以下代码进行展示：
> plot(pbc_rf,plots.one.page = TRUE)
图 7：两个事件的随机森林错误率
因此，我们现在已经看到了生存数据的随机森林构建和设置。
摘要
生存数据与典型的回归数据不同，不完整的观测值构成了挑战。由于数据结构完全不同，我们需要专门的技巧来处理不完整的观测值，为此，我们介绍了核心生存概念，如风险率和生存函数。然后，我们介绍了参数生存模型，这让我们对寿命分布的形状有一个大致的了解。我们甚至将这些寿命分布拟合到 pbc 数据集中。
我们还了解到参数设置可能非常严格，因此考虑了生存量估计的非参数方法。我们还展示了 Nelson-Aalen 估计量、Kaplan-Meier 生存函数和对数秩检验的实用性。参数风险回归模型以 Cox 比例风险回归模型为支持，并应用于 pbc 数据集。logrank 检验也可以帮助确定分裂标准，它也被视为随机森林的标准。生存树在前面章节中已展示。
在下一章中，我们将探讨另一种数据结构：时间序列数据。读者不需要熟悉时间序列分析，就可以理解其中应用的集成方法。
第十一章. 时间序列模型的集成
本书迄今为止开发的模型都处理了观测值相互独立的情况。海外游客的例子解释了一个时间序列，其中观测值依赖于先前观测到的数据。在简要讨论这个例子时，我们确定有必要开发时间序列模型。由于时间序列具有序列性质，时间戳可能以纳秒、秒、分钟、小时、天或月为单位显示。
本章将首先简要回顾时间序列自相关和偏自相关函数的重要概念，以及拟合模型评估指标。与分类和回归模型类似，有许多方法可用于分析时间序列数据。季节分解中的时间序列模型的一个重要类别包括 LOESS（STL）、指数平滑状态空间模型（ets）、Box-Jenkins（ARIMA）模型和自回归神经网络模型。这些将在下一节中进行讨论和说明。本章的最后部分将展示时间序列模型的集成。
本章将涵盖以下主要内容：
-
时间序列数据集
-
时间序列可视化
-
核心概念
-
时间序列模型
-
时间序列的 Bagging
-
集成时间序列模型
技术要求
在本章中，我们将使用以下 R 库：
-
forecast -
forecastHybrid
时间序列数据集
时间序列数据在结构上与之前讨论的数据不同。在第一章的海外游客部分 1 中，以及第一章的集成技术简介中，我们看到了时间序列数据的一瞥。在第二章的重抽样中，时间序列模型的重抽样被简要提及。分析时间序列数据时产生的复杂性在于，观测值不是独立的，因此我们需要指定依赖关系。Box 等人（2015 年）的《时间序列统计分析基准书》首次于 1970 年出版。Box 和 Jenkins 发明并普及的模型类别是流行的自回归积分移动平均，通常缩写为 ARIMA。这也常被称为 Box-Jenkins 模型。
表 1总结了二十一个时间序列数据集。长度列给出了序列的观测值/数据点的数量，而频率列给出了时间序列的周期性，其余六列是通过对数值对象应用汇总函数得到的总结。第一列，当然，给出了数据集在 R 中的名称，因此我们没有改变大小写。数据集中的观测值数量从 19 到 7980 不等。
但频率或周期性究竟是什么意思呢？在一个包含周期性的数据集中，相关的时间索引会被重复。例如，我们可能会有年度、季度、月度或周度数据，因此，在后两种情况下，频率将分别是4和12。频率不一定是整数，也可以是分数值。例如，致癌性测试将具有纳秒值。时间序列数据的总结就是将汇总函数应用于数值数据集的结果。因此，我们隐含地假设时间序列数据是数值的。总结中看到的变异也表明，不同的数据集将需要不同的处理。Tattar 等人（2017）的第十章中可以找到时间序列应用的简要介绍。
本章所使用的数据集的描述可以在下表中找到：
| Dataset | Length | Frequency | Minimum | Q1 | Median | Mean | Q3 | Maximum |
|---|---|---|---|---|---|---|---|---|
AirPassengers | 144 | 12 | 104.00 | 180.00 | 265.50 | 280.30 | 360.50 | 622.00 |
BJsales | 150 | 1 | 198.60 | 212.58 | 220.65 | 229.98 | 254.68 | 263.30 |
JohnsonJohnson | 84 | 4 | 0.44 | 1.25 | 3.51 | 4.80 | 7.13 | 16.20 |
LakeHuron | 98 | 1 | 575.96 | 578.14 | 579.12 | 579.00 | 579.88 | 581.86 |
Nile | 100 | 1 | 456.00 | 798.50 | 893.50 | 919.35 | 1032.50 | 1370.00 |
UKgas | 108 | 4 | 84.80 | 153.30 | 220.90 | 337.63 | 469.90 | 1163.90 |
UKDriverDeaths | 192 | 12 | 1057.00 | 1461.75 | 1631.00 | 1670.31 | 1850.75 | 2654.00 |
USAccDeaths | 72 | 12 | 6892.00 | 8089.00 | 8728.50 | 8788.79 | 9323.25 | 11317.00 |
WWWusage | 100 | 1 | 83.00 | 99.00 | 138.50 | 137.08 | 167.50 | 228.00 |
airmiles | 24 | 1 | 412.00 | 1580.00 | 6431.00 | 10527.83 | 17531.50 | 30514.00 |
austres | 89 | 4 | 13067.30 | 14110.10 | 15184.20 | 15273.45 | 16398.90 | 17661.50 |
co2 | 468 | 12 | 313.18 | 323.53 | 335.17 | 337.05 | 350.26 | 366.84 |
discoveries | 100 | 1 | 0.00 | 2.00 | 3.00 | 3.10 | 4.00 | 12.00 |
lynx | 114 | 1 | 39.00 | 348.25 | 771.00 | 1538.02 | 2566.75 | 6991.00 |
nhtemp | 60 | 1 | 47.90 | 50.58 | 51.20 | 51.16 | 51.90 | 54.60 |
nottem | 240 | 12 | 31.30 | 41.55 | 47.35 | 49.04 | 57.00 | 66.50 |
presidents | 120 | 4 | 23.00 | 46.00 | 59.00 | 56.31 | 69.00 | 87.00 |
treering | 7980 | 1 | 0.00 | 0.84 | 1.03 | 1.00 | 1.20 | 1.91 |
气体 | 476 | 12 | 1646.00 | 2674.75 | 16787.50 | 21415.27 | 38628.50 | 66600.00 |
uspop | 19 | 0.1 | 3.93 | 15.00 | 50.20 | 69.77 | 114.25 | 203.20 |
太阳黑子 | 2820 | 12 | 0.00 | 15.70 | 42.00 | 51.27 | 74.93 | 253.80 |
表 1：R 中的时间序列数据集
AirPassengers
AirPassengers数据集包含 1949 年至 1960 年每月的国际航空公司乘客总数。每月计数数字以千为单位。在十二年的时间里，每月数据累积了 144 个观测值。由于我们在多年的月份中有多于一个的观测值，因此可以从数据中捕捉到旅客计数的季节性方面。这在 Box 等人（2015）中得到了普及。
二氧化碳
二氧化碳时间序列数据与大气中二氧化碳的浓度相关。浓度以百万分之一（ppm）表示，该数据集报告在 1997 年 SIO 压力计摩尔分数初步尺度上。这个时间序列在 1959-97 年期间按月捕捉。在co2的帮助页面上，指出 1964 年 2 月、3 月和 4 月的缺失值是通过在 1964 年 1 月和 5 月之间的值进行线性插值获得的。
uspop
美国人口普查（以百万为单位），uspop，在 1790 年至 1970 年间的十年人口普查中被记录。这些数据以小时间序列数据集的形式提供，因此它只包含 19 个数据点。该数据集未捕捉到季节性。
气体
气体时间序列数据包含澳大利亚每月的天然气产量。这里提供的数据是 1956-95 年期间的数据。因此，我们有 476 个观测值。这个数据集来自forecast包。
汽车销售
汽车销售数据来自亚伯拉罕和莱德尔特（1983）。更多信息，请参阅书中第 68 页的表 2.7。销售和广告数据可用期为 36 个月。在这里，我们还有每周广告支出金额的额外信息。这是额外变量可用性的第一个实例，需要专门处理，我们将在本章后面进一步探讨。数据可在Car_Sales.csv文件中找到，并在代码包中提供。
austres
austres时间序列数据集包括 1971 年 3 月至 1994 年 3 月期间澳大利亚居民的季度数据。
WWW 使用情况
WWWusage时间序列数据集包括通过服务器连接到互联网的用户数量。数据以一分钟的时间间隔收集。时间序列值收集了 100 个观测值。
可视化提供了宝贵的见解，我们将绘制一些时间序列。
时间序列可视化
时间序列数据的主要特征是观测值是在固定的时间间隔内进行的。时间序列值（y轴）与时间本身（x轴）的图表非常重要，并提供了许多结构性的见解。时间序列图不仅仅是带有时间作为x轴的散点图。时间是单调递增的，因此在时间序列图中比散点图中的x轴有更多的意义和重要性。例如，线条可以连接时间序列图中的点，这将指示时间序列的路径，而在散点图中这样的连接是没有意义的，因为它们会散布在各个地方。路径通常表示趋势，并表明序列将向哪个方向移动。时间序列的变化很容易在图表中描绘出来。我们现在将可视化不同的时间序列。
plot.ts函数是这里可视化方案的核心。首先使用windows函数调用一个适当大小的外部图形设备。在 Ubuntu/Linux 中可以使用X11函数。接下来，我们在AirPassengers数据集上运行plot.ts函数：
> windows(height=100,width=150)
> plot.ts(AirPassengers)
以下图表显示了多年间每月乘客数量的增加，平均情况如下：
图 1：每月航空公司乘客计数
我们可以看到，在 12 个月的周期之后，似乎有一个模式在重复，这表明了月份间的季节性趋势。如果能得到一个图表，我们选择第一年，查看月份的图表，然后转到下一年，并展示下一年每月的数据并可视化，以此类推，直到显示完整的数据集，那就太好了。一个名为plotts的函数实现了这种描述的图表，其结构在图 3中给出：
图 2：时间序列频率图函数
plotts函数现在被应用于AirPassengers数据集。该函数位于配套章节代码包的Utilities.R文件中，并在C11.R文件的开始处使用source命令调用。数据包含 12 年的数据，因此结果的时间序列图将包含 12 条曲线。图表的legend需要比通常的区域更大，因此我们将它绘制在图表的右侧。所需的操作通过par、mar和legend函数完成，如下所示：
> par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
> plotts(AirPassengers)
> legend("topright", inset=c(-0.2,0), "-",
+ legend=c(start(AirPassengers)[1]:end(AirPassengers)[1]),
+ col=start(AirPassengers)[1]:end(AirPassengers)[1],lty=2)
我们现在可以清楚地看到以下图中的季节性影响：
图 3：AirPassengers数据集的季节性图
每月乘客计数在二月和十一月达到低点。从二月到七月，每月乘客计数稳步增加，八月保持在相似的水平，然后急剧下降至十一月。十二月和一月可以看到轻微的增加。因此，季节性图提供了更多的见解，它们应该与 plot.ts 函数互补使用。Hyndman 的预测包包含一个名为 seasonalplot 的函数，它实现了与这里定义的 plotts 函数相同的结果。
澳大利亚居民数据集 austres 将在下一部分进行介绍。我们将使用 plotts 函数和图例来增强显示效果：
>plot.ts(austres)
>windows(height=100,width=150)
>par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
>plotts(austres)
>legend("topright", inset=c(-0.2,0), "-",
+ legend=c(start(austres)[1]:end(austres)[1]),
+ col=start(austres)[1]:end(austres)[1],lty=2)
下面的图是澳大利亚居民数量的季度时间序列：
图 4：澳大利亚居民数量的季度时间序列
图 4 和图 5 的季节性图有什么不同？当然，我们寻找的是除了平凡的月度和季度周期性之外的不同之处。在图 5 中，我们可以看到，尽管月度居民数量在季度中有所增加，但这几乎不是一个季节性因素；它似乎更多的是趋势因素而不是季节性因素。因此，与 AirPassengers 数据集相比，季节性贡献似乎较小。
接下来将可视化二氧化碳浓度的时序图。我们在 co2 数据集上使用 plot.ts 函数：
>plot.ts(co2)
运行 plot.ts 函数的结果是下一个输出：
图 5：莫纳罗亚大气二氧化碳浓度可视化
在二氧化碳浓度的时序图中，季节性影响很容易看到。季节性图可能提供更多的见解，我们将在下一步使用 plotts 函数：
>windows(height=100,width=150)
>par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
>plotts(co2)
>legend("topright",inset=c(-0.2,0),
+ "-",
+ legend=c(c(start(co2)[1]:(start(co2)[1]+3)),". . . ",
+ c((end(co2)[1]-3):end(co2)[1])),
+ col=c(c(start(co2)[1]:(start(co2)[1]+3)),NULL,
+ c((end(co2)[1]-3):end(co2)[1])),lty=2)
plotts 和 legend 函数的使用方式与之前相同。程序的结果显示在 图 7 中，我们可以清楚地看到时间序列显示中的季节性影响：
图 6：莫纳罗亚大气二氧化碳浓度的季节性图
练习：使用 forecast 包中的 seasonalplot 函数复制季节性图。这里有什么不同吗？
季节性是时间序列的一个重要方面。及早识别它对于选择适当的时间序列分析模型非常重要。我们将可视化另外三个时间序列数据集，UKDriverDeaths、gas 和 uspop：
>windows(height=100,width=300)
>par(mfrow=c(1,3))
>plot.ts(UKDriverDeaths,main="UK Driver Deaths")
>plot.ts(gas,main="Australian monthly gas production")
>plot.ts(uspop,main="Populations Recorded by the US Census")
图 7中的三个显示与迄今为止所见的不同。看起来之前拟合良好的模型在这里可能不会表现相似。我们看到UKDriverDeaths和gas数据集有很多可变性。对于UKDriverDeaths，看起来在 1983 年之后，死亡人数有所下降。对于gas数据集，我们可以看到在 1970 年之前有一个规律的季节性影响，在那之后，gas的生产力急剧上升。这可能是一些技术突破或其他现象性变化的迹象。可变性也增加了，并且几乎在整个时间范围内都不稳定。uspop显示出指数增长。
练习：视觉检查UKDriverDeaths和gas数据集是否存在季节性影响：
图 7：三个时间序列图：UKDriverDeaths、gas 和 uspop
核心概念和指标
时间序列数据的可视化在所有示例中都传达了类似的故事，从图 1到图 7。观察到的趋势和季节性表明，时间序列的未来值依赖于当前值，因此我们不能假设观察值之间是独立的。但这意味着什么呢？为了重申这一点，考虑更简单的uspop（美国人口）数据集，图 7的第三右面板显示。在这里，我们没有季节性影响。现在，考虑 1900 年的普查年份。下一次普查的人口肯定不会少于 1890 年，也不会远低于同年的同一数量。对于大多数时间序列也是如此；例如，如果我们正在记录一天的最高温度。在这里，如果今天的最高温度是 42°C，那么第二天最高温度将高度受这个数字的影响，并且几乎可以肯定地说，第二天的最高温度不会是 55°C 或 25°C。虽然很容易看出观察值是相互依赖的，但正式的规范本身也是一个挑战。让我们正式介绍时间序列。
我们将用表示在时间
观察到的时序。对于观察到时间 T 的时序，另一种表示方法是
。时序可以被视为在时间 t=1, 2, 3, …观察到的随机过程 Y。与时序过程
相关的是误差过程
。误差过程通常假设为具有零均值和一些常数方差的白噪声过程。误差过程通常被称为创新过程。请注意，时序
可能依赖于
过程中的过去值，以及误差
的值和误差过程的过去值
。值
也被称为
的第一滞后值，
是
的第二滞后值，依此类推。现在，如果观测值相互依赖，那么它们之间关系的指定就是最大的挑战。当然，我们无法在这里详细说明。然而，如果我们认为一阶滞后项是相关的，那么这里必须存在某种关系，我们可以获得一个散点图，其中
的观测值位于 y 轴上，一阶滞后项位于 x 轴上。
AirPassengers数据集的一阶滞后散点图如下获得：
>plot(AirPassengers[1:143],AirPassengers[2:144],
+ xlab="Previous Observation",
+ ylab="Current Observation")
使用索引将时序对象的类别更改为数值对象：
图 8：AirPassengers数据集的滞后图
在先前的图形显示中，我们可以清楚地看到滞后观测值之间存在（几乎）线性关系，因此模型可能的形式为。散点图能否帮助确定依赖的阶数？几乎不能！我们将获得
WWWusage数据集的下一个图，并查看滞后图：
>windows(height=200,width=200)
>par(mfrow=c(2,2))
>plot.ts(WWWusage)
>plot(WWWusage[1:99],WWWusage[2:100],
+ xlab="Previous Observation",
+ ylab="Current Observation",main="Lag-1 Plot"
+ )
>plot(WWWusage[1:98],WWWusage[3:100],
+ xlab="Previous Observation",
+ ylab="Current Observation",main="Lag-2 Plot"
+ )
>plot(WWWusage[1:97],WWWusage[4:100],
+ xlab="Previous Observation",
+ ylab="Current Observation",main="Lag-3 Plot"
+ )
以下是WWWUsage的滞后图：
图 9：WWWUsage的滞后图
第一阶滞后图可能会给人留下观测值相关的印象。然而，更高阶的滞后图几乎没有任何意义，回到第一阶滞后图只会造成更多的困惑。因此，我们需要一种更正式的方法来获得滞后阶数上的洞察。
对于理解时间序列数据中依赖性的本质，有两个有用的度量：自相关函数（ACF）和偏自相关函数（PACF）。正如其名称所示，ACF 是时间序列与其滞后值之间的相关性。PACF 的偏命名法解释了从滞后值中去除中间变量的影响。简单来说，滞后 3 的 PACF 将只包括第一个 ；第三个滞后变量
和
变量不允许影响 PACF。滞后-k ACF 定义为随机变量Y t与第 k 滞后变量Y t-k之间的相关性：
其中 是时间序列的方差。时间序列的偏自相关函数 PACF 是Yt与其第 k 滞后Yt-k之间的偏相关。无法深入探讨 PACF 概念的数学细节；读者可以参考 Box 等人（2015）以获取更多信息。基于n个观察值的样本 ACF 公式如下：
对于 PACF 的显式公式，我们建议读者参考网络上可用的文档，网址为mondi.web.elte.hu/spssdoku/algoritmusok/acf_pacf.pdf。
尽管公式看起来令人畏惧，但我们可以通过在两个数据集austres和uspop上简单地使用acf和pacf函数来轻松地解决这个问题：
>jpeg("ACF_PACF_Plots.jpeg")
>par(mfrow=c(2,2))
>acf(austres,main="ACF of Austres Data")
>pacf(austres,main="PACF of Austres Data")
>acf(uspop,main="ACF of US Population")
>pacf(uspop,main="PACF of US Population")
>dev.off()
RStudioGD
2
我们将保持 ACF 和 PACF 图的解释更加简单。在 ACF 和 PACF 图中，重要的指导原则是水平蓝色线条。任何超出两条线的滞后 ACF 和 PACF 图都是显著的，而那些在限制范围内的则是无关紧要的：
图 10：austres 和 uspop 数据集的 ACF 和 PACF 图
从图 10中我们可以看到，对于austres时间序列，我们需要扩展 ACF 图以包括更多的滞后。这是因为所有绘制的滞后都超出了水平蓝色线条。对于uspop时间序列，第一个时间序列的前四个滞后是显著的，其余的都在水平蓝色线条内。PACF 图可以以类似的方式进行解释。
ACF 和 PACF 图在 ARMA 模型识别中起着重要作用。即使对于 AR 的情况，这些图揭示了滞后信息，也可以用来指定时间序列作为神经网络输入向量的前几个值。
在许多实际问题中，我们还有额外的变量，我们可能将这些称为协变量时间序列或外生变量。让我们用表示协变量时间序列，其中
可能是一个标量或向量时间序列。我们采用惯例，只有
的当前和过去值会影响
，而
的未来值不会以任何方式影响
。也就是说，只有协变量的滞后值会有影响，而不是它们的先行值。在
Car Sales 数据集中，销售是我们感兴趣的时间序列，我们认为广告方面影响了销售；销售不可能解释广告金额！ccf 函数用于以下方式获取交叉相关系数：
>CarSales <- read.csv("../Data/Car_Sales.csv")
>summary(CarSales)
Sales Advertising
Min. :12.0 Min. : 1.0
1st Qu.:20.3 1st Qu.:15.8
Median :24.2 Median :23.0
Mean :24.3 Mean :28.5
3rd Qu.:28.6 3rd Qu.:41.0
Max. :36.5 Max. :65.0
>jpeg("CCF_Car_Sales_Advertising.jpeg")
>ccf(x=CarSales$Advertising,y=CarSales$Sales,
+ main="Cross Correlation Between Sales and Advertising")
>dev.off()
RStudioGD
2
下图展示了销售和广告之间的交叉相关性：
图 11：广告支出和汽车销售的交叉相关系数
我们应该关注正滞后值还是负滞后值？请注意，ccf 图并不是对称的，因此我们不能因为忽略了滞后值的符号而免责。在 R 终端运行 ?ccf，我们得到 ccf(x, y) 返回的滞后 k 值来自帮助文件，该文件估计了 x[t+k] 和 y[t] 之间的相关性。因此，正滞后是先行指标，而负滞后对我们来说更有兴趣。在这个例子中，只有的前置滞后（即
和
）是显著的。
我们在本节结束时简要讨论了准确度测量。与早期学习问题一样，我们将有一系列模型可供我们选择。这是下一节讨论的主题，我们需要相应地定义某些评估指标。让表示时间序列，由于使用某种模型，拟合值将是
。我们可以通过各种方法访问模型的准确性；一些准确性测量包括以下内容：
目前，我们不会关注模型。相反，我们将使用它作为主要工具来提取拟合值并帮助获取定义的指标。使用 subset 函数，我们将定义训练数据集，并使用 forecast 包中的 auto.arima 函数拟合模型。然后，使用 accuracy 和 forecast 函数，我们将获得不同的准确性：
>co2_sub <- subset(co2,start=1,end=443)
>co2_arima <- auto.arima(co2_sub)
>accuracy(forecast(co2_arima,h=25),x=co2[444:468])
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.0185 0.283 0.225 0.00541 0.0672 0.211 0.0119
Test set -0.0332 0.349 0.270 -0.00912 0.0742 0.252 NA
forecast函数是一个非常重要的函数。给定一个拟合的时序，它将提供所需未来期间的预测，并且准确度函数计算了七个不同标准所需的准确度。平均误差标准通常是无用的，对于几乎无偏的模型，其数值将接近 0。RMSE、MAE、MPE 和 MAPE 指标通常很有用，它们的值越低，模型拟合得越好。此外，训练集误差和测试集误差必须可以比较。如果它们彼此差异很大，那么该模型对预测来说就没有用了。在下一节中，我们将回顾一类有用的时序模型。
重要的时序模型
我们已经遇到了不同回归模型的一系列模型。时间序列数据带来了额外的复杂性，因此我们有更多的模型可供选择（或者更确切地说，是从中集成）。这里提供了一个重要模型的快速回顾。这里讨论的大多数模型都处理单变量时间序列，我们需要更多专业化的模型和方法来包含
。我们将从最简单可能的时间序列模型开始，然后逐步过渡到神经网络实现。
简单预测
假设我们拥有的数据，并且我们需要预测下一个h个时间点
。简单预测模型不需要任何建模练习或计算，它只是简单地返回当前值作为未来的预测，因此
。就这么简单。即使是这个简单的任务，我们也会使用 forecast 包中的简单函数，并要求它提供下一个
25个观察值的预测，其中h=25：
>co2_naive <- naive(co2_sub,h=25,level=c(90,95))
>summary(co2_naive)
Forecast method: Naive method
Model Information:
Call: naive(y = co2_sub, h = 25, level = c(90, 95))
Residual sd: 1.1998
Error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1 1.2 1.07 0.029 0.319 0.852 0.705
Forecasts:
Point Forecast Lo 90 Hi 90 Lo 95 Hi 95
Dec 1995 360 358 362 357 362
Jan 1996 360 357 362 356 363
Feb 1996 360 356 363 356 364
Oct 1997 360 350 369 348 371
Nov 1997 360 350 369 348 371
Dec 1997 360 350 370 348 371
如预期的那样，预测值在整个期间保持不变。它们可以很容易地可视化，准确度也可以很容易地按以下方式计算：
>plot(co2_naive) # Output suppressed
>accuracy(forecast(co2_naive,h=25),x=co2[444:468])
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.10 1.20 1.07 0.029 0.319 1.00 0.705
Test set 3.54 4.09 3.55 0.972 0.974 3.32 NA
自然会出现的疑问是简单预测是否有效。对此的回答可以通过其他方式提供。复杂和精密的模型总是声称具有优势。模型确实可能具有优势，但参考和基准应该是清晰的。简单预测提供了这个重要的基准。请注意，对于训练期，简单预测的准确度值是不同的，并且重要的是所提出的模型至少要比简单预测的指标更好。这就是简单预测的主要目的。
季节性、趋势和局部加权回归拟合
季节、趋势和局部加权回归是三个技术术语的组合，形成了 stl 模型。之前，我们在时间序列的视觉显示中看到，其中一些描绘了季节性效应，一些显示了趋势，一些显示了季节性和趋势的组合，还有一些是简单的非规律性时间序列。这些显示因此表明了底层现象的具体性质的存在或不存在。在本节中，我们将探讨如何利用时间序列的季节性和趋势部分。在 loess 中的 stl 模型的第三个组成部分完全没有解释。局部加权回归是一种非参数回归技术，代表局部多项式回归，它将加权最小二乘准则推广到 p 次多项式。局部加权回归方法还包括一个称为核的重要组件。核是一种平滑方法，但我们将不会对此进行过多细节的讨论。
Cleveland 等人（1990 年）提出了基于局部加权回归的季节性-趋势分解，该过程的全部细节可以从以下来源获得：www.scb.se/contentassets/ca21efb41fee47d293bbee5bf7be7fb3/stl-a-seasonal-trend-decomposition-procedure-based-on-loess.pdf。这篇论文易于阅读，直观且富有洞察力，读者应该真正地跟随它。stl 模型是一种分解季节性时间序列的滤波方法，分为三个部分：趋势、季节性和剩余部分。设为时间序列，我们用
、
和
分别表示趋势、季节性和剩余部分；那么我们将有如下：
请参阅 Cleveland 等人发表的论文以获取完整细节。
使用stats包中的stl函数，我们将AirPassengers数据分解如下：
>AP_stl <- stl(AirPassengers,s.window=frequency(AirPassengers))
>summary(AP_stl)
Call:
stl(x = AirPassengers, s.window = frequency(AirPassengers))
Time.series components:
seasonal trend remainder
Min. :-73.3 Min. :123 Min. :-36.2
1st Qu.:-25.1 1st Qu.:183 1st Qu.: -6.4
Median : -5.5 Median :260 Median : 0.3
Mean : 0.1 Mean :280 Mean : -0.2
3rd Qu.: 20.4 3rd Qu.:375 3rd Qu.: 5.9
Max. : 94.8 Max. :497 Max. : 48.6
IQR:
STL.seasonal STL.trend STL.remainder data
46 192 12 180
% 25.2 106.4 6.8 100.0
Weights: all == 1
Other components: List of 5
$ win : Named num [1:3] 12 21 13
$ deg : Named int [1:3] 0 1 1
$ jump : Named num [1:3] 2 3 2
$ inner: int 2
$ outer: int 0
>jpeg("STL_Decompose_AirPassengers.jpeg")
>plot(AP_stl)
>dev.off()
windows
2
>accuracy(forecast(AP_stl))
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.00498 11.2 8.29 -0.129 3.29 0.259 0.000898
执行前面的代码后，可以得到以下图形：
图 12：AirPassengers 的 STL 分解
在stl函数中通过s.window选项指定季节性是很重要的。从季节性图中，我们可以看到每个成分随时间增加。然而，我们可以清楚地看到乘客数量随时间变化的各个组成部分。尽管季节性部分在幅度上增加，但整个期间的模式保持不变。趋势显示线性增加，这表明序列前进的方向。很明显，季节性在这个背景下起着重要作用。
练习：之前已经提到，季节性似乎不是austres数据集分析中的一个有用的因素。使用stl分解并检查这一观察结果是否成立。
接下来将考虑一个更参数化的模型，形式为指数平滑模型。
指数平滑状态空间模型
基本指数平滑模型可以明确定义。用表示平滑因子，并初始化
。基本指数平滑模型定义为以下：
指数模型的详细信息可以在labs.omniti.com/people/jesus/papers/holtwinters.pdf找到。一个更通用的模型形式是指数平滑状态空间模型。在这里，模型在三个方面被指定，就像 stl 模型一样：误差成分、趋势成分，第三个是季节成分。在 forecast 包的ets函数中，成分可以具有加性效应，用A表示，它可以具有乘性效应，用 M 表示，或者可能要求自动选择（Z），并且这种指定适用于每个成分。效应可以指定为既不是加性的也不是乘性的，用字母 N 表示。因此，ets函数中的模型是用第一个字母表示误差成分，第二个字母表示趋势成分，第三个字母表示季节成分来指定的。因此，表示model="ZZZ"意味着三个成分都是自动选择的。model="AMZ"意味着误差成分是加性的，趋势是乘性的，季节成分是自动选择的，等等。Hyndman 等人（2008 年）提供了指数平滑方法详细情况的全面概述。接下来，我们将使用forecast包中的ets函数来拟合指数平滑模型：
>uspop_sub <- subset(uspop,start=1,end=15)
>USpop_ets <- ets(uspop_sub)
>summary(USpop_ets)
ETS(A,A,N)
Call:
ets(y = uspop_sub)
Smoothing parameters:
alpha = 0.8922
beta = 0.8922
Initial states:
l = 2.3837
b = 1.7232
sigma: 1.68
AIC AICc BIC
66.2 72.8 69.7
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 1.11 1.68 1.4 3.26 4.6 0.0318 -0.28
ets函数为误差和趋势成分拟合了加性误差，而选择不对季节因子添加任何内容。这是有意义的，因为uspop数据集没有季节成分。使用这个拟合模型，我们将预测 1940-70 年间的美国人口，使用accuracy函数计算训练和测试数据集的准确性，并与naive预测进行比较：
>forecast(USpop_ets,h=4)
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
1940 139 137 141 136 142
1950 156 151 160 149 162
1960 172 165 180 161 183
1970 189 178 200 173 205
>plot(forecast(USpop_ets,h=4))
>accuracy(forecast(USpop_ets,h=4),x=uspop[16:19])
ME RMSE MAE MPE MAPE MASE ACF1
Training set 1.11 1.68 1.40 3.259 4.60 0.165 -0.28
Test set 2.33 9.02 8.26 0.578 4.86 0.973 NA
>accuracy(forecast(naive(uspop_sub),h=4),x=uspop[16:19])
ME RMSE MAE MPE MAPE MASE ACF1
Training set 8.49 9.97 8.49 21.7 21.7 1.00 0.778
Test set 43.58 51.35 43.58 24.2 24.2 5.13 NA
下图展示了美国人口数据的指数平滑：
图 13：美国人口数据的指数平滑
准确性比较是与简单预测进行的，我们发现ets预测有显著的改进。因此，ets预测是有用的，我们可以使用它们进行未来预测。
接下来，我们将继续介绍流行的 Box-Jenkins/ARIMA 模型。
自回归积分移动平均（ARIMA）模型
Box 和 Jenkins 使用 ARIMA 模型对时间序列的分析和预测方法产生了变化。ARIMA 模型是更一般线性过程模型的一个特例，对于具有创新过程 的时间序列
，其表达式如下：
在这里， 是线性过程的系数。请注意，对创新过程的滞后值没有限制，我们确实意味着在这个线性过程模型中有无限项。在流行的自回归 AR(p)模型中，p 是 AR 模型的阶数。这可以用以下方式表示：
AR 模型可以证明是线性过程模型的一个特例。当时间序列用创新过程表示时，另一个有用的模型是阶数为 q 的移动平均 MA(q) 模型：
时间序列可能不仅依赖于过去值，还依赖于过去的误差，这种结构依赖性在自回归移动平均 ARMA(p,q) 模型中得以捕捉，其阶数为 (p,q)：
p 和 q 的阶数可以通过 表 2 中 ACF 和 PACF 相关的经验规则来确定：
表 2：ARMA 模型的 ACF 和 PACF
ARMA 模型与平稳时间序列数据配合良好，这里的平稳大致意味着序列的变异性在整个过程中保持恒定。这是一个限制性假设，对于许多时间序列现象，这个假设并不成立。在许多实际场景中，通过差分序列 可以获得平稳性，也就是说，我们可以考虑差分
。差分
是一阶差分，有时可能需要更高阶的差分。在大多数实际场景中，差分到 4 阶已被证明可以带来平稳性。差分的阶数通常用字母 d 表示，将 ARMA 模型应用于差分被称为自回归积分移动平均模型，或 ARIMA 模型。简写为
ARIMA(p,d,q)。
在本章中，我们多次遇到了季节成分，它通过季节 ARIMA 模型在 ARIMA 中得到解决。有动机的读者可以通过 Tattar 等人（2017）的第十章了解更多细节。我们在这里简单地补充说，我们还有大写字母（P, D, Q）的进一步符号表示季节参数，以及频率项。我们现在能够理解上一节末尾拟合的模型。co2_arima的准确性已经评估，我们现在将查看摘要：
>summary(co2_arima)
Series: co2_sub
ARIMA(2,1,2)(1,1,2)[12]
Coefficients:
ar1 ar2 ma1 ma2 sar1 sma1 sma2
0.033 0.250 -0.369 -0.246 -0.828 0.014 -0.750
s.e. 0.341 0.122 0.342 0.197 0.230 0.210 0.173
sigma² estimated as 0.0837: log likelihood=-73.4
AIC=163 AICc=163 BIC=195
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.0185 0.283 0.225 0.00541 0.0672 0.179 0.0119
最佳拟合的 ARIMA 模型阶数是(2,1,2)(1,1,2)[12]，这意味着季节频率为12（这是我们已知的事情），季节阶数（P,D,Q）为(1,1,2)，ARIMA 阶数（p,d,q）为(2,1,2)。正是这种差分阶数实现了平稳性。预测结果已获得并进行了可视化：
>jpeg("CO2_Forecasts.jpeg")
>plot(forecast(co2_arima,h=25))
下图显示了输出：
图 14：二氧化碳浓度预测
因此，我们已经为二氧化碳浓度水平拟合了 ARIMA 模型。
接下来，我们将探讨非线性神经网络时间序列模型。
自回归神经网络
神经网络之前已被用于分类以及回归问题。由于时间序列是依赖性观测值，我们需要调整神经网络的架构以包含这种依赖性。调整是为了允许时间序列的滞后值作为输入层的成员。其余的架构遵循与常规神经网络相同的结构。forecast中的nnetar函数代表神经网络自回归，p=选项允许时间序列的滞后值，我们将这些值应用于gas问题：
>gas_sub <- subset(gas,start=1,end=450)
>gas_nnetar <- nnetar(gas_sub,p=25,P=12,size=10,repeats=10)
>plot(forecast(gas_nnetar,h=26))
>accuracy(forecast(gas_nnetar,h=26),x=gas[451:476])
ME RMSE MAE MPE MAPE MASE ACF1
Training set 2 318 237 -0.127 1.78 0.148 -0.0879
Test set 5033 6590 5234 10.566 10.94 3.276 NA
图 15：使用自回归神经网络的气体预测
我们现在已经看到了自回归神经网络的实际应用。
一团糟
我们在本章开头简要介绍了七个数据集，并在表 1中提到了 21 个数据集。数据可视化提供了适度的洞察力，而准确度指标有助于分析模型的有用性。到目前为止，本节已经介绍了一系列模型，现在我们将把一切都搞乱。定义了一个get_Accuracy函数，该函数将拟合六个不同的时间序列模型。代表线性模型的LM尚未解释；TBATS 模型也没有解释。线性模型非常简单，因为时间索引被用作协变量。一般来说，如果一个时间序列有 T 个观测值，协变量向量简单地由值 1，2，3，…，T 组成。我们预计线性模型的表现会不佳。TBATS 模型在此不再详细解释，因此建议阅读一些额外的资料以获取更多相关信息。get_Accuracy模型将每个模型拟合到 21 个数据集，命名模型，并列出其在整个数据集上的性能。以下程序得到所需的结果：
>get_Accuracy<- function(ts){
+ tsname <- deparse(substitute(ts))
+ Acc_Mat <- data.frame(TSName = rep(tsname,6),Models=c(
+ "ETS","STL","LM","ARIMA","NNETAR","TBATS"),
+ ME=numeric(6),RMSE=numeric(6),MAE=numeric(6),
+ MPE=numeric(6), MAPE=numeric(6),MASE=numeric(6))
+ for(i in 1:nrow(Acc_Mat)){
+ Acc_Mat[1,3:8] <- accuracy(ets(ts)$fitted,ts)[1:6]
+ if(frequency(ts)>1) Acc_Mat[2,3:8] <- accuracy(ts-stl(ts,
+ frequency(ts))$time.series[,3],ts)[1:6] else
+ Acc_Mat[2,3:8] <- NA
+ Acc_Mat[3,3:8] <- accuracy(fitted(lm(ts~I(1:length(ts)))),ts)[1:6]
+ Acc_Mat[4,3:8] <- accuracy(auto.arima(ts)$fitted,ts)[1:6]
+ Acc_Mat[5,3:8] <- accuracy(fitted(nnetar(ts)),ts)[1:6]
+ Acc_Mat[6,3:8] <- accuracy(fitted(tbats(ts)),ts)[1:6]
+ }
+ Acc_Mat
+ }
> TSDF <- data.frame(TSName=character(0),Models=character(0),
+ Accuracy=numeric(0))
> TSDF <- rbind(TSDF,get_Accuracy(AirPassengers))
> TSDF <- rbind(TSDF,get_Accuracy(BJsales))
> TSDF <- rbind(TSDF,get_Accuracy(JohnsonJohnson))
> TSDF <- rbind(TSDF,get_Accuracy(LakeHuron))
> TSDF <- rbind(TSDF,get_Accuracy(Nile))
> TSDF <- rbind(TSDF,get_Accuracy(UKgas))
> TSDF <- rbind(TSDF,get_Accuracy(UKDriverDeaths))
> TSDF <- rbind(TSDF,get_Accuracy(USAccDeaths))
> TSDF <- rbind(TSDF,get_Accuracy(WWWusage))
> TSDF <- rbind(TSDF,get_Accuracy(airmiles))
> TSDF <- rbind(TSDF,get_Accuracy(austres))
> TSDF <- rbind(TSDF,get_Accuracy(co2))
> TSDF <- rbind(TSDF,get_Accuracy(discoveries))
> TSDF <- rbind(TSDF,get_Accuracy(lynx))
> TSDF <- rbind(TSDF,get_Accuracy(nhtemp))
> TSDF <- rbind(TSDF,get_Accuracy(nottem))
> TSDF <- rbind(TSDF,get_Accuracy(presidents))
In addition: Warning message:
In ets(ts) :
Missing values encountered. Using longest contiguous portion of time series
> TSDF <- rbind(TSDF,get_Accuracy(treering))
> TSDF <- rbind(TSDF,get_Accuracy(gas))
> TSDF <- rbind(TSDF,get_Accuracy(uspop))
> TSDF <- rbind(TSDF,get_Accuracy(sunspots))
> write.csv(TSDF,"../Output/TS_All_Dataset_Accuracies.csv",row.names=F)
前述代码块输出的表格如下：
| TSName | Model | ME | RMSE | MAE | MPE | MAPE | MASE |
|---|---|---|---|---|---|---|---|
AirPassengers | ETS | 1.5807 | 10.6683 | 7.7281 | 0.4426 | 2.8502 | 0.0164 |
AirPassengers | STL | -0.1613 | 11.9379 | 8.5595 | -0.0662 | 3.4242 | 0.5515 |
AirPassengers | LM | 0.0000 | 45.7362 | 34.4055 | -1.2910 | 12.3190 | 0.7282 |
AirPassengers | ARIMA | 1.3423 | 10.8462 | 7.8675 | 0.4207 | 2.8005 | -0.0012 |
AirPassengers | NNETAR | -0.0118 | 14.3765 | 11.4899 | -0.2964 | 4.2425 | 0.5567 |
AirPassengers | TBATS | 0.4655 | 10.6614 | 7.7206 | 0.2468 | 2.8519 | 0.0215 |
BJsales | ETS | 0.1466 | 1.3272 | 1.0418 | 0.0657 | 0.4587 | -0.0110 |
BJsales | STL | NA | NA | NA | NA | NA | NA |
BJsales | LM | 0.0000 | 9.1504 | 7.1133 | -0.1563 | 3.1686 | 0.9872 |
BJsales | ARIMA | 0.1458 | 1.3281 | 1.0447 | 0.0651 | 0.4601 | -0.0262 |
BJsales | NNETAR | -0.0001 | 1.4111 | 1.0849 | -0.0040 | 0.4798 | 0.2888 |
BJsales | TBATS | 0.1622 | 1.3566 | 1.0666 | 0.0732 | 0.4712 | -0.0113 |
JohnsonJohnson | ETS | 0.0495 | 0.4274 | 0.2850 | 1.0917 | 7.0339 | -0.2948 |
JohnsonJohnson | STL | -0.0088 | 0.1653 | 0.1080 | -0.5953 | 2.8056 | -0.4155 |
JohnsonJohnson | LM | 0.0000 | 1.6508 | 1.3287 | 22.6663 | 66.3896 | 0.6207 |
JohnsonJohnson | ARIMA | 0.0677 | 0.4074 | 0.2676 | 2.0526 | 6.5007 | 0.0101 |
JohnsonJohnson | NNETAR | 0.0003 | 0.3501 | 0.2293 | -0.6856 | 5.8778 | -0.0347 |
JohnsonJohnson | TBATS | 0.0099 | 0.4996 | 0.3115 | 0.9550 | 7.5277 | -0.1084 |
sunspots | ETS | -0.0153 | 15.9356 | 11.2451 | #NAME? | Inf | 0.0615 |
sunspots | STL | 0.0219 | 12.2612 | 8.7973 | NA | Inf | 0.18 |
sunspots | LM | 0 | 42.9054 | 34.1212 | #NAME? | Inf | 0.9196 |
sunspots | ARIMA | -0.0267 | 15.6006 | 11.0258 | NA | Inf | -0.0106 |
太阳黑子 | NNETAR | -0.0672 | 10.3105 | 7.6878 | NA | Inf | 0.0108 |
太阳黑子 | TBATS | -0.0514 | 15.5788 | 11.0119 | NA | Inf | -0.0013 |
表 3：21 个数据集上六个模型的准确度
整体信息与分类问题介绍章节中得到的结果相同。由于我们并不总是能够持续进行结果检查，因此将多个模型的结果结合起来，以传达一个更准确的整体故事是很有吸引力的。我们从这个简单的想法开始，即对指数时间序列模型进行袋装。
袋装和时间序列
在本节中，我们只将说明袋装技术用于 ETS 模型。袋装的主要目的是稳定预测或预报。在这里，我们将基于 Box-Cox 和基于 Loess 的分解进行袋装。使用 500 个这样的自助样本，我们将获得 ETS 的袋装模型：
>uspop_bagg_ets <- baggedETS(uspop_sub,bootstrapped_series =
+ bld.mbb.bootstrap(uspop_sub, 500))
>forecast(uspop_bagg_ets,h=4);subset(uspop,start=16,end=19)
Point Forecast Lo 100 Hi 100
1940 141 136 145
1950 158 150 165
1960 175 164 184
1970 193 178 204
Time Series:
Start = 1940
End = 1970
Frequency = 0.1
[1] 132 151 179 203
>plot(forecast(uspop_bagg_ets,h=4))
使用袋装方法是否有优势？我们可以通过置信区间快速检查这一点：
>forecast(uspop_bagg_ets,h=4)
Point Forecast Lo 100 Hi 100
1940 141 136 145
1950 158 150 165
1960 175 164 184
1970 193 178 204
>forecast(USpop_ets,h=4,level=99.99)
Point Forecast Lo 99.99 Hi 99.99
1940 139 133 146
1950 156 142 169
1960 172 150 194
1970 189 157 221
袋装 ETS 的置信区间明显更短，因此反映了方差的减少，这是袋装的主要目的：
图 16：美国人口袋装预测
准确度比较也容易在这里进行：
>accuracy(forecast(USpop_ets,h=4),x=uspop[16:19])
ME RMSE MAE MPE MAPE MASE ACF1
Training set 1.11 1.68 1.40 3.259 4.60 0.165 -0.28
Test set 2.33 9.02 8.26 0.578 4.86 0.973 NA
>accuracy(forecast(uspop_bagg_ets,h=4),x=subset(uspop,start=16,end=19))
ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
Training set 1.137 1.44 1.24 2.226 4.48 0.0283 0.563 NA
Test set -0.359 7.87 7.48 -0.995 4.63 0.1700 0.296 0.299
这里清楚地显示了集成同质基学习者的优势。接下来，我们将转向异质基学习者和它们的集成。
集成时间序列模型
forecastHybrid R 包提供了一个平台来集成异构时间序列模型。实现这一任务的主要功能是 hybridModel 函数。核心函数提供了称为 models 的选项。它接受一个最多六个字符的字符串，这些字符是模型的代表：a 代表 auto.arima 模型，e 代表 ets，f 代表 thetam，n 表示 nnetar，s 代表 stlm，最后，t 代表 tbats。因此，如果我们给 models 一个 ae 字符串，它将结合 ARIMA 和 ets 模型的结果。这可以在 co2 数据集上对不同时间序列模型的组合进行说明：
>accuracy(forecast(co2_arima,h=25),x=co2[444:468])
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.0185 0.283 0.225 0.00541 0.0672 0.211 0.0119
Test set -0.0332 0.349 0.270 -0.00912 0.0742 0.252 NA
>AP_Ensemble_02 <- hybridModel(co2_sub,models="ae")
Fitting the auto.arima model
Fitting the ets model
>accuracy(AP_Ensemble_02,h=25,x=co2[444:468])
ME RMSE MAE MPE MAPE ACF1 Theil's U
Test set 0.0258 0.271 0.219 0.00755 0.0653 0.00289 0.226
>AP_Ensemble_03 <- hybridModel(co2_sub,models="aen")
Fitting the auto.arima model
Fitting the ets model
Fitting the nnetar model
>accuracy(AP_Ensemble_03,h=25,x=co2[444:468])
ME RMSE MAE MPE MAPE ACF1 Theil's U
Test set 0.017 0.304 0.245 0.0049 0.073 0.282 0.25
>AP_Ensemble_04 <- hybridModel(co2_sub,models="aens")
Fitting the auto.arima model
Fitting the ets model
Fitting the nnetar model
Fitting the stlm model
>accuracy(AP_Ensemble_04,h=25,x=co2[444:468])
ME RMSE MAE MPE MAPE ACF1 Theil's U
Test set 0.0165 0.275 0.221 0.00478 0.066 0.209 0.226
>AP_Ensemble_05 <- hybridModel(co2_sub,models="aenst")
Fitting the auto.arima model
Fitting the ets model
Fitting the nnetar model
Fitting the stlm model
Fitting the tbats model
>accuracy(AP_Ensemble_05,h=25,x=co2[444:468])
ME RMSE MAE MPE MAPE ACF1 Theil's U
Test set 0.0123 0.267 0.216 0.00348 0.0645 0.153 0.22
虽然这里的集成讨论非常简短，但时间序列文献最近才开始适应集成技术。
练习：weights 和 errorMethod 选项对于将不同的时间序列模型组合在一起至关重要。探索章节中介绍和讨论的不同数据集的这些选项。
摘要
时间序列数据带来了新的挑战和复杂性。本章从介绍重要且流行的数据集开始。我们探讨了不同的时间序列及其复杂性。时间序列的可视化提供了深刻的洞察，时间序列图以及季节性图互补地用于清晰的想法和特定实施。时间序列的准确度指标各不相同，我们探讨了其中的一些。自相关函数（ACF）和偏自相关函数（PACF）的概念在模型识别中至关重要，季节性成分对于时间序列建模也同样重要。我们还看到，不同的模型可以表达不同的数据集，变化的程度类似于常规回归问题。时间序列的袋装（仅 ets）减少了预测的方差。在结论部分讨论了结合异构基学习器。下一章是结论章节。我们将总结前十一章的主要收获，并概述一些不足和进一步的研究范围。
