VAM an R package for Virtual Age Models


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.

The strengths of VAM are:

VAM has been developed by

both members of the laboratory LJK, of Univ. Grenoble Alpes.


Before being able to use VAM, R has to be install on your computer. Then, download the archive file corresponding to your operating system.

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:

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.

Loading R package 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.

AMC Abassador data set
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,

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.

ARA$_{\infty}$ (also called Kijima type II)

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.

ARA$_{1}$ (also called Kijima type I)

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")
QR (also called geometric process)

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")
GQR (also called extended geometric process)

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")
mixing GQR and ARA

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")


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)))
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.


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.

Right censoring in simulation
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.

Simulation of iid trajectories
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.

Maximum Likelihood Estimators
AMC_mle<-mle.vam(Time & Type ~ (ARAInf(0.6)|Weibull(1,3)),data=AMC_Amb)
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.


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.

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.

Bayesian estimators
AMC_bayes <- bayesian.vam( Time & Type ~ (ARA1(~Beta(1.08,0.108)) | Weibull(~NonInform(),~Gamma(32,0.097))),data=AMC_Amb)
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.


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.

Off Road Engines data set
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

We can then compute the MLE, for example in the case of

MLE for off road engines data set
ORE_mle<- mle.vam( System&Time&Type ~ (ARA1(0.5) | LogLinear(1,1e-5)) & (ARAInf(0.5)),data=OREngines)
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.

95% confidence intervals
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:

PM policy at intensity
simCMPM<-sim.vam(  ~ (ARA1(.9) | Weibull(.001,2.5)) & (ARAInf(.4) | AtIntensity(0.23)))
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.


Combined PM policies

With VAM, we can also combine several PM effects following different PM policies.

Combined 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)))
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.


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)) )


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.

Related Publications

The package has been presented in different conferences:

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.


If you want to share your VAM experience with the author of the package you can contact:

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.