Reliability of repairable systems with Non-Central Gamma frailty

Maintenance actions on industrial equipment are essential to reduce expenses associated with equipment failures. Based on a well-fitted model, it is possible, through the estimated parameters, to predict several functions of interest, such as the cumulative average and reliability functions. In this paper, a new frailty model is proposed to analyze failure times of repairable systems subject to unobserved heterogeneity actions. The Non-Central Gamma distribution is assumed to the frailty random variable effect. The class of minimal repair models for repairable systems is explored considering an approach that includes the frailty term to estimate the unobserved heterogeneity over the systems’ failure process. Classical inferential methods were used to parameter estimation and define the reliability prediction functions. A simulation study was conducted to confirm the properties expected in the estimators. Two real-world data known in literature were used to illustrate the estimation procedures and validate the proposed model as a viable alternative to those already established in the literature. The results obtained highlight the potential of our proposed approach, particularly for industries dealing with such systems, where unquantifiable factors may impact equipment failure times.


Introduction
The reduction of expenses with actions for the maintenance of machinery(ies) is one of the main issues investigated in many industries that has led to the development of a huge literature dedicated to the so-called reliability theory.With technological advances and the increasing complexity of systems, the study of their reliability has been of great interest, especially for repairable systems.When considering a system, be it a machine, electronic equipment or software, it is said to be repairable when after a failure its activity can be satisfactorily resumed through repair, without the need for its replacement as a whole, according to the definition given by Ascher & Feingold (1984).Otherwise, a system is non-repairable when after its single failure it is discarded.Therefore, after a failure occurs, it is possible for the repairable system to be restored to operation condition by a repair action.
The study of repairable systems has wide application in contexts where the presence of recurring events is identified, such as in industry, where consecutive system failures are often observed.Investments in proper maintenance and repair to extend the lifetime of systems are essential, because in many cases repair is more economically feasible than replacing one or more components.The main challenge when modeling data from repairable systems is how to explain the effect of a repair action taken immediately after a failure.It is generally assumed that repair actions are instantaneous and the repair time is negligible (de Toledo et al., 2015).
Among the usual models of repairable systems, the Minimal Repair (MR) model, known in the literature as As Bad as Old (ABAO), is widely used due to the speed of repair and the low associated cost.This repair focuses on correcting only the component originating the failure, consequently leaving the system in the same working condition as before the failure, according to assumptions discussed in different works such as Barlow & Proschan (1996), Wang (2002), and de Toledo et al. (2015).
The statistical modeling usually used for these systems is done by counting process, which is a stochastic process that counts the number of occurrences of an event of interest in the time interval [0, t].We can classify the counting process as Homogeneous Poisson Process (HPP) and Non-Homogeneous Poisson Process (NHPP), being that for the latter, the intensity function is not constant.A particular parametric NHPP case widely used in the literature and that will be used in this study is the Power Law Process (PLP) proposed by Crow (1975), which is attractive for the ease of interpretation of its parameters and simplicity of implementation.
In the literature we found some methods for modeling data from repairable systems that do not incorporate the heterogeneity between the systems, and with this it is not possible to measure the most fragile, i.e., those that are more prone to fail (D' Andrea et al., 2017;Vaupel et al., 1979).It is relevant to verify how the heterogeneity of each system can impact on the main measures being investigated, such as increased failure rates (Proschan, 1963).
Frailty models are an alternative to estimate the unobserved heterogeneity of systems and are characterized by incorporating a random effect, either additively or multiplicatively, introduced by Vaupel et al. (1979).In the recent reliability literature, there has been increasing interest in data analysis of repairable systems incorporating frailty, with different parametric and non-parametric distributions assumed for frailty.
The choice of distribution for an unobservable random variable (frailty) is a problematic part when we are in an approach that considers this term in the model.Incorrect specification of the frailty distribution can lead to errors in the estimates of the parameters of interest, as is exposed in Almeida et al. (2020).The study considered a non-parametric approach to modeling the frailty density.By using Bayesian methods, an a prior Dirichlet process mixture is assigned to frailty, which in turn is more flexible and provides consistent estimates for the PLP model, as well as insights into heterogeneity among systems (Ascher & Feingold, 1984).
In particular, Slimacek & Lindqvist (2016) present a new method that seeks to avoid numerical problems often encountered in the estimation of models that incorporate frailty, in which the proposal is to estimate the parameters of a NHPP process with unobserved heterogeneity that does not require parametric assumptions.In addition, examples are illustrated involving the PLP compared to the Gamma frailty model and the model that does not consider the frailty term.
On the other hand, according to Wienke (2010), a parametric specification of the frailty distribution is a matter of mathematical convenience.Therefore, assigning a mathematically tractable distribution to the frailty variable with adequate properties can result in suitable models to estimate unobserved heterogeneity.In this sense, due to its mathematical simplicity, the Gamma distribution has traditionally been used as a frailty distribution in works in the areas of survival and reliability analysis, as also discussed by Wienke (2010).
In order to propose a distribution for the frailty term that is mathematically convenient as the Gamma distribution, which is used in most studies, Tomazella (2003) concluded that the inverse Gaussian distribution has properties equally simple to the Gamma distribution.Additionally, simulation studies were performed to infer on the parameters of interest and a Bayesian approach was proposed considering multiplicative and additive frailty models.
In a recent work, to deal with life insurance data, Onchere et al. (2021) propose a more flexible distribution for the frailty term, the Non-central Gamma distribution.According to these authors, this distribution can represent time-varying frailties and fits well to real life assurance data.In addition, they considered the Weibull, log-logistic and log-normal baseline risk functions and their results indicate that the Non-central Gamma-Weibull frailty model is adequate for the data in question.
Although the work of Onchere et al. deals with non-recurring events, the idea of a distribution capable of representing time-varying frailties is particularly attractive for models with recurring events.It is reasonable to think about the existence of unobservable factors that can influence the occurrence of recurring events differently in different periods of time.Furthermore, the Noncentral Gamma distribution is flexible and also presents mathematical facilities such as a closed form for the Laplace Transform.
In view of the above, it is possible to realize the relevance of incorporating the frailty term with a different distribution in the modeling of repairable systems, seeking for more flexible assumptions that better fit different scenarios.In this way, this paper will explore the models for repairable systems, submitted to minimal repairs with the effect of frailty, where times will be modeled by the PLP process, while for frailty the Non-central Gamma distribution will be considered.So far the approach with this frailty is unknown in the literature in a perspective that considers recurrent events and this is a relevant contribution of our work.
This paper is divided as follows.Section 2 reviews the Non-central Gamma distribution, models for repairable systems under minimal repair, and frailty modeling.In Section 3, the Non-central Gamma frailty model is introduced.In Section 4, some properties are verified through a simulation study.In Section 5, applications to real data are made.Finally, in Section 6, final considerations are raised.

Background
In this section we present the theoretical background used as a basis for our proposed Noncentral Gamma (NCG) frailty model for repairable systems.Initially, in Section 2.1, we define the NCG distribution and outline some of its key features.In Section 2.2, we present the main ideas of modeling failure times for systems subject to minimal repairs.Concluding our review in Section 2.3, we provide a brief characterization of multiplicative frailty models.

Non-central Gamma distribution
According to Knüsel & Bablok (1996) and de Oliveira & Ferreira (2013), the NCG distribution can be defined as a mixture of a discrete distribution and a continuous distribution, where the first distribution is a Poisson and the second is a central Gamma.Onchere et al. (2021), in turn, define the NCG distribution as a convolution generated by the sum of Gamma distributions with Poisson weights.Summarizing the ideas and definitions presented by these authors, the NCG distribution can be defined in a general way as a mixture of a Gamma(a + m, b) with weights Poisson(φ).In this way, let Z be a random variable with a NCG distribution, the probability density function (pdf ) of Z is defined by where φ > 0 is called the non-centrality parameter and a ≥ 0 and b > 0 comes from the Gamma distribution.
The NCG distribution is generally asymmetrical and its skewness is controlled by the noncentrality parameter.Figure 1 shows different behaviors of the density function of the NCG distribution for different parameter values.In this figure, parameter a has a fixed value of zero, since this will be an assumption adopted later in this work.A relevant characteristic of the NCG distribution is that it has a closed form for the Laplace Transform, which is given by As a direct consequence, we can easily obtain the expectation and variance of the NCG distribution, which are respectively given by E The calculations for obtaining the Laplace Transform, the expectation and the variance of the NCG distribution are presented in Appendix 1.

Review of minimal repair models for repairable systems
Minimal repair models are already widely known and discussed in the reliability literature.The tool of theses studies is the time elapsed until the occurrence of a failure, which can be represented by a random variable T. In the context of repairable systems, failure events are recurrent, which leads to obtaining a non-decreasing sequence of positive random variables 0 < T 1 ≤ T 2 ≤ . . .which represent the failure times.
Let N(t) be a number of observed failures in the interval [0, t].Under the MR assumption, {N(t), t ≥ 0} is considered a NHPP, that is, a counting process with independent increments, N(0) = 0 and completely characterized by a non-constant intensity function λ(t) (Gilardoni & Colosimo, 2007).Furthermore, for any t, N(t) has a Poisson distribution with mean Λ(t), where Λ(t) = E[N(t)] = t 0 λ(u)du is the cumulative intensity function.Assuming that the intensity function λ(t) can be described by a parametric model with a vector of parameters µ, the failure time history of one or multiple identical systems under MR models can be used to estimate this parameter vector.In a classical inference approach, the likelihood method can be used to obtain the maximum likelihood estimators for the parameters (Tomazella, 2003).Assume that k identical repairable systems are under observation and let 0 < t i,1 < t i,2 < • • • < t i,n i be their failure times, where t i,j is the j-th failure time of the i-th system, with i = 1, . . ., k and j = 1, . . ., n i .Note that n i is the total number of observed failures of the i-th system and consider that for each system i there is still an observation truncation time t * i with t i,n i ≤ t * i .So, the general likelihood function for the MR model considering the k identical systems is given by Furthermore, the intensity function can be used to compute predictors of the reliability function.Let T i,n = t i,n be the last observed failure time of the system i and t = T i,n+1 -t i,n the time until the next failure.Considering the parametric MR model with vector of parameters µ, the general reliability at time t for the system i is given by = exp - Note that the function (4) expresses the probability that the system will run without failure for a period of time greater than t after the last failure time t i,n , given the failure history of the process.
All the information about MR models that we present here is known in the reliability literature and, particularly, the estimation and prediction of reliability are general for whatever the parametric form of λ(t) is.Among the possible forms, we will present and use in this work the PLP.Introduced by Crow (1975), the PLP is a parametric model widely used in the literature for NHPP, in which its intensity and cumulative intensity functions are respectively given by where η > 0 is the scale parameter and β > 0 is the shape parameter.The PLP is commonly assumed for MR models since its parameters have a physical interpretation in the following sense: the parameter η indicates the time during which exactly one failure is expected to occur, that is, E[N(η)] = 1, while the parameter β indicates whether the system is deteriorating (β > 1) or improving (β < 1) (Dias De Oliveira et al., 2013).
To obtain the likelihood and reliability prediction functions considering the PLP parametric model, just replace λ(t) and Λ(t) in the equations ( 3) and ( 4) by the intensity and cumulative intensity functions given in (5).In this case, the model vector parameter will be µ = (β, η).This procedure and respective observations can be seen in

Frailty models
As stated in the introduction of this work, frailty models incorporate in the intensity function a latent random variable Z that captures random effects.This variable is drawn from a distribution f (z; θ) and, in this study, it will be included in the intensity function in a multiplicative manner.In this sense, frailty models are considered as natural extensions of the Cox model.Disregarding the presence of covariates, a multiplicative frailty model for repairable systems is described by where λ 0 (.) is the baseline intensity function and z is the frailty term.
The choice of f (z) is an important point to consider and some authors have already discussed the impacts of adopting different distributions for the frailty term, e.g.Hougaard (1984) and Tomazella (2003)).Any continuous or discrete distribution with positive support, finite mean and unknown finite variance can be chosen as a frailty distribution.However, some conditions must be satisfied for the baseline intensity function to be identifiable, such as E[Z] = 1 and Var(Z) = θ (Elbers & Ridder, 1982).These restrictions lead to the following interpretations of the frailty term: if z i > 1, the failure intensity is increasing; if z i < 1, the failure intensity is decreasing; and if z i = 1, the model reduces to the baseline intensity function.

