---

2019 年 5 月

第 34 卷，第 5 期

[机器学习]

对预测性维护使用生存分析

通过Zvi Topol |2019 年 5

几年前，我引入了生存分析的基础知识并介绍了如何实现非参数化的算法调用 Kaplan Meier C# (msdn.com/magazine/dn630650)。现在，我将看看生存分析，特别是在两个更高级的方法可方便地使用两个常用机器学习平台，Spark 机器学习库 (MLLib) 上并h2o.ai，同时支持哪些 Azure HDInsight。我将与正在进行的示例使用预测性维护用例。

预测性维护的工业物联网

工业物联网 (IIoT) 的主要思想是连接计算机、 设备、 传感器和工业设备和组织内的应用程序，并不断地从所有中收集数据，如系统错误和机器遥测这些做的目的是分析和操作对此数据以优化运营效率。

预测性维护的目标是准确地预测在一台计算机时或其任意组件将失败。如果可以这样做，您可以在预测此类故障发生之前执行维护。这是不执行任何维护，直到发生故障，在这种情况下机器或组件之前将不可失败固定的如果确实为可修复比效率更高。此类非计划停机时间很可能开销非常大。

预测性维护也很多有效比也可能是考虑的因为不必要的维护可能应用的频繁执行预防性维护。

示例和我将使用的数据是在该示例的改编的版bit.ly/2J4WnbN。该示例包含 100 台生产计算机，与在计算机之间不相互依赖关系。每台计算机是四个可能的模型之一。

计算机的数据包括故障、 维护操作和传感器的遥测，以及有关模型和使用年限 （年） 的计算机信息的历史记录。在.csv 文件可从前面所述的资源会提供此数据。我还将提供转换的数据文件 (comp1_df.csv) 的"生存分析就绪"，并将介绍如何更高版本上执行的转换。

在原始的示例中的每台计算机具有四个不同的组件，但我将重点介绍仅在一个组件上。组件可以在故障时之前, 主动维护或后未能修复它保留。

生存分析

在有关生存分析上一篇文章，我介绍了重要的基本概念，我将使用并扩展这篇文章中。建议您阅读本文来熟悉这些概念，包括的生存和危险函数、 审查和非参数化 Kaplan Meier (KM) 估计器。

在本文中，我将介绍如何扩展 KM 估算器，以包含 covariates 或可能会产生影响上生存，或者，在这种情况下，在机组件失败的变量 （也称为功能） 的概念。在示例中，我将使用与 covariates 机器模型、 机年龄和机器遥测，并使用生存回归模型来估计此类 covariates 对计算机发生故障的影响。 

这一概念估计的 covariates 效果的目标变量，在这种情况下故障、 危险率或生存概率的时间、 不是唯一的生存分析和一般情况下是回归模型的基础。

当生成的统计模型，请参阅 covariates 的三个主要的数据类型： 分类、 序号和连续。分类数据类型是这些类型可分为几个离散的类别。在这里，机器模型是分类数据类型，有四种不同的机模型。序号数据类型是分类数据类型具有一些有意义的顺序。例如，从 1 到 10，10 是最有趣，另一个电影的评级最少。最后，持续的数据类型是那些表示连续数值。这些是机遥测读数在这里，这是连续数字 （在这种情况下，每小时），在某些时候采样。

确定要使用的数据类型和方法后, 应编码为 covariates 的各种数据类型。通常情况下，对于回归模型，连续变量自然会编码为连续 covariates，而分类数据类型将需要某种形式的编码。我将在本文中使用，此类编码的常用选项是在其中，N 类别分类的数据类型 N-1 covariates 创建，并由其他人设置其特定的变数，并且所有值为零表示类别 i。通过将所有 covariates 都设置为零表示第 n 个类别。这通常是很适合回归模型与显式定义的基线，其中所有 covariates 可以都是等于零。这也是 R 编程语言使用分类变量或因素进行编码的格式。

Categoricals 此编码具有意义的部分或全部 covariates，若要设置为零的简单解释。但是，对于连续数据类型设置为零的某些变数始终可能有意义。例如，如果变数表示计算机高度或宽度，设置为零的变数将无意义的因为在现实中没有此类计算机。

一种方法解决此问题是使用 mean 居中的连续 covariates，其中对于给定的变数，训练数据集其平均值来代表其值减去。然后，当将该转换后的变数设置为零，则相当于设置为其平均值的原始变数。此技术称为"mean 居中"，我将使用在此处机年龄和遥测 covariates。

请务必记住，按照此转换，您应始终使用 mean 居中的 covariates 作为模型的输入。将该回归模型应用到新的测试数据集时，这也是如此。

