Chapter 15: Density and Hazard Estimation

HAZEZ

Performs nonparametric hazard rate estimation using kernel functions. Easy-to-use version of HAZRD.

Required Arguments

XNOBS by m matrix containing the raw data, where m = 1 if ICEN = 0, and m = 2 otherwise.   (Input)

IRT — Column number in X containing the times of occurrence of the events.   (Input)

ST — Vector of length NOBS containing the times of occurrence of the events, sorted from smallest to largest.   (Output)
Vector ST is obtained from matrix X and is used as input to routine HAZST.

JCEN — Vector of length NOBS containing the sorted censor codes.   (Output)
Censor codes are sorted corresponding to the events ST(i), with censored observations preceding tied failures. Vector JCEN is obtained from the censor codes in X and is used as input to routine HAZST.

ALPHA — Optimal estimate for the parameter α.   (Output)

BTA — Optimal estimate for the parameter β.   (Output)

K — Optimal estimate for the parameter k.   (Output)

VML — Optimal value of the criterion function.   (Output)
VML is the “modified likelihood”.

H — Vector of length 5 * NOBS containing the constants needed to compute the k-th nearest failure distance and the observation weights.   (Output)
H is used as input to routine HAZST.

Optional Arguments

NOBS — Number of observations.   (Input)
Default: NOBS = size (X,1).

LDX — Leading dimension of X exactly as specified in the dimension statement in the calling program.   (Input)
Default: LDX = size (X,1).

ICEN — Censoring option.   (Input)
If ICEN = 0, then all of the data is treated as exact data with no censoring. For
ICEN > 0, column ICEN of X contains the censoring codes. A censoring code of 0 means an exact event (failure). A censoring code of 1 means that the observation was right censored at the event time.
Default: ICEN = 0.

IPRINT — Printing option.   (Input)
If IPRINT = 1, the grid estimates and the optimized estimates are printed for each value of k. Otherwise, no printing is performed.
Default : IPRINT = 0.

NMISS — Number of missing (NaN, not a number) values in X.   (Output)

FORTRAN 90 Interface

Generic:          CALL HAZEZ (X, IRT, ST, JCEN, ALPHA, BTA, K, VML, H [,…])

Specific:                             The specific interface names are S_HAZEZ and D_HAZEZ.

FORTRAN 77 Interface

Single:           CALL HAZEZ (NOBS, X, LDX, IRT, ICEN, IPRINT, ST, JCEN, ALPHA, BTA, K, VML, H, NMISS)

Double:                              The double precision name is DHAZEZ.

Description

Routine HAZEZ is an implementation of the methods discussed by Tanner and Wong (1984) for estimating the hazard rate in survival or reliability data with right censoring. It uses the biweight kernel,

and a modified likelihood to obtain data-based estimates of the smoothing parameters α, β, and k needed in the estimation of the hazard rate. For kernel K(x), define the “smoothed” kernel
Ks(xx(j)) as follows:

where djk is the distance to the k-th nearest failure from x(j), and x(j)is the j-th ordered observation (from smallest to largest). For given α and β, the hazard at point x is given by:

where N = NOBS, δi is the censor code (0 = failed, 1 = censored) for the i-th ordered observation, and wi is the weight of the i-th ordered observation (given by 1/(Ni + 1)). The hazard may be estimated via routine  after the smoothing parameters have been obtained

Let

The likelihood is given by:

where Π denotes product. Since the likelihood, as specified, will lead to degenerate estimates, Tanner and Wong (1984) suggest the use of a modified likelihood. The modification consists of deleting the observation xi in the calculation of h(xi) and H(xi) when the likelihood term for xi is computed. For a given k, α and β can then be estimated via the usual optimization techniques.