NCG frailty model for repairable systems
In this section we present the proposed NCG frailty model for multiple repairable systems under MR.In addition, we present the likelihood function used to estimate the model parameters as well as the reliability prediction function.
As mentioned earlier, to avoid identifiability problems in the model with frailty following a Z distribution, the model requires that its expectation is equal to 1 and its variance is constant.To guarantee this assumption when defining that the frailty Z has an NCG distribution and following the idea presented in Onchere et al. (2021), it is necessary to assume that the parameter a in the equation ( 1) is equal to zero.Furthermore, considering Var[Z] = θ the parameters b and φ of equation (1) will be reparametrized as b = θ/2 and φ = 2/θ, so that the expectation is such that E[Z] = 1.In this case, the pdf of the frailty Z in equation ( 1) can be rewritten as Consequently, the Laplace Transform given in equation ( 2) can be rewritten as .
Obtaining a closed form for the Laplace Transform is a particularly interesting result for the frailty proposal with NCG distribution, especially for the calculation of the likelihood function.This is because the likelihood function must be obtained in a non-conditional way to frailty Z and the integrals used in these calculations are not trivial.However, if Z has a distribution that assumes a closed form for the Laplace Transform, this problem is easily solved.
According to Wienke (2010), the unconditional reliability function is given by R(t) = L(Λ 0 (t)) and the intensity function unconditional the frailty term is given by λ(t) = -λ 0 (t) Therefore, the unconditional reliability function for the NCG frailty model is given by and the unconditional intensity function and the unconditional cumulative intensity function for the NCG frailty model are respectively given by where λ 0 (t) and Λ 0 (t) are the baseline intensity function and the baseline cumulative intensity function, respectively.The functions in ( 7) can be interpreted as functions of population intensity and cumulative population intensity that depend directly on the baseline functions and the variance of the assumed NCG frailty.The failure intensity function in (7) can assume unimodal forms and not just an increasing or decreasing form (as the PLP intensity function without frailty, for example), which is especially useful for modeling situations where the failures can be concentrated in one period and occur less frequently in other periods.
Having defined the frailty model functions, we now want to use observed data to estimate the model parameters.Consider k repairable systems under observation for a period t * i , i = 1, . . ., k so that t i,1 , t i,2 , • • • , t i,k are the failure times of the i-th system and t * i is the truncation time of the observation of this system.Let us assume that these systems undergo minimal repairs after each failure, so on the classical approach the population parameters can be estimated through the maximum likelihood method.
Let us assume that each failure process follows a PLP whose failure intensity and cumulative intensity functions are given by (5) and that there is an unobserved heterogeneity between the failure times modeled by the NCG distribution given in (6).Therefore, the vector of parameters to be estimated in the model is µ = (β, η, θ), where β and η are the parameters of the PLP model and θ is the variance of the NCG frailty.In this case, the likelihood function in (3) can be rewritten as and consequently, the log-likelihood function is given by where N = k i=1 n i is the total number of observed failures across all systems.Given the complexity of the expression (9) and its partial derivatives, it is not possible to obtain the maximum likelihood estimators analytically.Therefore, numerical methods can be used to circumvent this problem and thus obtain the estimates.The confidence intervals of the estimated parameters can be constructed from the asymptotic theory based on the properties of the maximum likelihood estimators and the Normal distribution.
Finally, based on the reliability predictor function for the classical MR model presented in (4), we can define the reliability predictor function for our proposed NCG frailty model.Given the maximum likelihood estimates ( β, η, θ) of the NCG frailty model parameters and let t = T i,n+1t i,n be the time until the last observed failure time t i,n to the next failure time T i,n+1 , the reliability function is given by where the cumulative intensity function Λ(t) is the cumulative intensity function of the NCG frailty model given in (7).

