| Title: | Estimating Remaining Useful Life with Linear Mixed Effects Models |
|---|---|
| Description: | Provides tools for estimating the Remaining Useful Life (RUL) of degrading systems using linear mixed-effects models and creating a health index. It supports both univariate and multivariate degradation signals. For multivariate inputs, the signals are merged into a univariate health index prior to modeling. Linear and exponential degradation trajectories are supported (the latter using a log transformation). Remaining Useful Life (RUL) distributions are estimated using Bayesian updating for new units, enabling on-site predictive maintenance. Based on the methodology of Liu and Huang (2016) <doi:10.1109/TASE.2014.2349733>. |
| Authors: | Pedro Abraham Montoya Calzada [aut, cre, cph] (ORCID: <https://orcid.org/0009-0002-3497-210X>), Rogelio Salinas Gutiérrez [aut, cph] (ORCID: <https://orcid.org/0000-0002-1669-4460>), Silvia Rodríguez-Narciso [aut, cph] (ORCID: <https://orcid.org/0000-0001-5429-5914>), Netzahualcóyotl Castañeda-Leyva [aut, cph] (ORCID: <https://orcid.org/0000-0001-9414-3923>) |
| Maintainer: | Pedro Abraham Montoya Calzada <[email protected]> |
| License: | GPL-3 |
| Version: | 1.0.1 |
| Built: | 2026-06-05 07:01:15 UTC |
| Source: | https://github.com/abraham1011/degradr |
Given a fitted "healthindex" object, this function constructs the univariate health index for new multivariate sensor data by applying the stored projection (weights and offsets).
compute_healthindex(model, data)compute_healthindex(model, data)
model |
An object of class |
data |
A data frame with new sensor readings over time. Must include the columns |
This function applies the projection learned in the first stage of fit_healthindex to new data:
A data frame (tibble) with the columns:
unit |
Unit identifier. |
t |
Time index. |
x |
Constructed health index at each |
fit_healthindex for learning the health index and mixed-effects model,
predict_rul for RUL prediction based on the fitted model.
library(degradr) library(dplyr) # Load example data data(train_FD001) data(test_FD001) data <- train_FD001 %>% select(unit,t,T24,T50,P30, Nf,Ps30,phi, NRf, BPR,htBleed, W31, W32) %>% mutate(across(c(P30,phi,W31,W32), ~ . * -1)) test <- test_FD001 %>% select(unit,t,T24,T50,P30, Nf,Ps30,phi, NRf, BPR,htBleed, W31, W32) %>% mutate(across(c(P30,phi,W31,W32), ~ . * -1)) # Fit a health index model (exponential trajectory of degree 2) model <- fit_healthindex(data = data, type = "exponential", degree = 2, r = 0.8) # Construct the health index on new data using stored weights/offsets hi_new <- compute_healthindex(model = model, data = test) head(hi_new)library(degradr) library(dplyr) # Load example data data(train_FD001) data(test_FD001) data <- train_FD001 %>% select(unit,t,T24,T50,P30, Nf,Ps30,phi, NRf, BPR,htBleed, W31, W32) %>% mutate(across(c(P30,phi,W31,W32), ~ . * -1)) test <- test_FD001 %>% select(unit,t,T24,T50,P30, Nf,Ps30,phi, NRf, BPR,htBleed, W31, W32) %>% mutate(across(c(P30,phi,W31,W32), ~ . * -1)) # Fit a health index model (exponential trajectory of degree 2) model <- fit_healthindex(data = data, type = "exponential", degree = 2, r = 0.8) # Construct the health index on new data using stored weights/offsets hi_new <- compute_healthindex(model = model, data = test) head(hi_new)
Test data from a real-world degradation process involving the clogging of gas filters. The dataset includes right-censored lifetimes due to a preventive maintenance policy. Failure is defined when pressure differential exceeds 600 Pa.
data(filter_test)data(filter_test)
A data frame with 50 degradation trajectories. Variables include:
Time
Differential pressure across the filter in Pascal (Pa).
Unique identifier for each filter test unit. Integer from 1 to 50.
Remaining Useful Life in hours. Represents the time left until failure (threshold = 600 Pa).
Adapted from:
Hagmeyer, S., Mauthe, F., & Zeiler, P. (2021). Creation of Publicly Available Data Sets for Prognostics and Diagnostics Addressing Data Scenarios Relevant to Industrial Applications. International Journal of Prognostics and Health Management, 12(2). doi:10.36001/ijphm.2021.v12i2.3087
See full description in the dataset repository on Kaggle by "Prognostics @ HSE".
Training data from a real-world degradation process involving the clogging of gas filters. The dataset includes right-censored lifetimes due to a preventive maintenance policy. Failure is defined when pressure differential exceeds 600 Pa.
data(filter_train)data(filter_train)
A data frame with 49 degradation trajectories. Variables include:
Time
Differential pressure across the filter in Pascal (Pa).
Unique identifier for each filter test unit. Integer from 1 to 50.
Adapted from:
Hagmeyer, S., Mauthe, F., & Zeiler, P. (2021). Creation of Publicly Available Data Sets for Prognostics and Diagnostics Addressing Data Scenarios Relevant to Industrial Applications. International Journal of Prognostics and Health Management, 12(2). doi:10.36001/ijphm.2021.v12i2.3087
See full description in the dataset repository on Kaggle by "Prognostics @ HSE".
Fits a health index–based degradation model by projecting multivariate sensor signals into a univariate health index and modeling its evolution using a linear mixed-effects model.
fit_healthindex(data, type = "exponential", method = "lm", degree = 2, phi = NULL, r = 0.5)fit_healthindex(data, type = "exponential", method = "lm", degree = 2, phi = NULL, r = 0.5)
data |
A data frame containing sensor readings over time. Must include the columns |
type |
Model type. Either |
method |
Estimation method. Either |
degree |
Degree of the polynomial model. Default is |
phi |
Initial degradation level for non-defective units.
Used in the exponential model as a fixed offset to ensure that the logarithmic transformation is valid and interpretable.
If |
r |
parameter that controls the relative importance of the threshold variance and the weighted residual sum of squares in the index-fitted degradation model. |
This function implements a two-stage modeling strategy. In the first stage, a univariate health index is constructed as a weighted linear combination of the input signals, using correlation-based shrinkage. In the second stage, the resulting health index is modeled over time with a linear mixed-effects model (on the log scale for exponential models).
The exponential model uses a log transformation of x - phi, where phi ensures positivity and interpretability. The phi parameter can be supplied by the user or estimated automatically.
The resulting object stores both the projection (health index definition) and the fitted model used for RUL prediction.
Returns an object of class "healthindex", which contains:
index |
A list with components: weights |
model |
A fitted mixed-effects model of the health index over time. |
Liu, K. and Huang, S. (2016). Integration of Data Fusion Methodology and Degradation Modeling Process to Improve Prognostics. IEEE Transactions on Automation Science and Engineering, 13(1), 344–354.doi:10.1109/TASE.2014.2349733
library(degradr) library(dplyr) # Load example data data(train_FD001) data(test_FD001) data <- train_FD001 %>% select(unit,t,T24,T50,P30, Nf,Ps30,phi, NRf, BPR,htBleed, W31, W32) %>% mutate(across(c(P30,phi,W31,W32), ~ . * -1)) test <- test_FD001 %>% select(unit,t,T24,T50,P30, Nf,Ps30,phi, NRf, BPR,htBleed, W31, W32) %>% mutate(across(c(P30,phi,W31,W32), ~ . * -1)) # Fit a health index model (exponential trajectory of degree 2) model <- fit_healthindex(data = data, type = "exponential", degree = 2, r = 0.8) rul <- predict_rul(data = test, model = model) head(rul)library(degradr) library(dplyr) # Load example data data(train_FD001) data(test_FD001) data <- train_FD001 %>% select(unit,t,T24,T50,P30, Nf,Ps30,phi, NRf, BPR,htBleed, W31, W32) %>% mutate(across(c(P30,phi,W31,W32), ~ . * -1)) test <- test_FD001 %>% select(unit,t,T24,T50,P30, Nf,Ps30,phi, NRf, BPR,htBleed, W31, W32) %>% mutate(across(c(P30,phi,W31,W32), ~ . * -1)) # Fit a health index model (exponential trajectory of degree 2) model <- fit_healthindex(data = data, type = "exponential", degree = 2, r = 0.8) rul <- predict_rul(data = test, model = model) head(rul)
Fits a linear or exponential mixed-effects model of degree \(p\) for the degradation process.
fit_model(data, type = "exponential", method = "lm", degree = 2, phi = NULL)fit_model(data, type = "exponential", method = "lm", degree = 2, phi = NULL)
data |
A data frame with three columns: |
type |
Model type. Either |
method |
Estimation method. Either |
degree |
Degree of the polynomial model. Default is |
phi |
Initial degradation level for non-defective units.
Used in the exponential model as a fixed offset to ensure that the logarithmic transformation is valid and interpretable.
If |
This function fits a linear or exponential polynomial mixed-effects model of degree p to degradation data collected over time from multiple units. The model captures both fixed effects (population-level degradation trends) and random effects (unit-specific deviations).
The exponential model applies a logarithmic transformation with an offset parameter phi. The offset phi can be provided or automatically estimated from the data.
At least two distinct units are required to estimate random effects. The degree parameter controls the polynomial order for the time terms in both fixed and random effects.
Returns a list with the estimated model and prior distributions.
Liu, K. and Huang, S. (2016). Integration of Data Fusion Methodology and Degradation Modeling Process to Improve Prognostics. IEEE Transactions on Automation Science and Engineering, 13(1), 344–354.doi:10.1109/TASE.2014.2349733
library(degradr) # Load example data sets data(filter_train) data(filter_test) # Show the original column names colnames(filter_train) # Rename the columns to match the expected format: t, x, unit colnames(filter_train) <- c("t", "x", "unit") colnames(filter_test) <- c("t", "x", "unit", "RUL") # Plot the training set plot_degradr(data = filter_train, D = 600) # Fit an exponential mixed-effects model of degree 1 model <- fit_model(data = filter_train, type = "exponential", degree = 1) # Predict the remaining useful life (RUL) for the test units, # assuming a fixed failure threshold D = 600 predict_rul(data = filter_test, model = model, D = 600)library(degradr) # Load example data sets data(filter_train) data(filter_test) # Show the original column names colnames(filter_train) # Rename the columns to match the expected format: t, x, unit colnames(filter_train) <- c("t", "x", "unit") colnames(filter_test) <- c("t", "x", "unit", "RUL") # Plot the training set plot_degradr(data = filter_train, D = 600) # Fit an exponential mixed-effects model of degree 1 model <- fit_model(data = filter_train, type = "exponential", degree = 1) # Predict the remaining useful life (RUL) for the test units, # assuming a fixed failure threshold D = 600 predict_rul(data = filter_test, model = model, D = 600)
Generates a line plot of degradation signals over time for each unit, optionally overlaying a failure threshold line. This function is useful for visualizing degradation paths across multiple components or systems.
plot_degradr(data, D = NULL)plot_degradr(data, D = NULL)
data |
A data frame with three columns: |
D |
Optional numeric value indicating the failure threshold. |
The function is designed to work with degradation datasets where each row represents an observation of a unit at a particular time. The plot shows how the degradation variable x evolves over time t for each unit. This is especially useful for visual inspection before model fitting or threshold analysis.
Returns a ggplot object that can be further customized or directly printed.
library(degradr) # Load example data sets data(filter_train) data(filter_test) # Show the original column names colnames(filter_train) # Rename the columns to match the expected format: t, x, unit colnames(filter_train) <- c("t", "x", "unit") plot_degradr(data = filter_train, D = 600)library(degradr) # Load example data sets data(filter_train) data(filter_test) # Show the original column names colnames(filter_train) # Rename the columns to match the expected format: t, x, unit colnames(filter_train) <- c("t", "x", "unit") plot_degradr(data = filter_train, D = 600)
Estimates the Remaining Useful Life (RUL) for one or more partially observed degradation signals based on a previously fitted linear or exponential mixed-effects model.
predict_rul(data, model, D = NULL, upper = NULL)predict_rul(data, model, D = NULL, upper = NULL)
data |
A data frame with columns |
model |
An object of class |
D |
(Optional) Critical degradation threshold. If provided, it will be used to compute the RUL via a fixed-threshold model. If |
upper |
Optional upper bound for the search interval when solving for the quantiles of the RUL distribution. If |
This function applies Bayesian updating to compute the posterior distribution of the degradation model parameters for each unit, conditional on its observed signal. Then, it computes the Remaining Life Distribution (RLD) and returns the estimated Remaining Useful Life.
It supports both linear and exponential degradation models, matching the formulation used in fit_model. The posterior updating follows the methodology of Liu and Huang (2016).
A data frame with one row per unit and the following columns:
Unit identifier.
Estimated RUL.
Liu, K. and Huang, S. (2016). Integration of Data Fusion Methodology and Degradation Modeling Process to Improve Prognostics. IEEE Transactions on Automation Science and Engineering, 13(1), 344–354. doi:10.1109/TASE.2014.2349733
library(degradr) # Load example data sets data(filter_train) data(filter_test) # Show the original column names colnames(filter_train) # Rename the columns to match the expected format: t, x, unit colnames(filter_train) <- c("t", "x", "unit") colnames(filter_test) <- c("t", "x", "unit", "RUL") # Plot the training set plot_degradr(data = filter_train, D = 600) # Fit an exponential mixed-effects model of degree 1 model <- fit_model(data = filter_train, type = "exponential", degree = 1) # Predict the remaining useful life (RUL) for the test units, # assuming a fixed failure threshold D = 600 predict_rul(data = filter_test, model = model, D = 600)library(degradr) # Load example data sets data(filter_train) data(filter_test) # Show the original column names colnames(filter_train) # Rename the columns to match the expected format: t, x, unit colnames(filter_train) <- c("t", "x", "unit") colnames(filter_test) <- c("t", "x", "unit", "RUL") # Plot the training set plot_degradr(data = filter_train, D = 600) # Fit an exponential mixed-effects model of degree 1 model <- fit_model(data = filter_train, type = "exponential", degree = 1) # Predict the remaining useful life (RUL) for the test units, # assuming a fixed failure threshold D = 600 predict_rul(data = filter_test, model = model, D = 600)
Evaluates the cumulative probability that the Remaining Useful Life (RUL) of a unit is less than or equal to a specified time t. The computation is based on a fitted degradation model and the observed degradation signal for the unit.
prul(t, data, model, D = NULL)prul(t, data, model, D = NULL)
t |
Time at which to evaluate the RUL cumulative distribution function. |
data |
A data frame with columns |
model |
An object of class |
D |
Optional critical degradation threshold. If provided, a fixed-threshold model is used; otherwise a random-threshold model is assumed. For exponential models the threshold is automatically transformed to the log scale. |
For a fixed threshold model (D supplied), the function computes the Remaining Life Distribution (RLD) using the specified failure threshold. If D is NULL, the distribution is computed under a random-threshold formulation based on the training data.
Numeric value between 0 and 1 giving .
library(degradr) library(dplyr) # Load example data data(train_FD001) data(test_FD001) data <- train_FD001 %>% select(unit,t,T24,T50,P30, Nf,Ps30,phi, NRf, BPR,htBleed, W31, W32) %>% mutate(across(c(P30,phi,W31,W32), ~ . * -1)) test <- test_FD001 %>% select(unit,t,T24,T50,P30, Nf,Ps30,phi, NRf, BPR,htBleed, W31, W32) %>% mutate(across(c(P30,phi,W31,W32), ~ . * -1)) # Fit a health index model (exponential trajectory of degree 2) model <- fit_healthindex(data = data, type = "exponential", degree = 2, r = 0.8) #Probability that the run length will be less than or equal to 86 cycles head(prul(t = 86, data = test, model = model))library(degradr) library(dplyr) # Load example data data(train_FD001) data(test_FD001) data <- train_FD001 %>% select(unit,t,T24,T50,P30, Nf,Ps30,phi, NRf, BPR,htBleed, W31, W32) %>% mutate(across(c(P30,phi,W31,W32), ~ . * -1)) test <- test_FD001 %>% select(unit,t,T24,T50,P30, Nf,Ps30,phi, NRf, BPR,htBleed, W31, W32) %>% mutate(across(c(P30,phi,W31,W32), ~ . * -1)) # Fit a health index model (exponential trajectory of degree 2) model <- fit_healthindex(data = data, type = "exponential", degree = 2, r = 0.8) #Probability that the run length will be less than or equal to 86 cycles head(prul(t = 86, data = test, model = model))
Returns quantiles of the Remaining Life Distribution for one or more units based on their observed degradation signals and a fitted model.
qrul(prob = 0.05, data, model, D = NULL, upper = NULL)qrul(prob = 0.05, data, model, D = NULL, upper = NULL)
prob |
Probability at which to evaluate the quantile of the RUL distribution. |
data |
A data frame with columns |
model |
An object of class |
D |
Optional critical degradation threshold. If |
upper |
Optional upper bound for the numerical search when computing quantiles. If |
For each unit in data the function computes the posterior distribution of the model parameters and evaluates the specified quantile of the Remaining Life Distribution. Units that result in computational errors return NA.
A data frame with one row per unit containing:
unit |
Unit identifier. |
RUL |
The requested quantile of the RUL distribution. |
library(degradr) library(dplyr) # Load example data data(train_FD001) data(test_FD001) data <- train_FD001 %>% select(unit,t,T24,T50,P30, Nf,Ps30,phi, NRf, BPR,htBleed, W31, W32) %>% mutate(across(c(P30,phi,W31,W32), ~ . * -1)) test <- test_FD001 %>% select(unit,t,T24,T50,P30, Nf,Ps30,phi, NRf, BPR,htBleed, W31, W32) %>% mutate(across(c(P30,phi,W31,W32), ~ . * -1)) # Fit a health index model (exponential trajectory of degree 2) model <- fit_healthindex(data = data, type = "exponential", degree = 2, r = 0.8) head(qrul(prob = 0.05, data = test, model = model))library(degradr) library(dplyr) # Load example data data(train_FD001) data(test_FD001) data <- train_FD001 %>% select(unit,t,T24,T50,P30, Nf,Ps30,phi, NRf, BPR,htBleed, W31, W32) %>% mutate(across(c(P30,phi,W31,W32), ~ . * -1)) test <- test_FD001 %>% select(unit,t,T24,T50,P30, Nf,Ps30,phi, NRf, BPR,htBleed, W31, W32) %>% mutate(across(c(P30,phi,W31,W32), ~ . * -1)) # Fit a health index model (exponential trajectory of degree 2) model <- fit_healthindex(data = data, type = "exponential", degree = 2, r = 0.8) head(qrul(prob = 0.05, data = test, model = model))
Ground truth values of the Remaining Useful Life (RUL) for each engine unit in the FD001 subset of the C-MAPSS dataset. These values correspond to the last observed cycle in the test set and are used for evaluation purposes in prognostics models.
data("RUL_FD001")data("RUL_FD001")
A data frame with 100 observations on the following variable:
RULNumeric vector indicating the true Remaining Useful Life (in cycles) for each unit in the test set.
This dataset is part of the C-MAPSS (Commercial Modular Aero-Propulsion System Simulation) benchmark and is used as ground truth for performance evaluation of predictive maintenance models, particularly for models estimating Remaining Useful Life (RUL) under the FD001 operating condition scenario.
Saxena, A., Goebel, K., Simon, D., & Eklund, N. (2008). Damage propagation modeling for aircraft engine run-to-failure simulation. In 2008 International Conference on Prognostics and Health Management (pp. 1–9). IEEE. doi:10.1109/PHM.2008.4711414
data(RUL_FD001)data(RUL_FD001)
Truncated time series data from turbofan engine degradation simulations, generated using the C-MAPSS (Commercial Modular Aero-Propulsion System Simulation) model. This test dataset contains operational sensor data up to a point before failure, intended for validating prognostic algorithms that estimate Remaining Useful Life (RUL).
A data frame with multiple observations (rows) on the following 24 variables (columns):
unitEngine unit number (identifier)
tTime in cycles
T2Total temperature at fan inlet (°R)
T24Total temperature at LPC outlet (°R)
T30Total temperature at HPC outlet (°R)
T50Total temperature at LPT outlet (°R)
P2Pressure at fan inlet (psia)
P15Total pressure in bypass-duct (psia)
P30Total pressure at HPC outlet (psia)
NfPhysical fan speed (rpm)
NcPhysical core speed (rpm)
eprEngine pressure ratio (P50/P2)
Ps30Static pressure at HPC outlet (psia)
phiRatio of fuel flow to Ps30 (pps/psi)
NRfCorrected fan speed (rpm)
NRcCorrected core speed (rpm)
BPRBypass Ratio
farBBurner fuel-air ratio
htBleedBleed Enthalpy
Nf_dmdDemanded fan speed (rpm)
PCNfR_dmdDemanded corrected fan speed (rpm)
W31HPT coolant bleed (lbm/s)
W32LPT coolant bleed (lbm/s)
Key characteristics of this test dataset:
Simulates progressive degradation in the High Pressure Compressor (HPC) module
Time series are truncated prior to failure (true RUL values not included)
Includes realistic measurement noise and unit-to-unit variability
Saxena, A., Goebel, K., Simon, D., & Eklund, N. (2008). Damage propagation modeling for aircraft engine run-to-failure simulation. In 2008 International Conference on Prognostics and Health Management (pp. 1–9). IEEE. doi:10.1109/PHM.2008.4711414
data(test_FD001)data(test_FD001)
Run-to-failure simulation data for aircraft turbofan engines generated using C-MAPSS (Commercial Modular Aero-Propulsion System Simulation). This dataset represents engine degradation in the High Pressure Compressor (HPC) module under varying operational conditions.
data("train_FD001")data("train_FD001")
A data frame with multiple observations (rows) on the following 24 variables (columns):
unitEngine unit number (identifier)
tTime in cycles
T2Total temperature at fan inlet (°R)
T24Total temperature at LPC outlet (°R)
T30Total temperature at HPC outlet (°R)
T50Total temperature at LPT outlet (°R)
P2Pressure at fan inlet (psia)
P15Total pressure in bypass-duct (psia)
P30Total pressure at HPC outlet (psia)
NfPhysical fan speed (rpm)
NcPhysical core speed (rpm)
eprEngine pressure ratio (P50/P2)
Ps30Static pressure at HPC outlet (psia)
phiRatio of fuel flow to Ps30 (pps/psi)
NRfCorrected fan speed (rpm)
NRcCorrected core speed (rpm)
BPRBypass Ratio
farBBurner fuel-air ratio
htBleedBleed Enthalpy
Nf_dmdDemanded fan speed (rpm)
PCNfR_dmdDemanded corrected fan speed (rpm)
W31HPT coolant bleed (lbm/s)
W32LPT coolant bleed (lbm/s)
The data was generated for the Prognostics and Health Management (PHM) 2008 data challenge. Each engine unit starts with some initial wear and progresses to failure as efficiency and flow parameters degrade exponentially. The failure criterion is when the health index (calculated from stall margins and EGT) reaches zero. The dataset includes sensor measurements taken at cruise conditions.
Saxena, A., Goebel, K., Simon, D., & Eklund, N. (2008). Damage propagation modeling for aircraft engine run-to-failure simulation. In 2008 International Conference on Prognostics and Health Management (pp. 1–9). IEEE. doi:10.1109/PHM.2008.4711414
data(train_FD001)data(train_FD001)