The Covariate Order Method for Nonparametric Censored Exponential Regression

Bo Henry Lindqvist  Jan Terje Kvaløy
bo@math.ntnu.no  jtk@math.ntnu.no

Department of Mathematical Sciences Norwegian University of Science and Technology N-7491 Trondheim

1. Introduction and Summary

In reliability engineering a first analysis of failure data is often based on the assumption of independent and identically exponentially distributed lifetimes. An often more realistic assumption is to assume that the lifetimes Z are exponentially distributed with hazard rate λ varying from unit to unit according to individual differences, for example due to varying environmental conditions. In many cases these differences can be represented by an observable covariate (vector) x, leading to a model where λ = λ(x) is a function of x. Estimating the hazard rate function λ(x) from observed data is called exponential regression. If some of the observations are right censored we have censored exponential regression. In this case there is a (possibly random) censoring time C such that only T = min(Z, C) and δ = I(Z ≤ C) are observed. Thus in addition to the observed time T we know whether there was a failure (δ = 1) or a censoring (δ = 0). Assume that we have observations from n independent units, where in addition to T, δ we observe the corresponding covariate vector X. This gives the data (T_{1}, δ_{1}, X_{1}), …, (T_{n}, δ_{n}, X_{n}).

As a simple example of the use in reliability studies, suppose that times Z between failures of some equipment are assumed to be exponentially distributed with hazard rate λ(x) depending on covariates such as age, duration of previous repair, type of repair etc. Censoring could here be used to represent for example planned maintenance.

Lindqvist et al (1988) considered lifetimes Z of subsea valves for oil/gas production. A typical covariate vector would then be x = (x_{1}, x_{2}, x_{3}), where x_{1} is the gas/oil ratio, x_{2} is the CO_{2}- content and x_{3} is the speed of oil through the valve. Censoring did here correspond to planned testing and maintenance operations (for example due to governmental regulations). A parametric model was used for the hazard rate λ(x), namely λ(x) = exp(β_{0} + β_{1}x_{1} + β_{2}x_{2} + β_{3}x_{3}). This is an example of a generalized linear model and the parameters can be estimated by ordinary maximum likelihood using standard computer programs. A seemingly more realistic model would be to replace the linear terms β_{i}x_{i} by general functions g_{i}(x_{i}) and then estimate the functions g_{i} by some nonparametric technique like spline smoothing, kernel smoothing or local polynomial smoothing. The purpose of this paper is to present a new nonparametric method for estimation in such exponential regression models. For the case of a one-dimensional covariate, the method is based on arranging the observations according to increasing order of the covariate, and using kernel estimation on a specific point process defined by the corresponding cumulative observation times. The method is then generalized to models with vector valued covariates.

2. Local Likelihood Methods

Before presenting the covariate order method we briefly review some well known competing methods. Let the model and data be as outlined in the previous section, but assume now that X is a one-dimensional covariate. Thus, for given X = x, Z is assumed to be exponentially distributed with an unknown hazard rate λ(x), that is f_{Z}(t|x) = λ(x) exp(−λ(x)t). Further, C is distributed according to some unknown censoring distribution f_{C}(t|x) which may depend on x, and C is assumed to be independent of Z given X.

Suppose first that λ(x) ≡ λ, i.e. the hazard rate does not depend on the covariate. Then the likelihood of the data is given by ∏_{i=1}^{n} λ^{δ_{i}} exp(−λT_{i}) so the log-likelihood is

∑_{i=1}^{n}(δ_{i} log λ − λT_{i})  (1)

This is maximized by λ̂ = ∑ δ_{i}/ ∑ T_{i}, which is the usual “occurrence/exposure ratio” estimator for λ.

\overline{}0  S_{1}  S_{2}  S_{k}  s ^{\frac{1}{n} T_{1}}^{\frac{1}{n} T_{2}}^{\frac{1}{n} T_{3}}  ^{\frac{1}{n} T_{n}} …

Figure 1: Construction of covariate order process.