Simulation study
In this section, we conduct a simulation study to assess the efficiency and consistency performance of the maximum likelihood estimator (MLE) defined by maximizing the log-likelihood function in the model presented in (8).The main objective is to show that the asymptotic properties of the proposed maximum likelihood estimators are observed in samples of a simulated population of failure times with frailty effects induced by the NCG distribution.
To evaluate the consistency and efficiency of the MLE of the NCG frailty model with PLP base intensity we used three metrics, based on similar studies carried out in the reliability literature: the mean relative error (MRE), the square root of the mean standard error (RMSE) and the coverage probability (CP) of the confidence interval (CI).These three metrics are respectively defined by where: • M is the number of samples to be simulated; • µ represents one of the parameters to be estimated (β, η or θ) and μm is its estimate in the m-th sample; • a m = μm -1.96 × SE(μ m ) is the lower bound of the 95% CI of the m-th estimate; • b m = μm + 1.96 × SE(μ m ) is the upper bound of the 95% CI of the m-th estimate; • I is the indicator function and SE(μ m ) is the standard error of the m-th estimate.
By this approach, it is expected that as we increase the sample sizes of observed failure times, the MRE will asymptotically approach 1, the RMSE will approach 0 and the CP will reach an adequate percentage of coverage of the nominal values of the parameters by the asymptotic CI (here we are adopting 95%).The increase in the sample will occur by increasing the number of observed systems.In this case, the number of failures observed in each system is random, however, it is reasonable to expect that the more systems observed, the more failure times are recorded.
To generate the failure times for each system, the population unconditional cumulative intensity function of the NCG frailty model defined in (7) will be used.More specifically, let t j and t j+1 be two consecutive failure times such that t j+1 = t j + x.As described in Brito et al. (2022), the cumulative distribution function F(x) for the elapsed time X between two consecutive failures is F(x) = 1 -P[N(t j + x) -N(t j ) = 0].Under the MR assumption, the counting process N(t) is an NHPP and as discussed in Section 2.2, the cumulative distribution function can be rewritten as where Λ(t) is the cumulative intensity function of the NCG frailty model, given in (7).Thus, to generate each process failure time, simply generate the time between two consecutive failures as a solution to the equation F(x) = u, where u is an observation of the Uniform(0, 1) distribution.This procedure will be repeated until the desired number of failures are generated or until the predefined truncation time is reached.Furthermore, considering k distinct s, this entire procedure will be repeated independently k times.
The algorithm below describes the process of generating failure times for k systems, defining a truncation time t * of the observations for each system.
Different combination scenarios of model parameters were considered for the simulation study, as well as different sample sizes based on the number of systems.Due to the unimodal characteristic of the NCG frailty model's failure intensity function, for each parameter scenario the truncation time was defined as the time t * such that λ(t * ) < 0.001, since the model loses influence in the generation of failure times when the intensity asymptotically approaches zero.Two fixed values were defined for the parameter β = (0.8, 1.2), to consider the two cases where the PLP baseline intensity function is increasing and decreasing.For each fixed β, three values were defined for the parameter θ = (0.05, 0.1, 0.2), that include different variances for the assumed NCG frailty.For each combination of parameters β and θ, a parameter η was defined to complete the scenario.Here, we will present the results of six different scenarios that represent the aforementioned model variations.
For each defined scenario, 1,000 samples were generated following the previously defined algorithm and, for each sample, the β, η and θ estimates of the parameters β, η and θ, respectively, were obtained.From there, the previously defined MRE, RMSE and CP metrics were obtained and the results are summarized and shown in Figures 2 and 3. Additionally, all results are reported in Tables 5-10 in Appendix 2.
From Figures 2 and 3, it is possible to see that, in general, the expected behaviors of the three metrics presented are verified in all determined scenarios.It is noticeable that increasing the sample size by increasing the number of observed systems results in CP increasingly close to the nominal value of 95% and RMSE increasingly close to zero for all parameters.It is clear and expected that specifically for the parameter η the RMSE value converges more slowly to zero, given the magnitude of the nominal values of this parameter; however, observing the MRE graphs, the expected convergence of this parameter to the nominal value is verified.In scenarios where there is a greater variance of frailty (θ = 0.2), it can be seen from the MRE graphs that the parameter θ can be slightly underestimated by the model; however, the interval estimates encompass the true nominal values as observed by the high percentage coverage on the CP graphs.
We can conclude through the analysis of the simulation studies that the estimators for the proposed model perform well, especially for large samples.