Estimates for α and β are computed as follows. For a given β, a closed form solution is available for α. The problem is thus reduced to the estimation of β. To estimate α and β, a grid search is first performed. Experience indicates that if the initial estimate of β from this grid search is greater than exp(6), then the modified likelihood is degenerate because the hazard rate does not change with time. In this situation, β should be taken to be infinite, and an estimate of α corresponding to infinite β is computed directly. When the estimate of β from the grid search is less than exp(6) (approximately 400), a secant algorithm is used to optimize the modified likelihood. The secant algorithm is said to have converged when the change in β from one iteration to the next is less than 0.00001. Additionally, convergence is assumed when the value of β becomes greater than exp(6). This corresponds to an infinite β with a degenerate likelihood.

A grid of k-values is used to find the optimum value of the likelihood with respect to k. The grid is determined by HAZEZ and consists of at most 10 points. The starting value in the grid is the smallest possible value of k. An increment of 2 is then used to obtain the remaining grid points.

For each grid value, the modified likelihood is optimized with respect to α and β. That grid point, which leads to the smallest likelihood, is taken to be the optimal k.

Comments

1.         Informational errors

Type Code

4         6                  All observations are missing (NaN, not a number) values.

4         7                  There are not enough failing observations in X to continue.

2.         The grid values in the initial grid search are given as follows: Let
β* = − 8, − 4, − 2, − 1, − 0.5,0.5,1, and 2, and

            For each value of β, VML is computed at the optimizing β. The maximizing β is used to initiate the iterations.

3.         If the initial β* is determined from the grid search to be less than −6, then it is presumed that β is infinite, and an analytic estimate of α based upon infinite β is used. Infinite β corresponds to a flat hazard rate.

Programming Notes

1.         Routine HAZST may be used to estimate the hazard on a grid of points once the optimal values for α, β and k have been found. (The user should also consider using routine HAZRD, which allows for more options than HAZEZ.)

2.         Routine HAZEZ assumes that censored observations precede failed observations at tied failure/censoring times.

3.         The secant algorithm iterates on the transformed parameter β* = exp(−β). This assures a positive β, and it also seems to lead to a more desirable grid search. All results returned to the user are in the original parameterization.

Example

The following example is illustrated in Tanner and Wong (1984), and the data are taken from Stablein, Carter, and Novak (1981). It involves the survival times of individuals with nonresectable gastric carcinoma. Only those individuals treated with radiation and chemotherapy are used.

 

      USE HAZEZ_INT

      USE UMACH_INT

      USE WRRRN_INT

      USE WRIRN_INT

 

      IMPLICIT   NONE

      INTEGER    ICEN, IPRINT, IRT, LDX, NOBS

      PARAMETER  (ICEN=2, IPRINT=1, IRT=1, LDX=45, NOBS=45)

!

      INTEGER    JCEN(NOBS), K, NMISS, NOUT

      REAL       ALPHA, BTA, H(5*NOBS), ST(NOBS), VML, X(NOBS,2)

!

      DATA X/17, 42, 44, 48, 60, 72, 74, 95, 103, 108, 122, 144, 167, &

          170, 183, 185, 193, 195, 197, 208, 234, 235, 254, 307, 315, &

          401, 445, 464, 484, 528, 542, 567, 577, 580, 795, 855, 882, &

          892, 1031, 1033, 1306, 1335, 1366, 1452, 1472, 36*0, 9*1/

!

      CALL HAZEZ (X, IRT, ST, JCEN, ALPHA, BTA, K, VML, H, ICEN=ICEN, &

                 IPRINT=IPRINT, NMISS=NMISS)

!

      CALL UMACH (2, NOUT)

      WRITE (NOUT,99999) NMISS

99999 FORMAT (/' NMISS = ', I4/)

      CALL WRRRN ('ST', ST, 1, NOBS, 1)

      CALL WRIRN ('JCEN', JCEN, 1, NOBS, 1)

      CALL WRRRN ('H', H, NOBS, 5, NOBS)

      END