这些 covariates 和某种形式的生存目标变量 （我将讨论很快） 后的数据值被编码为 covariates，然后执行生存回归模型，并指定将对生存/时间事件的此类 covariates 影响的模型.

生存格式和功能设计到数据的转换

若要使用我将介绍生存回归模型，你的数据需要具有至少两个字段： 感兴趣 （此处，计算机发生故障） 字段和一个布尔值，该值指示是否审查发生的事件的时间戳。（在这里，审查描述或指定时间之前发生任何失败的情况。在示例中，维护预防性的方式，而不是作为对失败的响应被视为进行审查。

我将讨论生存回归模型具有进行以简化其数学派生的各种假设。一些这些假设可能不能够保持在这里，但仍可用于将应用到此示例中建模的生存。

生存分析宣传资料是非常丰富和许多高级生存回归模型和技术已开发的地址，放宽某些这些假设。你可以阅读更多有关此类模型和在书中，"统计分析的故障时间数据"Kalbfleisch 和 Prentice (Wiley-Interscience，2002年)，在技术bit.ly/2TACdLR。

我将进行完全执行计算机组件的每个维护操作将重置该组件的假设，并因此可以独立处理。这样，它将在两种类型的时间间隔上使用生存回归 (中所示图 1):

• 失败和前面的维护操作 （事件时间） 之间的间隔。
• 后面的维护操作 （审查） 之间的间隔。

图 1 生存表示形式的计算机故障

在每个间隔图 1开头维护操作。第一种类型的 X，表示失败，而第二个类型结尾 O，表示另一个维护操作之前失败 （这是实质上是预防性维护操作），这在本例中意味着审查的观察间隔结束。

因此，原始数据需要转换成此格式与两个必填字段。"Time_to_event"字段表示小时，直至发生故障或下一个维护的时间。"事件"字段设置为一个失败，并为零，以便在发生故障之前执行维护操作。

它是经常需要执行其他转换上 covariates，这通常称为"特征工程。" 此过程的目的是生成 covariates 具有更好地预测能力。例如，可以创建另一个将在失败前 10 个小时内计算方面的压力的平均值的变数。有许多不同的选项为函数和可能的时间范围内创建此类 covariates，并且有几个工具可用来帮助自动执行此过程，如开源 Python 包 tsfresh (tsfresh.readthedocs.io/en/最新)。

现在，我将讨论两个生存回归模型： Cox 成比例的危险模型 （或 Cox p H 模型） h2o.ai 和 Spark MLLib 中可用的韦伯加速故障时间模型中可用。

Cox 成比例的危险回归

前面曾提到危险函数在时间 t 的对象或个人时间 t 处处于活动状态，确定事件速率。用于预测性维护的示例中，可以被描述为下一个小时，为给定的时间 t 和组件 1 故障尚未自其上次维护的发生位置的所有计算机中失败的可能性。更高的危险率意味着遇到故障的风险更高。Cox p H 回归估计对以下模型由指定的危险速率 covariates 的影响：

在这里，h （t） 在时间 t 处是危险函数、 h0(t) 在时间 t 处为基线危险、 Xi 变量是不同 covariates 和相应的测试版是对应于 covariates 的系数 (将详细介绍一些更高版本)。基线危险为危险时所有 covariates 都都等于零。请注意，这与密切相关中其他回归模型，如在线性或逻辑回归的截距。

根据此模型中，covariates 和生存时间之间没有直接关系。此模型称为半参数，因为时间 t 处的危险速率是参数化的乘法组件和基线危险速率，估计的数据并不具有参数化的关闭窗体的函数。

此模型就是成比例的危险模式的原因是因为它允许您进行比较的两个危险函数比率。在这种情况下，给定估计的模型下，是两个不同的数据点之间的比率：

基线危险速率取消和危险生成比率是仅的系数和 covariates 函数并再次不依赖于时间。这是密切相关于逻辑回归，其中估计可能的日志。此外，Cox p H 回归模型不直接指定生存函数，并侧重于提供的危险函数比例的比率上的信息。因此，主要用于了解留存能力，对 covariates 的影响而不是直接估计生存函数。

Cox ph 指定模型、 系数和非参数化基线危险可以估计使用各种技术。一种常用方法是分部的最大可能性估计 （也在 h2o.ai 中使用）。

下面的代码段是 R 脚本，以在运行的 mean 居中 covariates （机器遥测和年龄） 上使用 h2o.ai Cox p H 模型和分类变数机器模型估计值：