Application to real data
In this section we present two applications revisiting two known datasets in the reliability literature.The idea is to verify whether in these data it is possible to identify the presence of unobserved heterogeneity in the failure times of the systems and to verify whether the proposed NCG frailty model is adequate for such data.
In both examples the parameters were obtained using the maximum likelihood method through the likelihood and log-likelihood functions defined in ( 8) and ( 9) respectively, and the confidence intervals were obtained using the asymptotic theory of the normal distribution, as mentioned in Section 3. In addition, we also used the well-known Gamma frailty model to compare the fits obtained with the proposed model in both datasets and we used the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) as selection criteria between these two models (for more details on these criteria, see, e.g., Burnham & Anderson (2004)).

Sugarcane harvester data
The first application deals with a set of failure times of 9 Brazilian sugarcane harvesters.This dataset was originally presented by Verssani (2018) and a subset of these harvesters was used in the frailty analysis under the assumption of minimal repair presented in D'Andrea et al. (2021).Specifically, the times of failure of the cutting blade system of these harvesters during the period of high sugarcane harvest are observed between the years 2014 and 2015, for a period of 200 days.Note that the truncation time of observations for all of them is the same equal to 200 days.
Table 1 shows the observed failure times of each harvester.Note that the truncation time of observations for all of them is the same equal to 200 * days.8), numerical and computational methods can be used to maximize this function and obtain the desired MLEs.This maximization procedure was performed using the software R (R Core Team, 2021), specifically using its optim function.Despite the algebraic complexity of the likelihood function, making it impossible to obtain the MLE analytically, this is not a major problem from a computational point of view, so these estimates were easily obtained.
As stated earlier, we also used the Gamma frailty model for repairable systems under minimal repair by D'Andrea et al. (2017) to estimate their parameters and compare with the results obtained by our NCG frailty model.The results of the parameter estimates for the both models are shown in Table 2.In addition, this table lists the maximum values of the log-likelihood function of each model, as well as the AIC and BIC values.It is immediately possible to see the similarity between the estimates of β, η and θ for both models.In terms of comparing the models, the AIC and BIC values indicate that the NCG frailty model is slightly superior to the Gamma model for this data set under minimal repair, since they presented lower values.Although the difference in the information criteria for the two models is small, this result shows that our proposed model for NCG frailty is a good alternative to the Gamma model, and may even present slightly more parsimonious estimates.
The point estimate of parameter β indicates that the systems are deteriorating over time once β > 1, however this interpretation may not be true since the lower limit of CI is less than 1.Regarding the frailty, despite small, a value greater than zero was obtained for the variance of the frailty variable ( θ > 0), indicating the existence of unobserved heterogeneity between the systems and their failure times.In this sense, this result is in line with the other results obtained for this data set in the literature, indicating the existence of unobserved effects on the failure time of the harvesters, since the existence of frailty effects was captured.
Finally, we use the equation (10) to estimate the reliability prediction function at the last observed failure time for each harvester.Figure 4 shows the reliability prediction function for the NCG frailty model related to each harvester of this study.It is possible to notice that all harvesters have almost zero reliability after 50 days of the last failure, given the history of each failure process.Furthermore, it is possible to notice that the median expected time to the next failure (i.e., the time such that the expected reliability is 0.5) is close to the estimated η = 14.4 value, for all harvesters.

