VAM an R package for Virtual Age Models
Introduction
VAM (Virtual Age Models) is a free R package for imperfect maintenance times data analyzing. R is an open source programming language and software environment for statistical computing and graphics. The capabilities of R are extended through user-created package, such as VAM.
VAM takes as input data sets consisting of successive times of preventive and corrective maintenance performed on repairable systems, possibly left and right-censored.
It implements the main imperfect maintenance models, of the virtual age type.
It performs both frequentist and Bayesian parametric estimation.
It computes and plots reliability indicators such as virtual age, failure intensity, conditional reliability function, ...
It allows to simulate the successive maintenance times and types in order to perform Monte Carlo analyses or bootstrap methods.
The strengths of VAM are:
It is time efficient because, even if it is an R package, its algorithms are developed in C++ using the rcpp package.
It can be adapted to many different situations thanks to the fact that the implemented model is generic and specified by a formula defined by the user. The formula describes the data set, the intrinsic ageing model, the number, the type and the effect models of maintenances, and how the different types of preventive maintenance are planned. For example, the user can choose (with his own formula) to analyze a model corresponding to the simultaneous mixing of Corrective Maintenances (CM) following an Arithmetic Reduction of Age model with memory 3 (ARA3), periodic Preventive Maintenances (PM) following a Quasi Renewal (QR) model also called geometric process, renewing when the age of the system reaches a predetermined level. In the Bayesian context, the formula also specifies the prior distribution of the different parameters.
VAM has been developed by
- L. Doyen (),
- R. Drouilhet (),
both members of the laboratory LJK, of Univ. Grenoble Alpes.
Download
Before being able to use VAM, R has to be install on your computer. Then, download the archive file corresponding to your operating system.
Windows, VAM_0.3.2.zip.
Mac OS, VAM_0.3.1.tgz.
Finally, from R, install the VAM package using the downloaded archive file (for example using the install.packages command of R).
Oldest available versions of the VAM packages are:
- Windows, VAM_0.3.1.zip, VAM_0.3.0.zip, VAM_0.1.1.zip, VAM_0.1.0.zip, VAM_0.0.1.zip,
- Mac OS, VAM_0.3.0.tgz, VAM_0.1.1.tgz, VAM_0.1.0.tgz, VAM_0.0.2.tgz, VAM_0.0.1.tgz.
For experimented users, who are familiarized with rccp object programming, you can find the source code of VAM on GitHub.
Quick Start Tutorial
For detailed explanation, notice that the different commands of the VAM package are documented in classical R help.
Imperfect maintenance models
First, load the VAM package.
require(VAM)
R> require(VAM)
Let us begin with a first data set proposed by C.W. Ahn, K.C. Chae and G.M. Clark (1998: "Estimating parameters of the power-law process with two measures of time" Journal of Quality Technology 30(2), pp. 127–132). This data set has also been studied in a Bayesian context by M. Guilda and G. Pulcini (2006: "Bayesian analysis of repairable systems showing a bounded failure intensity" Reliability Enginering and System Safety 91, pp. 828–838) and F. Corset, L. Doyen and O. Gaudoin(2012: "Bayesian analysis of ARA imperfect repair models" Communications in Statistics - Theory and Methods 41, pp. 3915-3941). It corresponds to the 18th successive failure times of an AMC Ambassodor car owned by the Ohio state government. In VAM, maintenance durations are not considered, consequently those failure times also correspond to the successive Corrective Maintenance (CM) times.
data("AMC_Amb") AMC_Amb
R> data("AMC_Amb") R> AMC_Amb Time Type 1 202 -1 2 265 -1 3 363 -1 4 508 -1 5 571 -1 6 755 -1 7 770 -1 8 818 -1 9 868 -1 10 999 -1 11 1054 -1 12 1068 -1 13 1108 -1 14 1230 -1 15 1268 -1 16 1330 -1 17 1376 -1 18 1447 -1
If you click on the previous R code window, you can visualize the corresponding outputs. You can also directly copy the version without the outputs on your R session and execute those instructions on your own computer. On the outputs, you can notice that AMC_Amb is a two columns data frame. The first column, denoted Time, contains the successive event times and the second column, denoted Type represents the corresponding event types. In this case, all the types are equal to -1 since the system is only submitted to CM actions.
An R formula contains a ~ symbol. The right-hand part of the formula (after the ~) specifies the model to be used. And the left part describes the data to be explained with this model. To define an imperfect CM model, you have to specify the stochastic behavior of the new unmaintained system and the CM effect. For the time to failure of the new unmaintained system, classical distributions of survival analysis can be used. VAM proposed for example Weibull or LogLinear (equivalent to Gompertz) distributions. The basic assumptions on maintenance effect are known as As Bad As Old (ABAO) and As Good As New (AGAN). In the AGAN case, each maintenance is supposed to renew the system, that is to say, after a maintenance the system is equivalent to a new one with a (virtual) age equal to the time elapsed since the last (perfect) maintenance. In the ABAO case, maintenances effects are minimal and leaves the system in the same state as it was just before failure. These behaviors can be illustrated with the following VAM instructions.
AMC_AGAN<-model.vam(Time & Type ~ (AGAN()|Weibull(1,3)) ,data=AMC_Amb) plot(AMC_AGAN ,"v")
AMC_ABAO<-model.vam(Time & Type ~ (ABAO()|Weibull(1,3)) ,data=AMC_Amb) plot(AMC_ABAO ,"v")
You can click on the header of the previous R code session in order to visualize the different basic assumptions models and the corresponding VAM formulas. The previously used model.vam object allows to apply a model to a data set in order to plot different characteristics, such as virtual age (characterized by the "v" argument in the plot function). We can also plot, for example, the failure intensity (argument "i"), which corresponds to conditional failure rate given all the past of the maintenance process. Many other plots are available in VAM.
plot(AMC_ABAO ,"i")
plot(AMC_ABAO ,"i") curve(1*3*x^(3-1) ,add=TRUE ,col="red", lty="dashed" ,lwd=4)
plot(AMC_AGAN ,"i")
plot(AMC_AGAN ,"i") curve(1*3*x^(3-1) ,add=TRUE ,col="red", lty="dashed" ,lwd=4)
AMC_ABAO_LL<-model.vam(Time & Type ~ (ABAO()|LogLinear(0.1,8e-3)) ,data=AMC_Amb) plot(AMC_ABAO_LL ,"i") curve(0.1*exp(8e-3*x) ,add=TRUE ,col="red", lty="dashed" ,lwd=4)
AMC_AGAN_LL<-model.vam(Time & Type ~ (AGAN()|LogLinear(0.1,8e-3)) ,data=AMC_Amb) plot(AMC_AGAN_LL ,"i") curve(0.1*exp(8e-3*x) ,add=TRUE ,col="red", lty="dashed" ,lwd=4)
You can click on the heading of the previous R session and notice that,
the stars on the abscissa axis correspond to the failure (or equivalently CM) times,
up to the first failure, the intensity corresponds to the failure rate of the new unmaintained system (for example a Weibull hazard rate, that is to say a power function with the appropriate parameters corresponding to the one defined in the previous
model.vamformula:(1,3));the failure intensity is equal to the failure rate of the new unmaintained system (for example the failure rate of the Weibull distribution) at the virtual age time.
Obviously, reality falls between the extrem AGAN and ABAO cases: standard maintenance has an effect on the system without necessary renewing it. The VAM package considers that the maintenance effect can only affects the virtual age of the system. The following models for imperfect maintenance effect are implemented in VAM. You can click on each model to have more explanations and to see the corresponding examples of VAM instructions.
The ARA∞ model (ARAInf) has been introduced by Doyen and Gaudoin (2004: "Classes of imperfect repair models based on reduction of failure intensity or virtual age" Reliability Engineering and System Safety 84(1), pp. 45-56). It is equivalent to a Kijima type II model (1989: "Some results for repairable systems with general repair" Journal of Applied Probability 26, pp. 89-102) with deterministic constant effect. The effect of such maintenance is to reduce the virtual age proportionally to its value just before maintenance. The corresponding proportionality factor is a model parameter characterizing the maintenance efficiency. If it is equal to 0.5 the virtual age is reduced by a half after each maintenance.
AMC_ARAInf<-model.vam(Time & Type ~ (ARAInf(0.5)|Weibull(1,3)) ,data=AMC_Amb) plot(AMC_ARAInf ,"v")
AMC_ARAInf_0<-model.vam(Time & Type ~ (ARAInf(0)|Weibull(1,3)) ,data=AMC_Amb) plot(AMC_ARAInf_0 ,"v")
AMC_ARAInf_1<-model.vam(Time & Type ~ (ARAInf(1)|Weibull(1,3)) ,data=AMC_Amb) plot(AMC_ARAInf_1 ,"v")
You can click on the heading of the previous R session and notice that
maintenance is ABAO, if the ARA∞ maintenance efficiency parameter is equal to 0,
maintenance is AGAN, if the ARA∞ maintenance efficiency parameter is equal to 1.
The ARA1 model (ARA1) has been introduced by Doyen and Gaudoin (2004: "Classes of imperfect repair models based on reduction of failure intensity or virtual age" Reliability Engineering and System Safety 84(1), pp. 45-56). It is equivalent to a Kijima type I model (1989: "Some results for repairable systems with general repair" Journal of Applied Probability 26, pp. 89-102) with deterministic constant effect. The effect of such maintenance is to reduce the virtual age proportionally to the supplement of age accumulated since the last maintenance. With a maintenance efficiency parameter equals to 0.5, this is only the supplement increment of age accumulated since the last maintenance that is reduced by a half.
AMC_ARA1<-model.vam(Time & Type ~ (ARA1(0.5)|Weibull(1,3)) ,data=AMC_Amb) plot(AMC_ARA1 ,"v")
AMC_ARA1_0<-model.vam(Time & Type ~ (ARA1(0)|Weibull(1,3)) ,data=AMC_Amb) plot(AMC_ARA1_0 ,"v")
AMC_ARA1_1<-model.vam(Time & Type ~ (ARA1(1)|Weibull(1,3)) ,data=AMC_Amb) plot(AMC_ARA1_1 ,"v")
The ARAm model (ARAm) has been introduced by Doyen and Gaudoin (2004: "Classes of imperfect repair models based on reduction of failure intensity or virtual age" Reliability Engineering and System Safety 84(1), pp. 45-56). The idea of the ARA∞ model is to reduce the whole virtual age whereas in the ARA1 model only the last supplement of age is reduced. The natural idea of the ARAm model is to propose an intermediate effect on the m last times between maintenances.
AMC_ARA3<-model.vam(Time & Type ~ (ARAm(0.5|3)|Weibull(1,3)) ,data=AMC_Amb) plot(AMC_ARA3 ,"v")
AMC_ARA3_0<-model.vam(Time & Type ~ (ARAm(0|3)|Weibull(1,3)) ,data=AMC_Amb) plot(AMC_ARA3_0 ,"v")
AMC_ARA3_1<-model.vam(Time & Type ~ (ARAm(1|3)|Weibull(1,3)) ,data=AMC_Amb) plot(AMC_ARA3_1 ,"v")
The Quasi Renewal model (QR), also called geometric process (Y. Lam. 2007: "The Geometric Process and its Applications" World Scientic) assumes that the successive times between maintenances follow the same distribution as Y1, aY2, a2Y3, ..., where Y1, Y2, Y3, ... are independent and identically distributed random variables. We can prove that this assumption is equivalent to consider that maintenance effect reset the virtual age to 0 but also multiply the virtual age time speed by a certain parameter which characterizes the maintenance efficiency (this parameter is equal to 1/a). Consequently, with a maintenance efficiency parameter equals to 0.8, the slope of the virtual age is equal to:
1, up to the first maintenance time,
0.8, between the first and second maintenance times,
0.82, between the second and third maintenance times, ...
AMC_QR<-model.vam(Time & Type ~ (QR(0.8)|Weibull(1,3)) ,data=AMC_Amb) plot(AMC_QR ,"v")
AMC_QR_bis<-model.vam(Time & Type ~ (QR(1.2)|Weibull(1,3)) ,data=AMC_Amb) plot(AMC_QR_bis ,"v")
The geometric evolution of the times between maintenances in the QR model often seems too much explosive to be realistic from a practical point of view. That is why L. Bordes and S. Mercier (2013: "Extended geometric processes: Semiparametric estimation and application to reliability" Journal of the Iranian Statistical Society 12(1), pp. 1-34) have proposed the extended geometric process ,called GQR in the VAM package. The idea is simply to replace an in the QR model by af(n) where f() stands for a non-decreasing function. In the VAM package, the implemented values for f are indentity (GQR is then equivalent to QR), sqrt, log.
AMC_GQR_sqrt<-model.vam(Time & Type ~ (GQR(1.5|sqrt)|Weibull(1,3)) ,data=AMC_Amb) plot(AMC_GQR_sqrt ,"v")
AMC_GQR_log<-model.vam(Time & Type ~ (GQR(1.5|log)|Weibull(1,3)) ,data=AMC_Amb) plot(AMC_GQR_log ,"v")
ARA models have an immediate effect and reduce the virtual age from a certain quantity. GQR modify the wearing-out time speed but necessarily reset the virtual age. With GQR_ARAInf, GQR_ARA1 or GQR_ARAm, both effects can be mixed. Such maintenance effect models depend on two parameters. The first parameter refers to GQR efficiency and the second one to ARA efficiency.
AMC_GQRARAInf<-model.vam(Time & Type ~ (GQR_ARAInf(1.1,0.5|identity)|Weibull(1,3)) ,data=AMC_Amb) plot(AMC_GQRARAInf ,"v")
AMC_GQRARAm<-model.vam(Time & Type ~ (GQR_ARAm(0.6,0.3|log,3)|Weibull(1,3)) ,data=AMC_Amb) plot(AMC_GQRARAm ,"v")
Simulation
The sim.vam object creates a simulator that can be launched thanks to the simulate function.
simARAInf<-sim.vam(Time & Type ~ (ARAInf(0.6)|Weibull(1,3))) simData<-simulate(simARAInf,5) simData
R> simARAInf<-sim.vam(Time & Type ~ (ARAInf(0.6)|Weibull(1,3))) R> simData<-simulate(simARAInf,5) R> simData Time Type 1 1.190677 -1 2 1.487090 -1 3 2.014754 -1 4 2.901206 -1 5 2.971390 -1
The last simulated data set is memorized in the corresponding sim.vam object and then plots can directly be generated using sim.vam objects.
simulate(simARAInf,5) plot(simARAInf,"i",col="blue") simulate(simARAInf,3) plot(simARAInf,"i",col="red",add=TRUE)
Many different stopping criteria for data simulation are available in VAM. For example, the user can choose to stop the data generation at a specified time t instead of stopping at the nth maintenance time.
simulate(simARAInf,Time>(RightCensorship=3))
R> simulate(simARAInf,Time>(RightCensorship=3)) Time Type 1 1.195878 -1 2 1.721718 -1 3 1.926600 -1 4 2.020666 -1 5 3.000000 0
You can click on the heading of the previous R session and notice that the last maintenance time is equal to 3, i.e. the right censoring time. But the associated type is 0, since it is a censoring time.
Finally, VAM can also generate several independent and identical systems.
simulate(simARAInf,3,nb.system=2)
R> simulate(simARAInf,3,nb.system=2) System Time Type 1 1 0.9322119 -1 2 1 1.9418285 -1 3 1 2.2236236 -1 4 2 0.9823267 -1 5 2 1.2097444 -1 6 2 2.2021332 -1
You can click on the heading of the previous R session and notice that the generated data frame has now 3 columns. The first column System specifies for each event to which system it corresponds.
Parameters estimation
Maximum Likelihood Estimators
On practical situations, the different parameters of the model are generally unknown. But they can be estimated. The mle.vam object defines a Maximum Likelihood Estimation method for a data set, and the corresponding parameter estimation can be numerically computed using for example the coef function.
AMC_mle<-mle.vam(Time & Type ~ (ARAInf(0.6)|Weibull(1,3)),data=AMC_Amb) coef(AMC_mle)
R> AMC_mle<-mle.vam(Time & Type ~ (ARAInf(0.6)|Weibull(1,3)),data=AMC_Amb) R> coef(AMC_mle) [1] 2.120545e-09 3.582878e+00 2.457933e-01
The obtained values correspond to the parameters estimations for the distribution of the new unmaintained system (in this case Weibull), and after to the different maintenance efficiencies in the same order as they appear in the formula. The parameter values in the mle.vam formula are used for the first initialization of the numerical maximization method.
The estimated values of the parameters are memorized in the mle.vam object and then the plug-in estimators of the different plots can directly been obtained with the plot function. For example, we can compare the plug-in estimator of the cumulative failure intensity (or compensator of the failure process), argument "i", with the observed counting process.
plot(AMC_mle,"I")
The function logLik.mle.vam provides the corresponding value of the log-likelihood and its derivatives. That can be used, for example, to compute AIC criteria.
(AIC<-2*3-2*logLik.mle.vam(AMC_mle))
R> (AIC<-2*3-2*logLik.mle.vam(AMC_mle)) [1] 191.3556
Bayesian estimators
The bayesian.vam object defines a Bayesian estimation method. The parameters values, in the formula, have to be replaced by a prior distribution.
AMC_bayes <- bayesian.vam( Time & Type ~ (ARA1(~Beta(1.08,0.108)) | Weibull(~NonInform(),~Gamma(32,0.097))),data=AMC_Amb) coef(AMC_bayes,profile.alpha=TRUE) summary(AMC_bayes)
R> AMC_bayes <- bayesian.vam( Time & Type ~ (ARA1(~Beta(1.08,0.108)) | Weibull(~NonInform(),~Gamma(32,0.097))),data=AMC_Amb) R> coef(AMC_bayes,profile.alpha=TRUE) [1] 4.083145e-06 2.986909e+00 8.558140e-01 R> summary(AMC_bayes) Initial parameters (by MLE): 1.304e-07, 3.102, 0.8981 (Mean) Bayesian estimates: 4.083e-06, 2.987, 0.8558 (SD) Bayesian estimates: 2.48e-05, 0.46, 0.1101 (0.025-Quantile) Bayesian estimates: 1.212e-09, 2.127, 0.5487 (0.975-Quantile) Bayesian estimates: 3.065e-05, 3.93, 0.9776 (Number) Bayesian estimates: 70371, 56970, 13401 Metropolis-Hasting acceptation rates: 0.7819, 0.633, 0.1489
The additional argument profile.alpha=TRUE in the coef function means that the prior distribution of the scale parameter of the Weibull distribution, denoted α, is not used. In fact, it is often difficult to find arguments to justify the prior distribution of α. But, for any value of the other parameters, we can explicitly compute the value of α that maximizes the likelihood. That can be used to profile the estimation method such that it does not depend on α (for more details, see F. Corset, L. Doyen and O. Gaudoin 2012: "Bayesian analysis of ARA imperfect repair models" Communications in Statistics - Theory and Methods 41, pp. 3915-3941).
And we can easily plot the posterior distribution of the parameters and for example of the maintenance efficiency parameter (the 3rd parameter of the output of the coef function).
hist(AMC_bayes,3,main="Posterior distribution of CM effect paramter")
We can also plot the posterior estimator of the different plots provided in the VAM package. The 95% credibility band is also represented on the plots.
plot(AMC_bayes,"I")
Preventive Maintenance
Real case study
Let us now consider a data set constitutes of independent and identical systems submitted to both Corrective and Preventive Maintenance actions. This data set have been proposed by M.L.G. de Toledo and Co (2016: "Optimal periodic maintenance policy under imperfect repair: A case study on the engines of off-road vehicles" IIE Transactions 48(8), pp. 747-758) and corresponds to the successive PM and CM times of 141 diesel engines of trucks used in Brazilian mining.
data(OREngines) head(OREngines) tail(OREngines)
R> data(OREngines) R> head(OREngines) System Time Type 1 1 18315 1 2 1 32133 -1 3 1 50934 -1 4 2 16137 1 5 2 34722 1 6 2 53990 -1 R> tail(OREngines) System Time Type 255 137 9972 -1 256 138 5218 -1 257 139 4497 -1 258 140 3330 -1 259 141 2205 -1 260 141 2614 -1
You can click on the heading of the previous R session and notice that
the data frame has 3 columns, the first column specifies to which system each event refers to,
some events have a type 1 and then correspond to PM times.
We can then compute the MLE, for example in the case of
log-linear distribution for the first time to failure of the new unmaintained system,
ARA1 CM,
ARA∞ PM.
ORE_mle<- mle.vam( System&Time&Type ~ (ARA1(0.5) | LogLinear(1,1e-5)) & (ARAInf(0.5)),data=OREngines) (ORE_est<-coef(ORE_mle))
R> ORE_mle<- mle.vam( System&Time&Type ~ (ARA1(0.5) | LogLinear(1,1e-5)) & (ARAInf(0.5)),data=OREngines) R> (ORE_est<-coef(ORE_mle)) [1] 1.183449e-05 1.732860e-04 6.781989e-01 9.349510e-01
Consequently, PM are estimated to be efficient and to reduce of 93% the virtual age. CM are estimated to be less efficient and to reduce of 68% the supplement of age accumulated since the last maintenance.
Using the logLik.mle.vam function, we can compute asymptotic marginal confidence intervals for the parameters.
Hessian<-logLik.mle.vam(ORE_mle,with_value=FALSE,with_hessian=TRUE) HesInverse<-solve(-Hessian,tol=1e-24) B<-matrix(c(rep(-1,4),rep(1,4)),ncol=2,dimnames=list(c('alpha','beta','rhoMC','rhoMP'))) ORE_est+qnorm(1-0.05/2)*sqrt(diag(HesInverse))*B
R> Hessian<-logLik.mle.vam(ORE_mle,with_value=FALSE,with_hessian=TRUE) R> HesInverse<-solve(-Hessian,tol=1e-24) R> B<-matrix(c(rep(-1,4),rep(1,4)),ncol=2,dimnames=list(c('alpha','beta','rhoMC','rhoMP'))) R> ORE_est+qnorm(1-0.05/2)*sqrt(diag(HesInverse))*B [,1] [,2] alpha 7.580448e-06 1.608853e-05 beta 1.456661e-04 2.009060e-04 rhoMC 5.637583e-01 7.926395e-01 rhoMP 8.103849e-01 1.059517e+00
Then, we can also plot plug-in estimator, for example, of the intensity for the 24th system. Vertical dashed lines correspond to PM times.
plot(ORE_mle,"i",system.index = 24)
PM policy
For simulating the maintenance times of a system subjected to PM actions, the user has to define when PM times are done. That can been specified with a PM policy in the VAM formula. The available PM policies are:
periodic PM times, denoted
Periodic,PM at predetermined times, denoted
AtTimes,PM done when the failure intensity reaches a predetermined level, denoted
AtIntensity,PM done when the virtual age reaches a predetermined level, denoted
AtVirtualAge,PM done when the conditional failure intensity reaches a predetermined level, denoted
AtFailureProbability.
simCMPM<-sim.vam( ~ (ARA1(.9) | Weibull(.001,2.5)) & (ARAInf(.4) | AtIntensity(0.23))) simData=simulate(simCMPM,Time>(RightCensorship=270)) head(simData) tail(simData)
R> simCMPM<-sim.vam( ~ (ARA1(.9) | Weibull(.001,2.5)) & (ARAInf(.4) | AtIntensity(0.23))) R> simData=simulate(simCMPM,Time>(RightCensorship=270)) R> head(simData) Time Type 1 15.74846 -1 2 34.55304 1 3 42.70481 1 4 45.54287 -1 5 53.41083 1 6 54.42935 -1 R> tail(simData) Time Type 56 244.4090 1 57 252.5607 1 58 260.7125 1 59 263.3309 -1 60 267.1955 -1 61 270.0000 0
For example, with the previous instructions, PM are done when the failure intensity reaches the level 0.23.
plot(simCMPM,"i",to=simData$Time[20]) abline(h=0.23,col="purple",lty=3)
Combined PM policies
With VAM, we can also combine several PM effects following different PM policies.
simCMPM<-sim.vam( ~(ABAO() | Weibull(.001,2.7)) & (ARAInf(.6)+ARAInf(-.2)+AGAN() | Periodic(12,prob=c(0.6,0.4))*AtIntensity(0.6))) simulate(simCMPM,20)
R> simCMPM<-sim.vam( ~(ABAO() | Weibull(.001,2.7)) & (ARAInf(.6)+ARAInf(-.2)+AGAN() | Periodic(12,prob=c(0.6,0.4))*AtIntensity(0.6))) R> simulate(simCMPM,20) Time Type 1 8.042624 -1 2 12.000000 2 3 13.201639 -1 4 15.705331 -1 5 21.613929 3 6 24.000000 1 7 36.000000 1 8 48.000000 1 9 59.541799 -1 10 60.000000 2 11 61.366679 3 12 72.000000 2 13 73.937744 -1 14 77.415553 -1 15 80.586810 -1 16 83.014405 -1 17 83.253944 3 18 84.000000 2 19 92.352679 -1 20 96.000000 1
You can click on the heading of the previous R session and notice that there is 4 different maintenance types.
Type -1 corresponds to ABAO CM (stars on the abscissa axis on the following plot), which have no effect on the failure intensity.
Type 1 corresponds to efficient ARA∞ PM (vertical red dashed lines on the following plot), which reduce of 60% the virtual age.
Type 2 corresponds to harmful ARA∞ PM (vertical green dashed lines on the following plot), which increase of 20% the virtual age. ARA∞ PM (type 1 and 2) are done periodically every 12 unit of times. The fact that PM is efficient or harmfull, is determined randomly: the PM has 60% of chance to be efficient.
Type 3 corresponds to renewing of the system (vertical blue dashed lines on the following plot), which are done when the failure intensity reaches the level 0.6.
plot(simCMPM,"i") abline(h=0.6,col="purple",lty=3)
PM policy depending on a model different from the one used for simulation
In the PM policies AtIntensity, AtVirtualAge or AtFailureProbability, PM are done when an indictor (i.e. the failure intensity, the virtual age or the conditional probability of failure) reaches a predetermined level. The computing of this indicator depends on the model and their parameters. By default, VAM uses the same model for simulating failure times and computing the PM threshold indicator. But we can also use different models or different parameters values. In the following planPMmod corresponds to the model used to compute the failure intensity threshold: ARA∞ PM are done each time this failure intensity (the black one) reaches the level 0.35. The formula in simCMPM describes the model used to simulate failure times. It corresponds to the blue failure intensity on the following plot.
planPMmod <- model.vam(Time & Type ~ (ARA1(.87) | Weibull(.0015,2.6)) & (ARAInf(.44))) simCMPM<-sim.vam( ~ (ARA1(.9) | Weibull(.002,2.5)) & (ARAInf(.25) | AtIntensity(0.35,planPMmod)) ) simData<-simulate(simCMPM,30) update(planPMmod,simData) plot(planPMmod,"i",ylim=c(0,0.5)) plot(simCMPM,"i",col="blue",add=TRUE) abline(h=0.35,lty=3)
Help
For more detailed informations about VAM functionalities you can read the R help delivered with the VAM package.
Further developments
We plan to add to VAM the following functionalities.
Adding some user-friendly functionalities in order to specify the names of the different estimated parameters, to compute directly parameter confidence intervals, to facilitate the adding of parameter constraints in the likelihood maximization, to simulate the future of a given data set, to make available the values of the different plotted characteristics, ...
Taking into account left censoring.
Considering semi parametric estimation methods (non-parametric estimation of the first time to failure distribution and parametric estimation of the different maintenance efficiencies parameters).
Considering observed and non-observed heterogeneity.
...
Related Publications
The package has been presented in different conferences:
- L. Doyen et R. Drouilhet : "VAM, outil logiciel adaptatif et libre pour prendre en compte l’effet des maintenances et du vieillisement". In λμ 20 - 20ème Colloque National de Maîtrise des Risques et Sûreté de Fonctionnement, Saint-Malo, France, October 2016.
- L. Doyen et R. Drouilhet : "VAM, an R package for maintenance and aging models". In MMR 2017 - 10th International Conference on Mathematical Methods in Reliability, Grenoble, France, July 2017.
- L. Doyen et R. Drouilhet : "VAM an R package to take into account the effect of maintenance and ageing". In ENBIS-17 - 17th Annual Conference of the European Network for Business and Industrial Statistics, Naples, Italia, September 2017. (Invited Conference).
- L. Doyen, R. Drouilhet : "A generic framework for recurrent event data based on virtual age models and implemented in the R package VAM". In CMStatistics 2017 - the 10th International Conference of the ERCIM WG on Computational and Methodological Statistics, London, United Kingdom, December 2017. (Invited Conference).
We are also finishing a journal paper describing the mathematical tools used to develop MARS. And we are writing another journal paper describing a real case data study with VAM.
Contacts
If you want to share your VAM experience with the author of the package you can contact:
- L. Doyen (), Univ. Grenoble Alpes, laboratory LJK.
- R. Drouilhet (), Univ. Grenoble Alpes, laboratory LJK.
We are also interested by your comments and ideas for further developments. VAM has been extensively tested, but if you notice some problems do not hesitate to contact us.