# Parametric Regression Models Using Reversed Hazard Rates

Authors: Asokan Mulayath Variyath; P. G. Sankaran
Journal: Journal of Probability and Statistics
 (2014)

---

## Abstract

Proportional hazard regression models are widely used in survival analysis to understand and exploit the relationship between survival time and covariates. For left censored survival times, reversed hazard rate functions are more appropriate. In this paper, we develop a parametric proportional hazard rates model using an inverted Weibull distribution. The estimation and construction of confidence intervals for the parameters are discussed. We assess the performance of the proposed procedure based on a large number of Monte Carlo simulations. We illustrate the proposed method using a real case example.


---

## Body

## 1. Introduction

In survival studies, covariates or explanatory variables are usually employed to represent heterogeneity in a population. The main objective in such situations is to understand and exploit the relationship between lifetime and covariates. Regression models are useful in such contexts to assess the effect of covariates on lifetime. These models can be formulated in many ways and several types are in common use. Parametric regression models for lifetime involve specification for the distribution of a lifetimeT given a vector of covariates X. The most commonly used parametric model is the Weibull regression model, which satisfies the proportional relationship between hazard rate functions of the lifetimes of two subjects. The maximum likelihood technique is usually employed to find estimates of the parameters of the model. For more properties and applications of parametric regression models, one should refer to Lawless [1].In survival studies, there are many occasions where lifetime data are left censored. For example, baboons in the Amboseli Reserve, Kenya, sleep in the trees and descend for ageing at certain times of the day. Observers often arrive later in the day after this descent has occurred and on such days they can only ascertain that the descent took place before a particular time, so that the descent times are left censored (see [2]). On such occasions, a reversed hazard rate is more appropriate than a hazard rate to analyze lifetime data due to the fact that estimators of hazard rates are unstable when data are left censored. The reversed hazard rate of T is defined as
(1)λ(t)=limΔt→0P(t-Δt≤T∣T≤t)Δt.Introduced by Barlow et al. [3], the function λ(t) has been used in various contexts such as the estimation of distribution function under left censoring [1], defining a new stochastic order [4], characterization of lifetime distributions [5–7], studying ageing behavior [8, 9], evolving new repair and maintenance strategies [10, 11], the mixed proportional hazards model [12], and stress hybrid hazards model [13].Recently, Sengupta and Nanda [14] introduced the proportional reversed hazards model in a semiparametric setup. In the present work, we introduce a fully parametric regression model that satisfies the proportional reversed hazards property. The inverted Weibull distribution is employed as a lifetime model, which can be extended to any parametric model. A large number of simulation studies indicate that the proposed approach is performing well.The rest of the paper is organized as follows. In Section2, we introduce a parametric regression model using an inverted Weibull distribution. The proposed model has the property that the reversed hazard rate for the lifetime of pair of subjects is proportional. The estimation of the parameters of the model is discussed in Section 3. Simulation studies are conducted in Section 4 to assess the finite sample behavior of the estimators. The proposed model is applied to real life data in Section 5 to illustrate its utility. Finally Section 6 provides the major conclusions of the study.


## 2. Statistical Model