Output

 

                  *** GRID SEARCH FOR K =     2 ***

           ALPHA                   BETA                   VML

          65.157967            2980.958008            -266.543945

          32.434208              54.598148            -262.551147

          17.100269              20.085537            -263.100769

          11.402525              12.182494            -264.410187

           7.263529               7.389056            -267.502014

           4.452315               4.481689            -270.548523

           2.689497               2.718282            -335.407288

           1.633702               1.648721            -338.566162

           0.995799               1.000000            -519.939514


                  *** OPTIMAL PARAMETER ESTIMATES ***

           ALPHA                   BETA                   VML

          32.219337              53.968315            -262.550781


                  *** GRID SEARCH FOR K =     4 ***

           ALPHA                   BETA                   VML

          25.596716            2980.958008            -266.471558

          20.476425              54.598148            -262.893860

          13.995192              20.085537            -262.792755

          10.109113              12.182494            -262.573212

           6.883837               7.389056            -263.030121

           4.407142               4.481689            -265.238647

           2.690131               2.718282            -265.606293

           1.633339               1.648721            -266.822693

           0.993371               1.000000            -271.831390


                  *** OPTIMAL PARAMETER ESTIMATES ***

           ALPHA                   BETA                   VML

           8.530729               9.683726            -262.545593


                  *** GRID SEARCH FOR K =     6 ***

           ALPHA                   BETA                   VML

          16.828691            2980.958008            -266.729248

          14.840095              54.598148            -264.019409

          11.215133              20.085537            -262.844360

           9.013870              12.182494            -263.663391

           6.557755               7.389056            -263.283752

           4.330785               4.481689            -263.732697

           2.691744               2.718282            -264.613800

           1.633932               1.648721            -265.381866

           0.990891               1.000000            -266.242767


                  *** OPTIMAL PARAMETER ESTIMATES ***

           ALPHA                   BETA                   VML

          12.553377              28.178671            -262.529877


                  *** GRID SEARCH FOR K =     8 ***

           ALPHA                   BETA                   VML

          11.377748            2980.958008            -266.746185

          10.773529              54.598148            -265.469299

           8.766835              20.085537            -262.476807

           7.427887              12.182494            -263.109009

           5.916299               7.389056            -264.492432

           4.184323               4.481689            -264.289886

           2.656351               2.718282            -264.807617

           1.623750               1.648721            -265.270691

           0.989442               1.000000            -264.738403


                  *** OPTIMAL PARAMETER ESTIMATES ***

           ALPHA                   BETA                   VML

           8.522110              18.281288            -262.438568


                  *** GRID SEARCH FOR K =    10 ***

           ALPHA                   BETA                   VML

           8.689023            2980.958008            -267.026093

           8.412854              54.598148            -266.250366

           7.196551              20.085537            -263.192688

           6.207793              12.182494            -262.648376

           5.143391               7.389056            -264.274384

           3.934601               4.481689            -264.523193

           2.630993               2.718282            -264.877869

           1.611710               1.648721            -264.332581

           0.984530               1.000000            -263.920013


                  *** OPTIMAL PARAMETER ESTIMATES ***

           ALPHA                   BETA                   VML

           6.483376              13.956067            -262.589661


                  *** GRID SEARCH FOR K =    12 ***

           ALPHA                   BETA                   VML

           6.669007            2980.958008            -266.778259

           6.551508              54.598148            -266.347595

           5.933167              20.085537            -264.141174

           5.252526              12.182494            -262.516205

           4.471936               7.389056            -262.691589

           3.598284               4.481689            -263.914032

           2.557817               2.718282            -263.390106

           1.588307               1.648721            -263.879578

           0.973723               1.000000            -263.361908


                  *** OPTIMAL PARAMETER ESTIMATES ***

           ALPHA                   BETA                   VML

           4.923379               9.819798            -262.336670


                  *** GRID SEARCH FOR K =    14 ***

           ALPHA                   BETA                   VML

           5.668086            2980.958008            -266.747559

           5.595870              54.598148            -266.436584

           5.195685              20.085537            -264.737946

           4.685275              12.182494            -262.971497

           4.044650               7.389056            -262.288147

           3.335586               4.481689            -263.126434

           2.496436               2.718282            -262.891663

           1.585763               1.648721            -263.418976

           0.969140               1.000000            -263.164032


                  *** OPTIMAL PARAMETER ESTIMATES ***

           ALPHA                   BETA                   VML

           4.145060               7.966486            -262.260559


                  *** GRID SEARCH FOR K =    16 ***

           ALPHA                   BETA                   VML

           4.970138            2980.958008            -266.419281

           4.924928              54.598148            -266.199646

           4.663393              20.085537            -264.938660

           4.280633              12.182494            -263.266602

           3.741570               7.389056            -262.020355

           3.132969               4.481689            -262.401733

           2.421248               2.718282            -262.555817

           1.586469               1.648721            -262.426025

           0.967658               1.000000            -263.135101


                  *** OPTIMAL PARAMETER ESTIMATES ***

           ALPHA                   BETA                   VML

           3.639074               6.767537            -261.987305


                  *** GRID SEARCH FOR K =    18 ***

           ALPHA                   BETA                   VML

           4.578322            2980.958008            -266.804504

           4.543117              54.598148            -266.619690

           4.336464              20.085537            -265.541168

           4.019334              12.182494            -264.001404

           3.542742               7.389056            -262.540100

           2.990577               4.481689            -262.511810

           2.351537               2.718282            -262.633911

           1.584173               1.648721            -262.158264

           0.966332               1.000000            -262.868408


                  *** OPTIMAL PARAMETER ESTIMATES ***

           ALPHA                   BETA                   VML

           1.695147               1.769263            -262.118530


                  *** GRID SEARCH FOR K =    20 ***

           ALPHA                   BETA                   VML

           4.053934            2980.958008            -266.525970

           4.032835              54.598148            -266.401428

           3.905046              20.085537            -265.648315

           3.687815              12.182494            -264.401672

           3.304344               7.389056            -262.665924

           2.822716               4.481689            -262.080078

           2.252759               2.718282            -262.445251

           1.555777               1.648721            -261.772278

           0.955586               1.000000            -262.617645


                  *** OPTIMAL PARAMETER ESTIMATES ***

           ALPHA                   BETA                   VML

           1.540533               1.631551            -261.771484


              *** THE FINAL SOLUTION     (K =    20) ***

           ALPHA                   BETA                   VML

           1.540533               1.631551            -261.771484


