Efficient semiparametric estimation of time-censored intensity-reduction models for repairable systems

The rate reduction models have been widely used to model the recurrent failure data for their capabilities in quantifying the repair effects. Despite the widespread popularity, there have been limited studies on statistical inference of most failure rate reduction models. In view of this fact, this study proposes a semiparametric estimation framework for a general class of such models, called extended geometric failure rate reduction (EGFRR) models. Covariates are considered in our analysis and their effects are modeled as a log-linear factor on the baseline failure rate. Unlike the existing inference methods for the EGFRR models that assume the failure data are censored at a fixed number of failures, our study considers covariates and time-censoring, which are more common in practice. The semiparametric maximum likelihood (ML) estimators are obtained by care-fully constructing the likelihood function. Asymptotic properties including consistency and weak convergence of the ML estimators are established by using the properties of the martingale process. In addition, we show that the semiparametric estimators are asymptotically efficient. A real example from the automobile industry illustrates the usefulness of the proposed framework and


INTRODUCTION
Many physical systems such as vehicles, commercial products, and manufacturing equipment are repairable in the sense that they are repaired rather than replaced in the presence of failures.Reliability evaluation of these systems is mainly based on modeling of the recurrent failure processes.Due to the complicated effects of the repair actions on the system health, however, modeling of the repairable systems is not an easy task, and it is much more difficult than modeling of the nonrepairable systems.
In the literature, a multitude of models have been proposed to account for the repair effects in the recurrent failures.Most of these models can be defined based on the intensity process.Consider a repairable system observed over time t ≥ 0, and let S 1 , S 2 , … be the failure time sequence of the system.Denote N * (t) as the number of failures/repairs in the time interval [0, t].Define  (t) = {N * (s) ∶ 0 ≤ s ≤ t} to be the filtration, that is, the history of the process up to time t.The intensity process v(t) ≡ v(t;  t− ) of the counting process N * (t) is defined as the conditional probability that an event occurs in [t, t + dt) given the filtration  t− , that is, Let  0 (t) be the baseline intensity function, and Λ 0 (t) = ∫ t 0  0 (u)du be the corresponding baseline cumulative intensity function.The intensity processes of many imperfect maintenance models are based on a modification of the baseline intensity  0 .When the repair effect is assumed to be as good as new (AGAN), {N * (t), t ≥ 0} is a renewal process and the intensity process is given by v . If the repair is as bad as old (ABAO), {N * (t), t ≥ 0} is a nonhomogeneous Poisson process and the intensity process is the same as  0 , that is, v(t) =  0 (t).In practice, it is more common that the repair effect is between AGAN and ABAO, which is known as imperfect repair.In terms of modeling of the imperfect repair, Kijima's Type I and Type II models (Kijima, 1989) represent the first attempts by introducing the concept of virtual age.The Kijima's models have been further extended in many follow-up studies with applications to optimal decisions in preventive maintenance.For an overview of these extensions, see Finkelstein (2008), Zhang and Xie (2017) and Alaswad and Xiang (2017).The Kijima's models treat the repair effect as a shift on the original survival function with a virtual age.In the literature, models with a change on the failure rate after repair have also been widely used for imperfect repairs.Among them, the Geometric Failure Rate Reduction (GFRR) model proposed by Finkelstein (2008) is a flexible model with a simple intensity process (t) =  N * (t−)  0 ( t − S N * (t−) ) , where the parameter  > 0 characterizes the repair effect.The empirical studies in Doyen et al. (2017) and Syamsundar and Kumar (2017) have revealed that the GFRR model outperforms most existing imperfect repair models when fitting to real datasets.To mitigate the explosion problem, for which  N * (t−) becomes too large when  is far from 1 or when N * (t) is large, Dauxois et al. (2019) proposed an Extended GFRR (EGFRR) model with the intensity process v(t) =  (N * (t−)+1)  0 ( t − S N * (t−) ) , where (k) is a known nondecreasing function such that (1) = 0.This model implies that the hazard function of the kth interoccurence time equals  (k)  0 (⋅).The EGFRR model degenerates to the GFRR model when (k) = k − 1.For the EGFRR model, Dauxois et al. (2019) have developed a semiparametric inference procedure assuming that the failure data are censored at a fixed number of failures.
In most maintenance applications, however, the observed repair process is time-censored because of an end-of-study date (Ye & Ng, 2014), which is also called administrative censoring in biostatistics.This time-censoring assumption is widely used in recurrence data; see Pena et al. (2001), Pena et al. (2007) and Rahman et al. (2014).In addition, the covariates may also play a critical role in affecting the failure rates; see Cox (1972), Wang and Xu (2010) and Qu et al. (2016).In view of this fact, we propose a semiparametric inference procedure for the recurrent failure data by considering the covariates and the time-censoring scheme.That is, in each experiment, we can only observe the failures/repairs before a censoring or follow-up time, which might incorporate the covariates.In specific, the intensity process considered is where  is a p × 1 vector and X = (x 1 , x 2 , ..., x p ) ⊤ is a p-dimensional covariate, such as the use rate, brand, and other observable features that vary from unit to unit.In this study, the maximum likelihood (ML) estimators and their asymptotic properties are thoroughly investigated.In particular, semiparametric efficiency properties, which are important in knowing how much the estimators can be improved, are established for the proposed estimators.The semiparametric efficiency is rarely discussed in the literature of statistical inference for the repairable systems, and our study bridges this gap by examining the efficiency of both the parametric and the nonparametric parts of the estimators.
The rest of the paper is organized as follows.In Section 2 we derive the likelihood and the ML estimators based on the time-censored repair data with covariates.Section 3 establishes the asymptotic properties including the consistency, the weak convergence and the semiparametric efficiency for the ML estimators.Comprehensive simulations are conducted in Section 4, and a real example from the automobile industry is used to illustrate the proposed inference procedures in Section 5.At last, Section 6 concludes the paper.