Literature systems data
The second application deals with the failure times of the set of repairable systems presented in Mettas & Zhao (2005).In the original work, failure times of six systems with different observation truncation times are presented.However, for our application we will use the failure times of only three of these systems, since the other systems have a very low number of observed failures (one or two failures).
Table 3 shows the observed failure times of each of the systems.The times marked with " * " are the truncation times of the observation of each system.
The same procedures as in the previous example were repeated here.Table 4 shows the MLE obtained for the Gamma and NCG frailty models, as well as their estimates of maximum log-likelihood and AIC and BIC values.Once again, it is noticeable that the point estimates are very close for the    The 95% CIs for the MLE β, η and θ of the NCG frailty model are, respectively CI β = (1.04,3.38), CI η = (907.92,2763.25) and CI θ = (0.04, 0.41).In this case, we have more assertive conclusions about the estimates obtained.First, regarding the parameter β, both the point and interval estimates are greater than 1, indicating that the systems are deteriorating over time.Furthermore, the point and interval estimates of parameter θ are greater than zero, again capturing the existence of unobserved heterogeneity between systems and their failure times.
Once again, we use the equation (10) to estimate the reliability prediction function at the last observed failure time for each system in the studied dataset.The results of the estimated reliability prediction functions for the NCG frailty model are shown in Figure 5.The reliability probability for 2,000 days after the last failure is very close to zero for each system, given the entire history of its failure process.As in this case the estimate θ did not approach zero, the interpretation for the parameter η cannot be done as in the previous example, since the high variance of the frailty variable impacts the interpretation of this parameter in the PLP, as discussed previously.