library(h2o)
localH2O <- h2o.init()
inputFileName<-'comp1_df.csv'
df<-read.csv(inputFileName, header=TRUE, stringsAsFactors=TRUE)
df.hex <- as.h2o(df, key = "df.hex")
model <- h2o.coxph(x = c("age_mean_centered", "model","volt_mean_centered",
  "rotate_mean_centered","pressure_mean_centered", "vibration_mean_centered"),
  event_column = "event", stop_column ="time_to_event" ,training_frame = df.hex)
summary(model)

在撰写本文时，在 h2o.ai Cox p H 模型不可用来使用从 Python，因此提供了 R 代码。在提供了安装说明bit.ly/2z2QweL，或者，对于使用 Azure HDInsight 的 h2o.ai 处bit.ly/2J7nXp6。

运行代码段将生成的输出中所示图 2。

图 2 输出 Cox p H 回归

Surv(time_to_event, event) ~ model + volt_mean_centered + rotate_mean_centered +
    pressure_mean_centered + vibration_mean_centered + age_mean_centered
  n= 709, number of events= 192
                               coef      exp(coef) se(coef)    z      Pr(>|z|)   
model.model2                  -0.066955  0.935237  0.257424 -0.260    0.795   
model.model3                  -0.021837  0.978400  0.215614 -0.101    0.919   
model.model4                   0.308878  1.361896  0.227469  1.358    0.174   
volt_mean_centered             0.031903  1.032418  0.003990  7.995    1.33e-15 ***
rotate_mean_centered           0.001632  1.001633  0.001362  1.199    0.231   
pressure_mean_centered        -0.008164  0.991869  0.005768 -1.415    0.157   
vibration_mean_centered        0.018220  1.018387  0.013866  1.314    0.189   
age_mean_centered              0.004804  1.004815  0.013293  0.361    0.718   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
                                  exp(coef) exp(-coef) lower .95 upper .95
model.model2                      0.9352     1.0692    0.5647     1.549
model.model3                      0.9784     1.0221    0.6412     1.493
model.model4                      1.3619     0.7343    0.8720     2.127
volt_mean_centered                1.0324     0.9686    1.0244     1.041
rotate_mean_centered              1.0016     0.9984    0.9990     1.004
pressure_mean_centered            0.9919     1.0082    0.9807     1.003
vibration_mean_centered           1.0184     0.9819    0.9911     1.046
age_mean_centered                 1.0048     0.9952    0.9790     1.031
Rsquare= 0.094   (max possible= 0.941 )
Likelihood ratio test= 70.1  on 8 df,   p=4.69e-12
Wald test            = 70.19  on 8 df,   p=4.514e-12

需要注意的第一个重要的一点是估计的 covariates 系数。机模型变数编码为分类数据类型。此类别的基线是 model1，设置三个 covariates 编码其他三个机器模型 （model.model2、 model.model3 和 model.model4） 为零表示。每个变数获取自己的系数。了解如何解释这些系数非常重要。

如果应用的指数函数的系数机模型 covariates (exp(coeff) 输出中)，将看到该 model.model2 的值为 0.9352，而 model.model4 1.3619 的值。这意味着 model2 机具有危险率低于基线机器模型 （模型 1） 的危险速率 6.5 %model.model4 机具有 36.2%相比 model.model1 机具有高很多风险。换而言之，model.model4 机具有最高风险的故障，而 model.model2 机具有失败的最小的风险。因此，当按优先级排列的维护操作计算机的模型应需要考虑的一个重要因素。

所有其他 covariates 是 mean 居中连续 covariates。与它们相关联的系数的解释是，现在由指数 covariates 围绕其表示的给定危险比率。因此，通过增加一个单位 （保持固定的所有其他 covariates） 变数值，危险比率增加 （或减少） 的指数的系数 (类似的方式与分类变量)。因此，例如，通过增加一个单位电压，失败的风险会增加 3.2%。

此处提及的另一个要点涉及模型诊断技术。此类技术提供基础若要了解是否适当 （在此情况下，Cox p H 模型） 被视为模型。在这里，Rsquare 值 （介于零和一，越高越好） 是相对较低 (0.094) 和大部分 z-分数的系数不表示系数是具有统计学意义 （没有足够的证据来支持的它们不同于零）。这两个这些指示器会导致的结论是，没有改进，例如通过特征工程的空间。也有其他特定于应进行 Cox p H 模型的统计测试。有关详细信息，可以前面参考我提到过的生存分析宣传资料。

加速的故障时间模型

Spark MLLib 中的生存回归模型是加速故障时间 (AFT) 模型。此模型中直接指定一个生存函数中的某些理论的数学分布 （韦伯），并具有加速的故障时间属性。

AFT 模型定义，如下所示。假设对象的特征是使用 （线性） covariates 和系数：