DATA AND LIKELIHOOD
Consider a system whose failure/repair process {N * (t), t ≥ 0} is governed by the EGFRR model (1) with a covariate vector X.By letting  = log(), the intensity can be written as ) .
Denote T j as the interoccurence time between the (j − 1)th and the jth failures, and S j = ∑ j l=1 T l as the occurrence time of the jth failure, j = 1, 2, • • •.For notational convenience, we stipulate T 0 = 0 and S 0 = 0. Let C denote the censoring or follow-up time of the system.The censoring mechanism is supposed to be conditionally independent of the failure process, that is, C and N * (t) are independent conditional on X for any t ≥ 0. Define Ñ(t) = N * (t ∧ C) and Ỹ (t) = I{C ≥ t}, where a ∧ b = min{a, b}, and I{A} is the indicator function that equals 1 if A is true, and 0 otherwise.The observed failure data from the system is {X, Ñ(⋅), Ỹ (⋅)}.We are interested in estimating the true parameters, denoted as  0 ,  0 and Λ 0 , from the data.

2.1
The likelihood Based on the above problem setting, Ñ(s) is a submartingale with respect to the filtration  s with compensator Ã(s; is the time from the last failure to time u.However, the compensator is a complicated function of the model parameters since the integrand contains a random argument R(u).Therefore, decomposition of this submartingale may not be useful in parameter estimation.In this study, we follow Pena et al. (2001) and introduce a two-dimensional counting process where Z(u, t) = I{R(u) ≤ t}.Here, N(s, t) is the failure count in the time window [0, s] whose interoccurence times are shorter than t.It is readily seen that N(s, t) = N(s) when s < t.Let  s = {X, Ñ(u), Ỹ (u+) ∶ 0 ≤ u ≤ s} be the event history up to time s.The following lemma summarizes the decomposition results for N(s, t), which can be easily verified by proposition 31.1 in Pena et al. (2000).
Lemma 1.For each fixed t, consider the counting process {N(s, t) ∶ 0 ≤ s ≤ } indexed by s as defined above.Then there exists a unique increasing right-continuous  s -predictable process {A(s, t;  0 ,  0 , Λ 0 ) ∶ s ≤ } given by such that N(s, t) − A(s, t;  0 ,  0 , Λ 0 ), denoted by M(s, t;  0 ,  0 , Λ 0 ), is a right-continuous local martingale with respect to the filtration  s for fixed t.Here, In the lemma, Y (s, t; ) can be roughly understood as the weighted count of subintervals in [0, s ∧ C] whose lengths are longer than t, where [0, s ∧ C] is partitioned by the mesh defined by all failure times within this time interval.The lemma above enables us to write the likelihood function.According to corollary II.7.3 in Andersen et al. (2012), the likelihood based on the observation {X, where we stipulate 0 0 = 1.Based on the definition of N(s, t) and A(s, t; , , Λ), the log-likelihood function can be derived as where Λ{x} denotes the pointmass of the baseline cumulative intensity function Λ(x), which equals (x) when Λ(x) is differentiable (Van der Vaart, 2000, pp. 425).By proposition 1 in Pena et al. (2001), the first term of the right-hand side of the above display can be expressed as

ML estimator
For notational convenience, we let  = (,  ⊤ ) ⊤ .Based on the log-likelihood in (4), the score function for the parametric part  is where Y (k) represents the kth derivative of Y with respect to , k = 1, 2, … .In order to derive the ML estimator for the nonparametric part, for a fixed Λ and each h(⋅) ∈  ∞ ([0, ]), which is the space of uniform bounded functions, we consider the parametric submodel Λ ,h indexed by  in a neighborhood of zero given by Notice that Λ ,t is nonnegative and monotone increasing when  takes its value in a neighborhood of zero.Intuitively, an optimal nonparametric estimator should maximize the log-likelihood among submodels going through it.We will first show that the nonparametric part is unique under the submodel, and demonstrate the asymptotic properties of the estimator in the next section.For this submodel, we can calculate the score function by ) .
Now, suppose we observe n i.i.d.copies of {X, Ñ(⋅), Ỹ (⋅)}, denoted as {X, Ñi (⋅), Ỹi (⋅)}, i = 1, 2, • • • , n.Let P n be the empirical measure based on the n i.i.d.realizations.By fixing , we can solve where N(, dx) = ∑ n i=1 N i (, dx), as with Fleming and Harrington (2005).By substituting the expression (6) into (4), the profile log-likelihood function for  = (,  ⊤ ) ⊤ can be obtained as Here we multiply the log-likelihood by n −1 for notational convenience when we consider the asymptotic properties.The following proposition shows that  n is concave.7) is continuous and concave.

