The PHREG Procedure

Example 86.4 Firth’s Correction for Monotone Likelihood

(View the complete code for this example.)

In fitting the Cox regression model by maximizing the partial likelihood, the estimate of an explanatory variable X will be infinite if the value of X at each uncensored failure time is the largest of all the values of X in the risk set at that time (Tsiatis 1981; Bryson and Johnson 1981). You can exploit this information to artificially create a data set that has the condition of monotone likelihood for the Cox regression. The following DATA step modifies the Myeloma data in Example 86.1 to create a new explanatory variable, Contrived, which has the value 1 if the observed time is less than or equal to 65 and has the value 0 otherwise. The phenomenon of monotone likelihood will be demonstrated in the new data set Myeloma2.

data Myeloma2;
   set Myeloma;
   Contrived= (Time <= 65);
run;

For illustration purposes, consider a Cox model with three explanatory variables, one of which is the variable Contrived. The following statements invoke PROC PHREG to perform the Cox regression. The IPRINT option is specified in the MODEL statement to print the iteration history of the optimization.

proc phreg data=Myeloma2;
   model Time*Vstatus(0)=LOGbun HGB Contrived / itprint;
run;

The symptom of monotonicity is demonstrated in Output 86.4.1. The log likelihood converges to the value of –136.56 while the coefficient for Contrived diverges. Although the Newton-Raphson optimization process did not fail, it is obvious that convergence is questionable. A close examination of the standard errors in the "Analysis of Maximum Likelihood Estimates" table reveals a very large value for the coefficient of Contrived. This is very typical of a diverged estimate.

Output 86.4.1: Monotone Likelihood Behavior Displayed

The PHREG Procedure

Maximum Likelihood Iteration History
IterRidgeLog LikelihoodLogBUNHGBContrived
00-154.85799143840.00000000000.0000000000.000000000
10-140.69340526861.9948819671-0.0843185191.466331269
20-137.78416294161.6794678962-0.1090678882.778361123
30-136.97118977541.7140611684-0.1115642023.938095086
40-136.70789326061.7181735043-0.1122732485.003053568
50-136.61642648791.7187547532-0.1123697566.027435769
60-136.58352008951.7188294108-0.1123820797.036444978
70-136.57151527881.7188392687-0.1123837008.039763533
80-136.56711260451.7188405904-0.1123839179.040984886
90-136.56549479871.7188407687-0.11238394710.041434266
100-136.56489989131.7188407928-0.11238395011.041599592
110-136.56468107091.7188407960-0.11238395112.041660414
120-136.56460057601.7188407965-0.11238395113.041682789
130-136.56457096421.7188407965-0.11238395114.041691020
140-136.56456007071.7188407965-0.11238395115.041694049
150-136.56455606321.7188407965-0.11238395116.041695162
160-136.56455458891.7188407965-0.11238395117.041695572

Analysis of Maximum Likelihood Estimates
ParameterDFParameter
Estimate
Standard
Error
Chi-SquarePr > ChiSqHazard
Ratio
LogBUN11.718840.583768.66970.00325.578
HGB1-0.112380.060903.40530.06500.894
Contrived117.0417010800.00020.987425183399


Next, the Firth correction was applied as shown in the following statements. Also, the profile-likelihood confidence limits for the hazard ratios are requested by using the RISKLIMITS=PL option.

proc phreg data=Myeloma2;
   model Time*Vstatus(0)=LogBUN HGB Contrived /
      firth risklimits=pl itprint;
run;

PROC PHREG uses the penalized likelihood maximum to obtain a finite estimate for the coefficient of Contrived (Output 86.4.2). The much preferred profile-likelihood confidence limits, as shown in (Heinze and Schemper 2001), are also displayed.

Output 86.4.2: Convergence Obtained with the Firth Correction

The PHREG Procedure

Maximum Likelihood Iteration History
IterRidgeLog LikelihoodLogBUNHGBContrived
00-150.73611974940.00000000000.0000000000.0000000000
10-136.99339491422.0262484120-0.0865191381.4338859318
20-134.57965943641.6770836974-0.1091726042.6221444778
30-134.15729232171.7163408994-0.1111662273.4458043289
40-134.12296071931.7209210332-0.1120077263.7923555412
50-134.12283648051.7219588214-0.1121783283.8174197804
60-134.12283552561.7220081673-0.1121877643.8151642206

Analysis of Maximum Likelihood Estimates
ParameterDFParameter
Estimate
Standard
Error
Chi-SquarePr > ChiSqHazard
Ratio
95% Hazard Ratio Profile
Likelihood Confidence
Limits
LogBUN11.722010.583798.70080.00325.5961.76117.231
HGB1-0.112190.060593.42790.06410.8940.7941.007
Contrived13.815161.558125.99550.014345.3845.4066005.404