Final remarks
In this article, the main focus was to incorporate and estimate the frailty term in a model for repairable systems subject to unobserved heterogeneity, under the assumption of minimal repair.For the studied model, the likelihood functions and the reliability prediction functions were derived.The times were modeled by the PLP process, while the NCG distribution was considered as distribution for the frailty term.
A simulation study considering different scenarios was carried out to verify the quality of the obtained estimators and their behavior.It was found that the results obtained are in line with expectations, since the estimates generated were close to the actual values of the parameters predefined in the study, returning expected values for the analyzed metrics (MRE, RMSE and CP).
Two applications were made to real data.We verified that the model has better performance according to the considered evaluation criteria even if they were close to the alternative model that assumes the Gamma distribution for the frailty term.Furthermore, the parameter estimates are in agreement with the results found in the literature in applications for the same datasets.
Furthermore, in each application it was possible to predict the reliability of each system, based on their failure process history and the obtained estimates of the model parameters.These results show that the approach proposed in this work is effective and can be applied in different industry situations With this, the application in real data showed that the proposed approach is capable of providing useful information for decision making in complex systems, such as forecasting until the first failure.In summary, the results from the simulation and applications on real data corroborate the assertion that the proposed approach constitutes a valuable tool for the analysis of reliability in repairable systems with NCG distribution frailty.
Therefore, as up to the present moment the approach with this frailty was unknown in the literature in a framework of repairable systems.The results obtained from theoretical discussions, applications to real data and the simulation study show that the model proposed in this work is relevant and presents a modeling alternative that clearly identifies the existence of unobservable factors, such as heterogeneity over the systems' failure times.It is worth remembering that such factors directly interfere with the failure intensity function of the systems, which is closely linked to the useful life of the systems.
Finally, this study presents an alternative propose when we are dealing with non-central distributions.We emphasize here the importance of one more option for frailty distribution among the models already existing in the literature.The studies developed here are not exhaustive and we emphasize the importance of investigating other approaches, such as other methods of estimating the parameters of interest.
As a consequence, the first and second derivatives of the Laplace transform are respectively given by Finally, the expectation and variance of the random variable Z with NCG distribution can be easily obtained by derivatives of the Laplace Transform (see, e.g.Wienke, 2010), as follows:

Appendix 2. Simulation study results tables
In this Appendix we present the detailed results of the simulation study discussed in Section 4. Tables 5-7 below refer to the results summarized in Figure 2, while Tables 8-10 refer to the results summarized in Figure 3.

Figure 1 .
Figure 1.(a) NCG pdf for fixed b = 3 and variable φ values, and (b) NCG pdf for fixed φ = 50 and variable b values.Both cases with a = 0 fixed.

Figure 4 .
Figure 4.Estimated reliability functions at last failure times tn, for each harvester in the data set, under the fitted NCG frailty model.

Figure 5 .
Figure 5.Estimated reliability functions at last failure times tn, for each system in the data set, under the fitted NCG frailty model.

Table 1 .
Harvester cutting blade system failure times *Substituting the failure times t i,j observed in Table1into the likelihood function shown in (

Table 2 .
Estimation results for the frailty NCG and Gamma models applied to sugarcane harvesters data

Table 3 .
Systems failure times

Table 4 .
Estimation results for the frailty NCG and Gamma models applied to literature systems data