Let
θ⊤ n = (γ n , β⊤ n ) be a solution to the estimating equation    n (, ) = 0.Because of the concavity, existing convex optimization packages can be used to find θn numerically and stably.Then the estimator of the nonparametric part is Λn (t; γn , βn ).

Time-varying covariates
In this section, we consider the case of time-varying covariates.The intensity function is given by where X(t) is the time-dependent covariate.As we have shown in Section 2.2, according to Cox (1972), Lin (1991) and Huang et al. (2010), the log-likelihood based on the data {X(⋅), Ñ(⋅), Ỹ (⋅)} is given by where the weighted count Y (s, t; , , X) is Therefore, we can again consider the family of submodels which simply replaces e  ⊤ X Y (, x; ) by Y (, x; , , X) in ( 6).The parametric estimator can then be readily obtained by maximizing the profile likelihood function.For the purpose of concise presentation, we use X to denote the covariate vector in the remaining of the paper.In the presence of time-varying covariates, the derivatives of Y (s, t; , , X) with respect to  and  can replace the derivatives of e  ⊤ X Y (s, t; ) in all the theoretical findings.

Failure censoring
In this section, we verify that the derived ML estimators in Section 2.2 remain valid when the data are failure censored, which is the setting in Dauxois et al. (2019).The authors considered a single recurrent event with imperfect maintenances in their paper.Assume that the number of imperfect maintenances of a repairable system is censored at m, which is a positive integer-valued random variable upper bounded by a fixed given number N. Further assume that the distribution of m does not involve the parameters  and Λ.Based on the observation (T 1 , • • • , T m ), the likelihood function can be written as This likelihood is equal to (4) by letting C = S m , where S m is the occurrence time of the mth failure.Therefore, the estimators θn and Λ(⋅) given in Section 2.2 are also ML estimators under the failure censoring mechanism.In Section 4, extensive simulations will be conducted to compare the proposed ML estimators with the estimators in Dauxois et al. (2019) based on the failure censored data.

ASYMPTOTIC PROPERTIES
In this section, we study the asymptotic properties of the estimators derived in Section 2, including the consistency, weak convergence, and semiparametric efficiency.