此外假定该对象具有参数化生存函数 s(t) 并且，为由 s0(t)，（具有设置为零的所有 covariates) 的基线对象的生存函数。AFT 模型定义了 s(t) 和作为 s0(t) 之间的关系：

此定义中可以看到该模型调用加速故障时间模型的原因。这是因为生存函数包含一个快捷键因素，它是相乘生存时间 t covariates 的线性组合的指数函数。

某些 covariates，如年龄时，此类型的模型非常有用 （我的数据集，在计算机年龄），可能会造成单调加速或减速效果的生存/故障时间。

韦伯分布是指数分布的泛化，并且是参数化生存模型中常用的连续分布。有几个有关如何将其参数化的变体。在这里，我将使用以下两个参数韦伯分布版本为 t > = 0:

（也有三个参数的版本。） 分发的两个参数是由 k 和由 lambda 的小数位数的形状。粗略的类比是钟形分布具有特征平均值和标准偏差的方法。

回想一下，f(t) 由给定分发密度函数 f(t)、 危险函数 h （t） 和生存函数 s(t) 之间的关系 = h(t)s(t)。

以下是韦伯危险和生存函数：

与 Cox p H 模型不同的生存和危险函数完全指定，并且具有参数化表示形式。请参阅图 3并图 4为 k 和 lambda 的不同值的韦伯分布和生存函数的可视化效果。

图 3 韦伯分布形状的不同值的 K 和 Lambda 函数

图 4 韦伯生存函数形状的不同值的 K 和 Lambda

图 5说明了 AFT 模型 covariates 韦伯生存函数形状上的影响。

图 5 为韦伯生存概率函数加速故障时间

完成的 Spark MLLib 中的 AFT 韦伯模型的系数估计使用最大可能性估计算法。你可以了解有关如何在完成详细bit.ly/2XSauom，并查找在实现代码bit.ly/2HtJw0v。

与不同的 Cox p H 模型，其中仅 covariates 的系数报告 （以及一些诊断），估计结果中获得从估计韦伯 AFT 模型报表 covariates，以及特定参数的系数有关韦伯分布 — 截获和小数位数参数。我将介绍如何将转换到 k 和在某种程度的 lambda。

有关韦伯 AFT 实现 Spark MLLib 中的结果匹配韦伯 AFT 实现从流行的 R 库"生存"使用 survreg 函数的结果 (更多详细信息位于bit.ly/2XSxkw8)。

您可以运行下面的 R 脚本 AFT 韦伯模型估计 (代码是运行在本地安装的 Spark MLLi，但也可以在 HDInsight 上使用 Spark bit.ly/2u2U5Qf):

library(survival)
library(SparkR, lib.loc = c(file.path(Sys.getenv("SPARK_HOME"), "R", "lib")))
sparkR.session(master = "local[*]")
inputFileName<-'comp1_df.csv'
df<-read.csv(inputFileName, header=TRUE, stringsAsFactors=TRUE)
aftMachineDF<-suppressWarnings(createDataFrame(df))
aftMachineModel <- spark.survreg(aftMachineDF,Surv(time_to_event, event) ~ model +
  age_mean_centered +  volt_mean_centered + rotate_mean_centered +
  pressure_mean_centered + vibration_mean_centered)
summary(aftMachineModel)

该脚本会生成不包含其他信息仅估计的系数。就可以通过运行 survreg （因为结果匹配） 来获取此类信息：

library(survival)
machineModel<-survreg(Surv(time_to_event, event) ~ model + age_mean_centered +
  volt_mean_centered+rotate_mean_centered + pressure_mean_centered +
  vibration_mean_centered, df, dist='weibull')
summary(machineModel)

在这种情况下，R 脚本生成更复杂的输出中所示图 6。

图 6 输出韦伯 AFT 回归

survreg(formula = Surv(time_to_event, event) ~ model + age_mean_centered +
    volt_mean_centered + rotate_mean_centered + pressure_mean_centered +
    vibration_mean_centered, data = df, dist = "weibull")
                              Value        Std. Error    z            p
(Intercept)                 8.172991       0.119133    68.6040     0.00e+00
modelmodel2                 0.040289       0.154668    0.2605      7.94e-01
modelmodel3                 0.027225       0.129629    0.2100      8.34e-01
modelmodel4                  -0.163865      0.136382   -1.2015     2.30e-01
age_mean_centered           -0.000753      0.007960    -0.0946     9.25e-01
volt_mean_centered           -0.019731      0.002583   -7.6391     2.19e-14
rotate_mean_centered         -0.000767      0.000821   -0.9334     3.51e-01
pressure_mean_centered       0.005173       0.003496     .4795     1.39e-01
vibration_mean_centered      -0.008214      0.008391   -0.9789     3.28e-01
Log(scale)                   -0.508060      0.051963   -9.7773     1.41e-22
Scale= 0.602
Weibull distribution
Loglik(model)= -1710.3   Loglik(intercept only)= -1747.2
                           Chisq= 73.73 on 8 degrees of freedom, p= 8.9e-13