NMISS =    0


                               ST

     1        2        3        4        5        6        7        8

  17.0     42.0     44.0     48.0     60.0     72.0     74.0     95.0


     9       10       11       12       13       14       15       16

 103.0    108.0    122.0    144.0    167.0    170.0    183.0    185.0


    17       18       19       20       21       22       23       24

 193.0    195.0    197.0    208.0    234.0    235.0    254.0    307.0


    25       26       27       28       29       30       31       32

 315.0    401.0    445.0    464.0    484.0    528.0    542.0    567.0


    33       34       35       36       37       38       39       40

 577.0    580.0    795.0    855.0    882.0    892.0   1031.0   1033.0


    41       42       43       44       45

1306.0   1335.0   1366.0   1452.0   1472.0


                                     JCEN

 1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20

 0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0


21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40

 0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   1   1   1   1


41  42  43  44  45

 1   1   1   1   1


                       H

          1        2        3        4        5

 1    217.0    218.0      1.0     21.0      1.0

 2    192.0    193.0      1.0     21.0      0.5

 3    190.0    191.0      1.0     21.0      0.3

                      .

                      .

                      .

43   1173.0   1181.0     17.0     45.0      0.0

44   1259.0   1267.0     17.0     45.0      0.0

45   1279.0   1287.0     17.0     45.0      0.0



http://www.vni.com/
PHONE: 713.784.3131
FAX:713.781.9260