Consistency and asymptotic normality
Following the conventional notation in the large sample theory in survival analysis, we let , where X ⊗2 = XX ⊤ .Notice that S (1) n is the gradient of S (0) n and S (2) n is the Hessian matrix of S (0) n with respect to the parametric part (,  ⊤ ) ⊤ ∈ R p+1 .Further define and n (, t; , ) − s (k) (, t; , )|| → P 0 for every fixed  and  satisfying s (k) (, 0; , ) < ∞.Furthermore, we can obtain a result of uniform convergence with respect to  according to the following lemma.
Lemma 2. Consider a convex neighborhood  of  0 .Let M  and M  be the supremum of || and |||| in .Suppose sup t∈[0,] |X(t)| < M X , and the random variables defined below has finite second moment: where k = 0, 1, 2, and 3. Then We should remark that the conclusion in Lemma 2, that is, the last display, is usually treated as a regularity condition in the literature of survival analysis.Using the advanced tools of empirical process, the above lemma directly establishes these uniform convergence results, which will be useful in establishing the asymptotic properties of our estimators.Actually, under the conditions given in Lemma 2, it is easy to see that sup That is, s (k) (, t; , ) is uniformly bounded on [0, ] for k = 0, 1 and 2. In addition, notice that ES (0)  n (, ; , ) ≥ P {C ≥ , T 1 ≥ } ⋅ Ee (1)+ ⊤ X for every fixed  and .Therefore, if P{C ≥ } > 0 and Λ 0 () < ∞, then such s (0) (, t; , ) is uniformly bounded away from 0.
Before presenting the asymptotic results, we list the following regularity conditions.
(i) P{C ≥ } > 0 for the predetermined constant , and the baseline cumulative intensity function satisfies Λ 0 () < ∞. (ii) The covariate X is uniformly bounded, and there exists an M  >  such that the second moment of following random variables , where s (k) is the limit of S (k)  n in Lemma 2. The matrix Condition (i) is a mild requirement on the underlying intensity process, which is commonly adopted in recurrent data analysis.See Lin et al. (2000), Gao et al. (2019) and Huang et al. (2010).Condition (ii) is the sufficient condition for Lemma 2. Conditions (iii) are widely used in survival analysis, for example, see Andersen and Gill (1982), Lin et al. (2000) and Adekpedjou and Stocker (2015), among others.

Theorem 1. Under regularity conditions
Next, we investigate the asymptotic normality of the estimators.For notational convenience, define and By Theorem 1, Lemma 2, and a continuous mapping argument, it is readily seen that D n (⋅; γn , βn ) is a consistent estimator of D. Note that D n (t; , ) defined above is the derivative of Λn (t; , ) in ( 6) with respect to , and it will naturally appear in the estimator of the asymptotic variance of Λn (t; γn , βn ) when the delta method is applied.With the regularity conditions (i)-(iii), the weak convergence results of θn and Λn (t; γn , βn ) are stated below.
Theorem 2. Under regularity conditions (i)-(iii), we have (i) For the parametric part, √ n( θn −  0 ) is asymptotic normal with mean zero and variance which can be consistently estimated by where V n is given in Equation 8.
(ii) For the nonparameric part, { √ n( Λn (t; γn , βn ) − Λ 0 (t))} converges weakly to a mean-zero Gaussian process with covariance function and the covariance function can be consistently estimated by