Number of Newton-Raphson Iterations: 8
n= 709

继续之前来描述输出，我应提到韦伯参数化和 survreg Spark MLLib 中是稍微不同于我讨论了参数化。

转换是必需的可以按如下所示。表示参数报告 — 由 m 和小数位数由 s 截距，然后 k = 1/秒，lambda = exp(-m/s) 和每个应乘以系数 (-1/s)。没有 R 包调用可以自动执行此转换在模型上使用 ConvertWeibull 的 SurvRegCensCov 该 survreg 估计：

$vars
                                            Estimate        SE
lambda                                    1.260459e-06   8.642772e-07
gamma                                     1.662064e+00   8.636644e-02
modelmodel2                               -6.696297e-02  2.569595e-01
modelmodel3                               -4.524990e-02  2.155000e-01
modelmodel4                               2.723541e-01   2.268785e-01
age_mean_centered                         1.251958e-03   1.322780e-02
volt_mean_centered                        3.279500e-02   3.947495e-03
rotate_mean_centered                      1.274045e-03   1.365339e-03
pressure_mean_centered                    -8.598142e-03  5.807130e-03
vibration_mean_centered                   .365213e-02    1.391255e-02

在这里，gamma 等于 k 从上一韦伯参数化。(有关 SurvRegCensCov 的详细信息，请参阅bit.ly/2CgcSMg。)

使用 Cox p H 模型不同于给定估计的参数，则现在可以直接获取生存函数 （它是韦伯 AFT 生存函数），并使用它来预测任何 covariates 生存概率。假设在数据集中的第一个点一个新的数据点，您可以运行以下命令：

predict(machineModel, newdata=df[1,], type='quantile')

这将生成的时间 （以小时为单位） 的事件分位数 0.1 和 0.9 （默认值），如下所示：

807.967 5168.231

也就是说，给定的 （此处列出） 的第一个数据点 covariates，故障的概率为 10%或只是之前 807.967 小时遵循一个维护操作，并且故障的概率为 90%或只是之前 5168.231 小时以下维护操作：

model       age           volt_mean_centered                 rotate_mean_centered
 model3       18                  3.322762                          51.8113
pressure_mean_centered     vibration_mean_centered             age_mean_centered
 10.10773                         11.4267                           6.488011

此外可以使用参数"p"若要获取任何零和一; 之间的分位数的生存时间例如，添加参数"p = 0.5"将提供中间故障时间，即，第一个数据点，2509.814 小时后维护操作。

与 Cox p H 模型估计一样，survreg 的输出中的 p 列提供有关的统计意义的估计，不过在这种情况下将数字可以更好的系数 （较低 p-值） 的信息。仍然是特征工程此处为 Cox p H 模型前面描述的空间。

也很重要执行模型诊断结果，因为在 Cox p H 回归中，以确保韦伯 AFT 模型相比，例如，其他参数的模型的数据非常适合这种情况。虽然我不会介绍此过程，可以通过引用我之前提到的"生存分析"一书来了解有关它的详细信息。

总结

我在介绍使用预测性维护的 IIoT 作为有趣的示例示例采用两个生存回归模型中的可用h2o.ai和 Spark MLLib。我介绍了如何为 covariates 并转换为生存格式的时间系列数据通过编码变量建模生存分析框架中的计算机故障预测性维护问题。我还介绍了两个生存模型，以及如何将它们应用于数据之间的差异。最后，我简要提到过的结果和模型诊断解释。请务必注意，我只是简要的本主题中有吸引力且非常丰富，并且建议您浏览的详细信息。通过引用我的文章中提到的资料因此是一个起始点，这样做。

---

Zvi Topol一直在各种行业垂直领域，包括市场营销分析、 媒体和娱乐和工业物联网数据科学家就。他已传递，并会导致多个机器学习和分析项目，包括自然语言和语音接口、 认知搜索、 视频分析、 推荐器系统和市场营销决策支持系统。Topol 目前 MuyVentive llc，高级的分析 R 和 D 公司，并且可以到达zvi.topol@muyventive.com。

衷心感谢以下 Microsoft 技术专家对本文的审阅：James McCaffrey

---

在 MSDN 杂志论坛讨论这篇文章