LetT be a nonnegative random variable representing the lifetime of a subject with the distribution function F(t). Assume that the probability density function of T, f(t), exists. The reversed hazard rate of T given in (1) can be written as
(2)λ(t)=f(t)F(t).
Let X be p×1 vector of auxiliary information which may be time dependent. The proportional reversed hazard (PRH) model is defined by
(3)λ(t∣X)=λ0(t)g(β;X),
where λ0(t) is the baseline reversed hazard rate, g(·) is a nonnegative function of X and β,p×1 vector of regression parameters, and λ(t∣x) is the reversed hazard rate of T given the covariate X. The PRH model can be expressed in terms of the distribution function as
(4)F(t∣X)=(F0(t))g(β;X),
where F(t∣X) is the distribution function of T given X and F0(t) is the baseline distribution function in the absence of covariates. It should be noted that, for two subjects, the ratio of reversed hazard rates is independent of the time t. Semiparametric analysis of the model (2) is recently discussed in Sengupta and Nanda [14]. Our objective here is to carry out the parametric analysis of an inverted Weibull distribution under left censoring. When the lifetime random variable follows the inverted Weibull distribution, the baseline distribution function is given by
(5)F0(t)=e-γ/tα,t>0;α,γ>0.
The baseline reversed hazard rate of T is then obtained as(6)λ0(t)=γαtα+1.
Note that the baseline reversed hazard rate is decreasing as t increases. In the presence of the covariate X, we have
(7)λ(t∣X)=γαtα+1g(β;X).
We assume that g(β;X)=eβ′X so that
(8)F(t∣X)=e(-γ/tα)eβ′x,
with
(9)f(t∣X)=γαeβ′xtα+1e-(γ/tα)eβ′x.
Suppose that the lifetime random variable T is randomly left censored by Z. In practice, one could observe the vectors (Y,δ,X), where Y=max(T,Z) and δ=I(T=Y) with I(·) being the indicator function. Let (yi,δi,xi), i=1,2,…,n, be i.i.d. copies of (Y,δ,X). Then the likelihood function can be written as
(10)L(β;y)=∏i=1n‍f(yi∣xi)δiF(yi∣xi)1-δi.
Under the inverted Weibull distribution assumption, the likelihood function given in (10) is obtained as
(11)L(β,α,γ;y)∝∏i=1n‍(γαeβ′xiyiα+1e-(γ/yiα)eβ′xi)δi(e(-γ/yiα)eβ′xi)1-δi,
so that the log likelihood function is
(12)l(β,α,γ;y)=logc+∑i=1n‍δiβ′xi+∑i=1n‍δi(logγ+logα)-(α+1)∑i=1n‍δilogyi-γ∑i=1n‍eβ′xiyiα,
where c is a real constant independent of β, γ, and α. We maximize (12) to estimate the parameters β, α, and γ by equating the partial derivatives with respect to each parameter to zero as
(13)∂l∂β=∑i=1n‍δixi-γ∑i=1n‍eβ′xixiyiα=0,∂l∂α=∑i=1n‍δiα-∑i=1n‍δilogyi+γ∑i=1n‍eβ′xilogyiyiα=0,∂l∂γ=∑i=1n‍δiγ-∑i=1n‍eβ′xiyiα=0.
Since there is no closed form solution available for (13), we use numerical methods to estimate the parameters. The observed information matrix is given by
(14)I(β,α,γ)=[-∂2l∂β∂β-∂2l∂β∂α-∂2l∂β∂γ-∂2l∂α∂β-∂2l∂α2-∂2l∂α∂γ-∂2l∂γ∂β-∂2l∂γ∂α-∂2l∂γ2].
Note that the matrix (14) is of order (p+2)×(p+2). Under the standard regularity conditions, the vector of estimates (β^,α^,γ^) is asymptotically (p+2)-variate normal with mean vector (β,α,γ) and dispersion matrix I*-1, where I* is the Fisher information matrix obtained from I by taking the expected values of each entry.There are different algorithms available to estimate the parameters by solving the score equations or directly optimizing the likelihood function. The Newton-Raphson method is the most common method used to estimate since it is easy to determine the derivatives of the score equations. In this numerical iterative method, the initial values play a vital role due to the logarithm function. In the simulation studies given in Section4, we use the simplex method proposed by Neldar and Mead [15] to estimate the parameters. The simplex method is a simple method to use to estimate the parameters by maximizing the likelihood function where we do not need the derivatives of the function to be optimized.


## 3. Testing and Confidence Intervals forβ

Tests and interval estimates of parameters can be derived by the likelihood ratio test procedure. We are mainly interested in the regression parameterβ where the parameters θ=(α,γ) are normally considered as nuisance parameters.Let thep vector regression parameter be denoted as β=(β1,β2), where β1 and β2 are of vectors of sizes k and p-k, respectively, and θ is the other parameter in the model. We are interested in testing
(15)H0:β1=β10againstH1:β1≠β10,
where β10 is the specified regression parameter value. To test H0, we construct the likelihood ratio statistic
(16)Λ*=2l(β^1,β^2,θ^)-2l(β10,β^2,θ^),
where β^1, β^2, and θ^ are the maximum likelihood estimates under the full model. For a large value of n, Λ* follows the χk2 distribution under the null hypothesis.Alternatively, we can use the test statistic(17)Λ1=(β^1-β10)V11-1(β^1-β10),
where V11 can be obtained from V=I(β^,θ^)-1, which is partitioned as
(18)V=[V11V12V21V22].
Under the null hypothesis, Λ1 follows χk2 distribution.Assuming asymptotic normality, we can construct the100(1-α)% confidence interval for the individual regression parameter βj as
(19)β^j±zα/2s.e.(β^j),
where s.e.(β^j) can be obtained from V11.Another important problem is the selection of important covariates in the proportional reversed hazard models. Since we assume a parametric model, we can use variable selection methods such as the Akaike information criterion (AIC) and the Bayesian information criterion (BIC). To test adequacy of the parametric model, Cox-Snell residuals can be used, which is explained in Section5 with a case example.


## 4. Performance Analysis