Asymptotic efficiency of the estimators
We show that the estimators derived above are asymptotically efficient by deriving the efficient score functions for both the parametric and the nonparametric parts.For this purpose, we consider a generic point (, Λ) in the parameter space, where  = (,  ⊤ ) ⊤ .Consider the log-likelihood function (4) based on the observations from one system.The score function for the parameter  can be derived by taking derivative on the log-likelihood with respect to , which can be neatly written as To derive the tangent space for Λ, consider a parametric path   →  ,h through  with where h ∈  ∞ [0, ], the space of bounded functions over [0, ].For this parametric submodel, as we have already done in Section 2.2, the score function for the nonparametric part with respect to this submodel is where ) is a linear operator, and L 0 2 (P ,Λ ) is the space of mean-zero functions with finite second moments under P ,Λ .A similar result could be found in proposition 1 of Pena et al. (2001).
We begin with the efficiency of the parametric estimator θn .Notice that the efficient score for  is given by l1 (, Λ) = l1 (, Λ) − Π ,Λ l1 (, Λ), where Π ,Λ l1 (, Λ) is the orthogonal projection of l1 (, Λ) onto the closure of the linear space spanned by B ,Λ h for h ∈  ∞ [0, ].In the following, we show that B ,Λ h * is the projection of l1 (, Λ) onto the nuisance space, in which According to Lemma 2, h * ∈  ∞ [0, ] because s (1) (, x; , ) is uniformly bounded, and s (0) (, x; , ) is bounded from zero.Next, we show that h * satisfies the following orthogonal condition According to the score functions for the parameter part  and nonparametric part Λ in Equations 11 and ( 12), the left-hand side of the above display can be expressed as Use proposition 1 in Pena et al. (2001) to see that the first term in the expectation is and the second term is Notice that and Therefore, ( 14) is equal to With h * given in (13), it is readily seen that the above display equals 0 p+1 for any h ∈  ∞ [0, ].
That is, the function h * satisfies the orthogonal condition (14).and thus, the efficient score for  is and the efficient information is given by Compare ( 15) with the asymptotic covariance matrix of θn in Theorem 2 to see that the estimator θn is asymptotically efficient.
Next, we verify the asymptotic efficiency of the nonparametric estimator.Recall that the statistical model is  = {P ,Λ ∶  ∈ R p+1 , Λ ∈ } where  is the class of nondecreasing and right-continuous functions with Λ(0) = 0. Define a map  ∶   →  as (P ,Λ ) = Λ.According to theorem 25.48 of Van der Vaart (2000), verification of the asymptotic efficiency of Λ(⋅) at P ,Λ is equivalent to verifying that (i) the map  ∶   →  ∞ [0, ] is differentiable at P ,Λ in the sense of Van der Vaart (2000); and (ii) For each fixed t, the estimator Λ(t) is asymptotically efficient at P ,Λ for estimating Λ(t) for every t ∈ [0, ].
Next, we show that Λ(t) is asymptotically efficient at P ,Λ for estimating (P ,Λ )(t).Notice that Λ(t) can be considered as a parameter of interest for a fixed t.We use a similar argument when proving the efficiency of parametric part based on the differentiability of , which implies that (⋅)(t) is differentiable at P ,Λ .Following lemma 25.19 and relevant definitions in Van der Vaart (2000), it is obvious that the efficient influence function for Λ(t) is ψ,Λ (t).Thus, the efficient information for Λ(t) is given by Comparing the above display with the covariance function Σ Λ (t, t) stated in Theorem 2, we can conclude that the estimator Λn (t) for a fixed t with t ∈ [0, ] is asymptotically efficient.So far, we have verified all conditions in theorem 25.48 of Van der Vaart (2000).As a result, the nonparametric estimator Λn (⋅) is asymptotically efficient at P ,Λ .

SIMULATION
Simulations are conducted to assess the performance of the proposed semiparametric framework.Because the traditional method in Dauxois et al. (2019) is only applicable to a fixed number of failures, Section 4.1 first considers the failure censoring case for the comparison purpose.In Section 4.2, the time censoring setting is then considered and some modifications are made to the traditional method in order to make it applicable.At last, Section 4.3 considers the scenario with covariates.

Scenario without covariates (failure censoring)
In this subsection, we compare the proposed method with the traditional method in the case of failure censoring.The covariates are not considered here as the traditional method is not able to incorporate the covariate information.In specific, let T ij be the interoccurence time between the (j − 1)th and jth failures in the ith system, where i = 1, … , n and j = 1, … , n i .Define To avoid T inf > T sup , let In Dauxois et al. (2019), two estimators of  have been proposed and the one with better performance is considered here.In specific, the estimator of  is given by where N = max{n 1 , … , n n }, and ΛT j (t) is the empirical cumulative hazard function of the jth interoccurence time T j .The simulation settings are as follows.We set the baseline cumulative intensity function as Λ 0 (t) = t 2 , which corresponds to a Weibull distribution with shape parameter 2 and scale parameter 1.In addition, let the fixed number of the failures N = 5 and (j) = j − 1.To investigate the performance under a variety of scenarios, let  ∈ {0.1, 0.5, 1} and n ∈ {20, 50,100}.We first consider the accuracy in estimating the parameter .Based on 2000 replications, the estimated absolute biases and root mean square errors (rMSEs) under each setting are shown in Figure 1.As seen, the proposed ML estimator performs satisfactorily in all the scenarios and it outperforms the estimators by Dauxois et al. (2019).We then consider the estimation part.The  2, where the true curves are also plotted.Generally, the estimated cumulative intensity functions based on the proposed method match very well with the true curves in all the scenarios, while the estimated curves by the traditional method perform not so well when  = 0.1.