If it is known that λ(x) varies slowly with x, then the above estimator suggests the following local estimator at a point x_{0}: λ̂(x_{0}) = ∑_{i∈N(x_{0})} δ_{i}/ ∑_{i∈N(x_{0})} T_{i}, where N(x_{0}) is a set of indices i for which x_{i} is sufficiently near x_{0} (where we have a choice regarding membership requirements for N(x_{0})). A formal way of producing such estimators is through a local likelihood approach, where we estimate λ(x_{0}) at a fixed x_{0} by maximizing

∑_{i=1}^{n} (δ_{i} log λ(x_{0}) − λ(x_{0})T_{i}) \frac{1}{h}K (\frac{x_{0} − x_{i}}{h} )  (2)

where K(·) a positive kernel function with integral 1, usually a unimodal density function, symmetric about 0 (e.g. the standard normal density). Note that the local likelihood (2) is similar to (1), but that the ith term is weighted according to the distance between x_{i} and x_{0}. Maximizing (2) leads to the estimator λ̂(x_{0}) = ∑ w_{i}δ_{i}/ ∑ w_{i}T_{i} where w_{i} = K((x_{0} − x_{i})/h). The value of h determines the smoothness of the estimator and is chosen by for example cross validation (see Section 3).

Finally we shall briefly mention local polynomial exponential regression, where in (2) we replace λ(x_{0}) by a local polynomial β_{0} + β_{1}(x_{i} − x_{0}) + β_{2}(x_{i} − x_{0})^{2} + ⋯, then estimate the parameters β_{i} by maximizing (2) and finally put λ̂(x_{0}) = β̂_{0}.

3. The Covariate Order Method

Let the situation be as in the previous section. The covariate order method consists in first arranging the observations (T_{1}, δ_{1}, X_{1}), …, (T_{n}, δ_{n}, X_{n}) such that X_{1} ≤ X_{2} ≤ ⋯ ≤ X_{n}. Next, for convenience, we divide the observation times by n, the total number of observations. Then let the scaled observation times T_{1}/n, …, T_{n}/n, irrespectively if they are censored or not, be subsequent intervals on an artificial time axis s. Let points which are endpoints of an interval corresponding to an uncensored observation (δ = 1) be considered as events, occurring at times denoted S_{1}, S_{2}, . . .. This is visualized in Figure 1, for an example where the ordered observations are (T_{1}, δ_{1} = 1), (T_{2}, δ_{2} = 0), (T_{3}, δ_{3} = 1), …, (T_{n−1}, δ_{n−1} = 0), (T_{n}, δ_{n} = 1).

Formally S_{i} = ∑_{j=1}^{k(i)} T_{j}/n where k(i) = min{r| ∑_{j=1}^{r} δ_{j} = i}. Note that if λ(x) ≡ λ is constant in x, then the process S_{1}, S_{2}, … is a homogeneous Poisson process. Indeed, given X_{1}, …, X_{n} we can regard the process on the s-axis as putting into operation sequentially n units, with i.i.d. exponential lifetimes, and where replacement takes place at either failure or censoring of the unit. Then at each point in time there is a constant hazard λ of experiencing a failure.

If λ(x) varies reasonably smoothly and not too much with x, then intuitively the process S_{1}, S_{2}, … could be regarded as nearly a nonhomogeneous Poisson process, for which the intensity function can be estimated by ordinary kernel density estimation based on the points S_{1}, S_{2}, . . .. An estimate of λ(x) then results by transforming this kernel estimate in an appropriate way to be described next.

Define N(s) for s > 0 to be the number of events S_{i} on the s-axis in the interval (0, s]. Then the conditional intensity of the process S_{1}, S_{2}, … at a point s is nλ(X_{N(s)+1}). This is because at time s the unit “in operation” is the one numbered N(s) + 1. The basic idea of the covariate order method is to estimate this intensity from the process S_{1}, S_{2}, … by for example a kernel estimator, as motivated above, yielding the estimator ρ̂(s). Then finally we invert the relation nλ̂(X_{N(s)+1}) = ρ̂_{n}(s) to obtain an estimate λ̂(x) at given points x. The key here is the relationship between the process X_{1}, X_{2}, … on the “covariate axis” and the process S_{1}, S_{2}, … on the s-axis. A possible way of estimating a general relationship between points on the two axes is by plotting points (X_{1}, T_{1}/n), …, (X_{n}, ∑_{i=1}^{n} T_{i}/n) and drawing straight lines between successive points, in other words using the function