To assess the performance of the proposed method, we carried out a large number of simulations. We generated samples of size 100 from an inverted Weibull distribution, with different values of parameters asα=(0.5,1,1.5,2) and γ=(0.5,1,1.5,2). We considered a single covariate which is generated from Uniform (0,1) and a regression parameter assumed to be β=(0.5,1.1.5,2). We developed the censoring mechanism using the random data generated from the exponential distribution with the parameter λ. We choose the value of λ such that the percent of censored data is between 10 and 20 percent. We used the simplex method proposed by Neldar and Mead [15] to estimate the parameters. We repeated the study for 10000 times and computed the mean and standard deviation of the parameter estimates. The entire study was repeated for a sample size of 250. The summary of the parameter estimates is given in Table 1.Table 1
Averages and standard deviation (in brackets) of parameter estimates.
True parameter values
Estimate of parameters (n=100)
Estimate of parameters (n=250)
α
β
γ
α
^
β
^
γ
^
α
^
β
^
γ
^
1
1
1
1.0228
1.0189
1.0351
1.0080
1.0127
1.0096
(0.0924)
(0.4040)
(0.2486)
(0.0556)
(0.2434)
(0.1444)
1
0.5
1
1.0199
0.5104
1.0287
1.0075
0.5046
1.0115
(0.0906)
(0.3984)
(0.2389)
(0.0564)
(0.2489)
(0.1476)
1
1.5
1
1.0196
1.5321
1.0302
1.0083
1.5113
1.0121
(0.0877)
(0.4033)
(0.2413)
(0.0547)
(0.2457)
(0.1430)
1
2
1
1.0200
2.0348
1.0313
1.0090
2.0178
1.0113
(0.0862)
(0.4123)
(0.2384)
(0.0535)
(0.2530)
(0.1423)
1
1
0.5
1.0234
1.0260
0.5107
1.0085
1.0090
0.5049
(0.0977)
(0.4338)
(0.1347)
(0.0607)
(0.2708)
(0.0824)
1
1
1.5
1.0199
1.0191
1.5590
1.0081
1.0064
1.5225
(0.0860)
(0.3848)
(0.3553)
(0.0539)
(0.2381)
(0.2128)
1
1
2
1.0200
1.0223
2.0800
1.0081
1.0104
2.0301
(0.0843)
(0.3759)
(0.4708)
(0.0534)
(0.2361)
(0.2823)
0.5
1
1
0.5099
1.0277
1.0252
0.5040
1.0088
1.0109
(0.0459)
(0.3986)
(0.2400)
(0.0285)
(0.2452)
(0.1459)
1.5
1
1
1.5308
1.0156
1.0372
1.5142
1.0081
1.0121
(0.1330)
(0.4072)
(0.2480)
(0.0814)
(0.2474)
(0.1449)
2
1
1
2.0460
1.0229
1.0329
2.0160
1.0129
1.0100
(0.1794)
(0.4050)
(0.2444)
(0.1111)
(0.2512)
(0.1484)From Table1 we can see that the mean of the parameter estimates based on 10000 simulation is very close to the true parameter values and the standard deviation is also small. When the sample size increases the standard error of the estimates decreases and bias reduces. It should be noted that there is a slight positive bias in all cases, even though it is negligibly small. Since there are no comparable models based on reversed hazard rates, we did not perform any comparison studies.


## 5. An Example

We consider an extract of left censored data from an Australian twin study given in Duffy et al. [16]. The data consist of information on the age of appendectomy of monozygotic and dizygotic twins. There are observations with missing age at onset and therefore the data are left censored. The individuals having age at onset of less than 11 are left censored. The covariate, namely, Zygocity, has values from 1 to 6. This data set consists of 54 observations of which 15 are left censored. We use this data to illustrate the utility of the parametric reverse hazard rate model. Probability plotting and statistical test confirmed the distribution of data as inverted Weibull distribution. We use the simplex method to estimate the parameters. Since the parameter values are unknown and to avoid the effect of inappropriate initial values, we consider different initial values and choose the estimates which have maximum likelihood. Estimates of the parameters are α^=2.3940, β^=-0.0142, and γ^=444.0586. The 95% confidence interval for β indicates that the regression coefficient corresponding to Zygocity is not significantly different than zero; that is, the effect of Zygocity is negligible. This conclusion is also verified through the likelihood ratio test statistic value 0.0074 having a P value of 0.93.We use a Cox-Snell residual plot to assess the goodness of fit. The Cox-Snell residual is defined by(20)ri=-logF(yi)=γyiαeβxi,i=1,…,n.
If the model fits the data, then the residuals should have a standard exponential distribution, so that a hazard plot of residuals versus the Nelson-Aalen estimator of the cumulative hazard of the residuals will be a straight line with slope one. A plot of Cox-Snell residuals against the Nelson-Aalen estimates of the cumulative hazard rate of residuals is given in Figure 1, which shows that the fit is reasonably good.Figure 1
Plot of cumulative hazard rates of Cox-Snell residuals versus residuals.


## 6. Conclusions

Proportional reversed hazard rate models are more suitable for modeling the left censored lifetime data. In this paper, we proposed a parametric PRH model assuming that the lifetime data follows an inverted Weibull distribution. The estimation and hypothesis testing of the parameters of the model have been discussed in detail. The performance of the proposed model is assessed based on a large number of Monte Carlo simulations. Our simulations results clearly indicated that the proposed model is performing well. We applied the proposed model to a real life example to illustrate the utility of the method. Recently, Bayesian methodology was extensively employed in the analysis of lifetime data. The inference procedures of the proposed model by selecting an appropriate prior distribution are topics of research to be explored. The present work can be easily extended to any location-scale families of distributions.



---

*Source: 645719-2014-01-06.xml*