Scenario without covariates (time censoring)
In this subsection, we make a comparison between the proposed method and the traditional method under the time censoring mechanism, which is a more practical setting in real applications.In order to implement the traditional method, Equations 18-( 20) are again used.Because the data are time censored, this treatment essentially means the time between the last repair and the censoring time is discarded, and only the interoccurence times are utilized.
The censoring time C is assumed to follow the exponential distribution with mean 5 and the other settings are the same as those in Section 4.1.Figures 3 and 4, respectively, show the performance in estimating the parameter  and the baseline cumulative intensity function.As seen, the proposed method achieves high estimation accuracies in all the combinations of n and  and it outperforms the traditional method by a large margin.The poor performance of the traditional method is not surprising as its implementation requires a truncation of the original data and hence the information is not fully utilized.
To further validate our method, we then consider a different baseline cumulative intensity function Λ 0 (t) = √ t∕3, which corresponds to a Weibull distribution with shape parameter 0.5 and scale parameter 3, and the simulations results are shown in Figures 5 and 6.As expected, the proposed method is more accurate in estimating .More importantly, the nonparametric estimators based on the proposed method tally well with the true curves in all the scenarios, while the estimated curves by the traditional method deviate significantly from the true values.

Scenario with covariates
In this subsection, we investigate the performance of the proposed methods when the covariates are incorporated.We consider where ( 1 ,  2 ,  3 ) ⊤ = (0.5, −1, 0.1), and X 1 , X 2 , and X 3 , respectively, follow the standard normal distribution, the standard exponential distribution and the standard uniform distribution.All the other settings are the same as those in Section 4.1.
We first consider estimating the parameters  0 = (,  1 ,  2 ,  3 ) ⊤ .Figures 7 and 8, respectively, show the estimated absolute biases and rMSEs based on 2000 replications by considering n = 20, 50 and 100.As seen, all the parameters can be satisfactorily estimated and the estimation accuracy generally improves with the sample size n.On the other hand, the estimated baseline cumulative intensity function as well as the true values are plotted in Figure 9.It is observed that the nonparametric part can also be well estimated by the proposed method and when n ≥ 50, the estimated baseline cumulative intensity functions are almost identical to the true ones.

REAL DATA APPLICATION
In this section, we demonstrate the proposed models and methods by a real example from the automobile industry.The automobile manufacturer considered here is one of the largest automakers in the world and it has maintained a comprehensive warranty database which recorded detailed warranty claim data of different types of automobiles.The warranty claim data contain valuable information about product quality and reliability, and they can be used to, for example, give early warnings of faulty designs, predict future claims and warranty costs, and compare reliability with competing products.To achieve these, it is important to estimate the reliability of the components wisely.
As an illustrative example, the component window regulator is considered in this study.From the warranty claim database, we extract the repair records of the window regulator installed in 69 cars of the same model sold in a district, and the data are illustrated in Figure 10.Here, the red circle and the blue cross represent the date of sale and the date of each repair, respectively.To illustrate the methods, the end of study date is set as January 1, 2013, which is denoted by the red vertical line in Figure 10.Among the 69 cars, 54 cars have at least one failure record of the window regulator while the remaining 15 cars sold within the study period have no repair records.Based on the use rates and the sales dates, it is readily seen that all the cars are still within the two-dimensional warranty limit (i.e., 24 months or 60,000 miles warranty, whichever comes first) on the data-freeze date, and the censoring is caused by administrative censoring.In addition, it is observed from Figure 10 that the component failure time generally decreases with each subsequent repair, which indicates the effects brought by the imperfect repair.Thus, we use the EGFRR model to fit this dataset.In the dataset, the mileage information is available and it reflects the actual use rate of the car.Hence, it is likely that the component failure rate depends on the rate at which mileage accumulates on a car.In this study, the value of usage rate, which is defined as mileage per day, was calculated for each car with an observed mileage value prior to the end-of-study date.Let (j) = j − 1.Based on the estimation procedures in Section 2, we obtain that γ = 0.867 and β = 0.997, and the respective 95% confidence intervals are (0.391, 1.347) and (0.344, 1.650), based on the large-sample normal approximation.These results indicate that the imperfect repair and the use rate have adverse effects on the component failure rate.In addition, Figure 11 shows the estimated baseline cumulative hazard function and the corresponding baseline distribution function of the first interoccurence time, that is, the time between the sale date and the first repair, as well as the 95% confidence intervals based on the asymptotic results in Section 3. Note that the baseline cumulative hazard function of the first interoccurence time is essentially equivalent to the baseline cumulative intensity function Λ 0 (⋅).For comparison purpose, we have also applied the classical Cox proportional hazards model to fit the first interoccurence times.This is because there is no repair effect on the first failure times, and a satisfactory baseline cumulative hazard/distribution function estimator of the first interoccurence time should match well with that from the Cox model.In particular, the baseline cumulative hazard function based on the Cox model can be obtained by using the famous Breslow estimator and then the baseline distribution function can be readily obtained (Lin, 2007).As seen from Figure 11, the estimated baseline cumulative hazard function and distribution function of the first interoccurence time by the proposed method tally reasonably well with the curves by the Cox model, indicating an adequate fit of the EGFPP model.In addition, it is observed that the 95% pointwise confidence intervals by the proposed method are narrower than those based on the Cox model.This is reasonable as more failure data are utilized in the proposed method and hence more accurate estimation of the component lifetime distribution is expected.