s̃(x) = \frac{1}{n} ∑_{i=1}^{J−1} T_{i} + \frac{x − X_{J−1}}{X_{J} − X_{J−1}} \frac{1}{n}T_{J}, X_{J−1} < x ≤ X_{J}  (3)

We are now in position to write down an expression for the covariate order estimator of λ(x). Let K(·) be a positive kernel as defined earlier, and let h_{s} be a smoothing parameter which is either constant or varying along the s-axis. Assume that h_{s} → 0 as n → ∞ for all s. Further assume that there is a sequence h_{n} such that h_{s} ≥ h_{n} for all s, n where nh_{n} → ∞ as n → ∞. Then the estimator

λ̂(x) = ρ̂(s̃(x)) ≡ \frac{1}{nh_{s}} ∑_{i=1}^{k} K (\frac{s̃(x) − S_{i}}{h_{s}} )  (4)

can be shown to be a uniformly consistent estimator of λ(x). Here k = ∑_{i=1}^{n} δ_{i} is the number of uncensored observations.

In practice the estimated correspondence function (3) may be replaced by more sophisticated estimators, improving on the smoothness of the estimator (4).

To avoid the estimator λ̂(x) to be seriously downward biased near the endpoints special care must be taken at the boundaries. Viewed only as a problem on the s-axis the estimator (4) is simply density estimation on the s-axis, and techniques for handling boundary problems in density estimation can be adopted. A common technique is to reflect the data points around both endpoints, see for example Silverman (1986).

In order to choose the “best” value of the smoothing parameter h_{s}, a likelihood cross-validation criterion has been used, see for example Simonoff (1996).

4. Several Covariates

The mentioned methods for nonparametric exponential regression can, at least in principle, be directly generalized to multidimensional covariates with, say, m components. For the local methods the generalization is in principle immediate, by using m-dimensional kernels. For the covariate order method it is less clear, however, how to generalize to multidimensional x.

In any case direct generalizations will soon be beaten by the so called “curse of dimensionality”. A better approach in many cases is to reduce the dimension of the problem by assuming some structure on the covariate space. Such an approach is also more appropriate for generalizing the covariate order method.

One way to proceed is to assume that the hazard rate can be written in the form of the generalized additive model

λ(x) = exp(α + g_{1}(x_{1})+, …, +g_{m}(x_{m})),  (5)

where x = (x_{1}, …, x_{m}), and where g_{1}(·), …, g_{m}(·) are unspecified smooth functions. These functions can be estimated iteratively by using the one dimensional covariate order method repeatedly. The algorithm below is similar to the backfitting algorithm of Hastie and Tibshirani (1990).

Let ĝ_{i}^{(d)}(x_{i}) be the estimate of g_{i}(x_{i}) after the dth iteration of the algorithm, and define ĝ_{i}^{(0)}(x_{i}) = 0. Then the algorithm can be written as follows:

1. Repeat 2-3 for d = 1, …, D where |α̂^{(D−1)} − α̂^{(D)}|/α̂^{(D−1)} < ϵ.

2. Calculate α̂^{(d)} = log (∑_{i=1}^{n} δ_{i}/ ∑_{i=1}^{n} T_{i} exp(ĝ_{1}^{(d−1)}(X_{i1}) + … + ĝ_{m}^{(d−1)}(X_{im}))).

3. Repeat 4-5 for j = 1, …, m.

4. Transform T_{1}, …, T_{n} to L_{1j}, …, L_{nj} where L_{ij} = T_{i} exp(α̂^{(d)} +ĝ_{1}^{(d)}(X_{i1})+⋯+ĝ_{j−1}^{(d)}(X_{i,j−1})+ĝ_{j+1}^{(d−1)}(X_{i,j+1})+⋯+ĝ_{m}^{(d−1)}(X_{im})), i = 1, …, n.