CONCLUSIONS AND REMARKS
In this paper, we have proposed the semiparametric inference framework for the time-censored recurrent failure data with covariates under the EGFRR model.By careful projecting the score functions to the tangent spaces, we have shown that both the parametric and nonparametric estimators are asymptotically efficient.Extensive simulations and real examples have illustrated that the proposed framework outperforms the existing methods and it is useful to extract valuable reliability information from the recurrent failure data in presence of time censoring and covariates.
There are several research topics worth exploration.The first potential topic is about the exponent function (j).In this paper as well as many other existing studies, (j) is pre-given with the constraints that it is nondecreasing and (1) = 0.In practice, a more plausible way to determine (j) is to base on the recurrent failure data.One feasible solution is to first estimate (j) based on the empirical process and then verify the asymptotic properties of the estimator of (j).In addition, the  2 test may be used to test the goodness-of-fit of (j).Nevertheless, more thorough investigation is needed to look into the selection of (j).Another possible research direction is to consider the statistical inference of other popular imperfect repair models with time-censoring, such as the models in Finkelstein (2008), Wu and Zuo (2010) and Alaswad and Xiang (2017).If a similar martingale process can be constructed, the proposed method in this paper can be extended to those models and the asymptotic properties may be similarly established.

ACKNOWLEDGMENTS
We are grateful to the editors and two referees for their insightful comments that have led to a substantial improvement to an earlier version of the paper.Wang and Ye were supported by the National Science Foundation of China (72071138) and Singapore MOE AcRF Tier 2 grant (R-266-000-143-112).

Proof of Lemma 1
We use a similar proof as that in Pena et al. (2001).According to the definition of N(⋅, t), for each fixed t the compensator is given by A(s, t;  0 ,  0 , Λ

Proof of Lemma 2
We only show the proof for k = 0 here, since the proof for k = 1 and 2 is almost the same as that for k = 0. Consider the class of functions  = {g t, (X, Ñ, Ỹ ) = e  ⊤ X Y (, t; ) ∶ t ∈ [0, ],  ∈ } indexed by t and .It is easy to see that )2 , which is finite because X has bounded support, and E [ Y (0) (, 0; M  ) 2 ] < ∞.Further, the map t  → g t, (X, Ñ, Ỹ ), when treated as a function of t, is right-continuous and monotone decreasing.Thus, we can use dominated convergence theorem to see that Eg 2 t, is right continuous and monotone decreasing.Therefore, for a given  > 0, we first fix  ∈  and consider the partition 0 2, … , m are brackets of L 2 (P)-size , and the total number of brackets can be chosen of order (1∕) 2 .Next, we allow  to change.Notice that Let ε = ∕ √ E ̇g2 .For a fixed , we first find an ε-net  1 , … ,  q of  with respect to the Euclidean metric in R p+1 .For each j ≤ q, we follow the above procedure at the beginning of the proof to find brackets of L 2 (P)-size  for the functions {g t, j ∶ t ∈ [0, ]}.Enlarge each of the above brackets to [g (t k −), j − ε ̇g, g t k−1 , j + ε ̇g].All of these brackets have L 2 (P)-size 3.Repeat the above procedure for all j ≤ q.Here, q is not larger than the -covering number of , and is of order (1∕) p+1 .For each fixed j ≤ q, the -bracketing number for {g t, j ∶ t ∈ [0, ]} is of order (1∕) 2 .Therefore, the total number of -brackets under the L 2 (P) norm is of order (1∕) p+3 .This shows that the class of functions  is P-Donsker.