5. Calculate the estimate g̃_{j}^{(d)}(x_{j}) from (L_{1j}, δ_{1}, X_{1j}), …, (L_{nj}, δ_{n}, X_{nj}) using the covariate order method, and let ĝ_{j}^{(d)}(x_{j}) = g̃_{j}^{(d)}(x_{j}) − ∑_{i=1}^{n} g̃_{j}^{(d)}(X_{ij})/n.

The motivation for the algorithm is the observation that if Z is exponentially distributed with parameter exp(α+g_{1}(x_{1})+…+g_{m}(x_{m})), then Z exp(α+g_{1}(x_{1})+…+g_{j−1}(x_{j−1})+g_{j+1}(x_{j+1})+…+g_{m}(x_{m})) is exponentially distributed with parameter exp(g_{j}(x_{j})). This explains step 4 and 5 in the algorithm. In step 5 all the estimated functions are normalized to have mean 0 to avoid possible problems caused by constant factors jumping back and forth between the estimated functions. Finally the estimate α̂^{(d)} in step 2 is the maximum likelihood estimate based on the sample (L_{1}, δ_{1}), …, (L_{n}, δ_{n}) where L_{i} = T_{i} exp(ĝ_{1}^{(d−1)}(X_{i1}) + … + ĝ_{m}^{(d−1)}(X_{im}))

Notice that the algorithm can easily be extended to allow some of the g-functions to be parametric (for example for discrete covariates) and then estimate the parameters of these functions in step 2 of the algorithm using maximum likelihood.

5. Nonparametric Tests for Covariate Effect

It was noted in Section 3 that if λ(x) ≡ λ, i.e. there is no covariate effect, then the process S_{1}, S_{2}, … is a homogeneous Poisson process (HPP). This observation suggests that statistical tests for the null hypothesis of an HPP versus various non-HPP alternatives can be applied as tests for covariate effect in this situation. The further presumption that the process S_{1}, S_{2}, … can be considered as being approximately an NHPP when λ(x) depends on x, suggests adopting exisitng tests for HPP versus NHPP, see for instance Ascher and Feingold (1984) and references therein. A detailed presentation of this approach for testing covariate effects is given in Kvaløy (1999b).

6. Concluding Remarks

A comprehensive simulation study (reported in Kvaløy and Lindqvist, 1998) shows that the covariate order method is generally competitive when compared to the local methods reviewed in Section 2. In particular the covariate order method seems to be doing better than the other methods in “difficult” cases with few observations or large degree of censoring. Also, the covariate order method has a large flexibility due to the possibility of utilizing any method for density estimation.

It is finally noted that the covariate order method can be adapted for use in the Cox-regression setting. The idea here is that Cox-regression can be reduced to exponential regression by transforming the lifetimes by the estimated cumulative baseline hazard function. Thus our method is also a method for nonparametric estimation in Cox-models. This has been further explored by Kvaløy (1999a).

References

1. Ascher, H. and Feingold, H. (1984). Repairable Systems Reliability. New York: Marcel Dekker.

2. Hastie, T. J. and Tibshirani, R. J. (1990). Generalized Additive Models. London: Chapman & Hall.

3. Kvaløy, J. T. (1999a). Nonparametric estimation in Cox-models: Time transformation versus partial likelihood methods. Statistics No. 3/99, Department of Mathematical Sciences, Norwegian University of Science and Technology.

4. Kvaløy, J. T. (1999b). Covariate order tests for covariate effect. Statistics No. 7/99, Department of Mathematical Sciences, Norwegian University of Science and Technology.

5. Kvaløy, J. T. and Lindqvist, B. H. (1998). The covariate order method for censored exponential regression. Statistics No. 10/98, Department of Mathematical Sciences, Norwegian University of Science and Technology.

6. Lindqvist, B., Molnes, E. and Rausand, M. (1988). Analysis of SCSSV performance data. Reliability Engineering and System Safety, 20, 3-17.

7. Silverman, B. W. (1986). Density Estimation. London: Chapman & Hall.

8. Simonoff, J. S. (1996). Smoothing Methods in Statistics. New York: Springer Verlag.