Proof of Theorem 1
To establish the consistency of parametric estimator θn , we let () = { ∶ || −  0 || ≤ ,  ∈ R p+1 }, which is a compact set in R p+1 .Notice that to prove the weak consistency is equivalent to verifying that for each  > 0, P( θn ∈ ()) converges to 1. Without loss of generality, we consider an arbitrary  such that () ⊂ .
First we are going to find a continuous function () such that sup ) where we denote  n (, ) by  n () for notational convenience.Define and let () be the expectation of l() under  0 ,  0 and Λ 0 .Then Since () is compact in R p+1 , we let K = () to derive (A3).Furthermore, the continuity of  can be easily derived from the continuity of  n and the uniformly convergence on every compact set in .
On the other hand, to show the consistency of the nonparametric part, we telescope using Λn (⋅;  0 ,  0 ) to get Because (γ n , β⊤ n ) ⊤ → P ( 0 ,  ⊤ 0 ) ⊤ and S (0) n (, x; , ) is a uniformly bounded continuous function with respect to , the second line above converges to 0 in probability.Combine the results to see that the nonparametric estimator Λn (⋅) is weakly consistent.For notational convenience, we denote this (1 + p + m)-dimensional vector by

Proof of theorem 2 Define
) .
Next, we use the Martingale CLT given in theorem 5. 3.5, Fleming and Harrington (2005) Absolute bias (aB) and root mean square error (rMSE) in estimating  based on the proposed method (new) and the traditional method (old) by considering n = 20, 50, and 100 under failure censoring [Colour figure can be viewed at wileyonlinelibrary.com]The true baseline cumulative intensity function (true) and the estimated ones by the proposed method (new) and the traditional method (old) by considering n = 20, 50, and 100 under failure censoring [Colour figure can be viewed at wileyonlinelibrary.com] mean values of 2000 estimated baseline cumulative intensity functions are shown in Figure Absolute bias (aB) and root mean square error (rMSE) in estimating  based on the proposed method (new) and the traditional method (old) by considering n = 20, 50, and 100 and Weibull shape 2 and scale 1 under time censoring [Colour figure can be viewed at wileyonlinelibrary.com]The true baseline cumulative intensity function (true) and the estimated ones by the proposed method (new) and the traditional method (old) by considering n = 20, 50, and 100 and Weibull shape 2 and scale 1 under time censoring [Colour figure can be viewed at wileyonlinelibrary.com] Absolute bias (aB) and root mean square error (rMSE) in estimating  based on the proposed method (new) and the traditional method (old) by considering n = 20, 50, and 100 and Weibull shape 0.5 and scale 3 under time censoring [Colour figure can be viewed at wileyonlinelibrary.com]The true baseline cumulative intensity function (true) and the estimated ones by the proposed method (new) and the traditional method (old) by considering n = 20, 50, and 100 and Weibull shape 0.5 and scale 3 under time censoring [Colour figure can be viewed at wileyonlinelibrary.com] Absolute bias in estimating  based on the proposed method by considering n = 20, 50, and 100 [Colour figure can be viewed at wileyonlinelibrary.com]Root mean square error in estimating  based on the proposed method by considering n = 20, 50, and 100 [Colour figure can be viewed at wileyonlinelibrary.com] The estimated baseline cumulative intensity function based on the proposed method by considering n = 20, 50, and 100 [Colour figure can be viewed at wileyonlinelibrary.com]Repair records of the window regulator in 69 cars.The red circle denotes the sale date, the blue cross denotes the repair date and the red vertical line denotes the censoring time [Colour figure can be viewed at wileyonlinelibrary.com] U R E 11 Estimated baseline cumulative hazard function and baseline distribution function of the first interoccurence time, as well as the 95% confidence intervals based on the proposed model and the Cox proportional hazards model [Colour figure can be viewed at wileyonlinelibrary.com] Piao Chen https://orcid.org/0000-0001-9714-450XR E F E R E N C E S How to cite this article: Wang, J., Chen, P., & Ye, Z. (2022).Efficient semiparametric estimation of time-censored intensity-reduction models for repairable systems.Scandinavian Journal of Statistics, 1-29.https://doi.org/10.1111/sjos.12564APPENDIX .