Example 66.4 Firth’s Correction for Monotone Likelihood

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 66.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 monotonity is demonstrated in Output 66.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 "Analyis of Maximum Likelihood Estimates" table reveals a very large value for the coefficient of Contrived. This is very typical of a diverged estimate.

Output 66.4.1 Monotone Likelihood Behavior Displayed
The PHREG Procedure

Maximum Likelihood Iteration History
Iter Ridge Log Likelihood LogBUN HGB Contrived
0 0 -154.8579914384 0.0000000000 0.000000000 0.000000000
1 0 -140.6934052686 1.9948819671 -0.084318519 1.466331269
2 0 -137.7841629416 1.6794678962 -0.109067888 2.778361123
3 0 -136.9711897754 1.7140611684 -0.111564202 3.938095086
4 0 -136.7078932606 1.7181735043 -0.112273248 5.003053568
5 0 -136.6164264879 1.7187547532 -0.112369756 6.027435769
6 0 -136.5835200895 1.7188294108 -0.112382079 7.036444978
7 0 -136.5715152788 1.7188392687 -0.112383700 8.039763533
8 0 -136.5671126045 1.7188405904 -0.112383917 9.040984886
9 0 -136.5654947987 1.7188407687 -0.112383947 10.041434266
10 0 -136.5648998913 1.7188407928 -0.112383950 11.041599592
11 0 -136.5646810709 1.7188407960 -0.112383951 12.041660414
12 0 -136.5646005760 1.7188407965 -0.112383951 13.041682789
13 0 -136.5645709642 1.7188407965 -0.112383951 14.041691020
14 0 -136.5645600707 1.7188407965 -0.112383951 15.041694049
15 0 -136.5645560632 1.7188407965 -0.112383951 16.041695162
16 0 -136.5645545889 1.7188407965 -0.112383951 17.041695572

Analysis of Maximum Likelihood Estimates
Parameter DF Parameter
Estimate
Standard
Error
Chi-Square Pr > ChiSq Hazard
Ratio
LogBUN 1 1.71884 0.58376 8.6697 0.0032 5.578
HGB 1 -0.11238 0.06090 3.4053 0.0650 0.894
Contrived 1 17.04170 1080 0.0002 0.9874 25183399

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 66.4.2). The much preferred profile-likelihood confidence limits, as shown in (Heinze and Schemper; 2001), are also displayed.

Output 66.4.2 Convergence Obtained with the Firth Correction
The PHREG Procedure

Maximum Likelihood Iteration History
Iter Ridge Log Likelihood LogBUN HGB Contrived
0 0 -150.7361197494 0.0000000000 0.000000000 0.0000000000
1 0 -136.9933949142 2.0262484120 -0.086519138 1.4338859318
2 0 -134.5796594364 1.6770836974 -0.109172604 2.6221444778
3 0 -134.1572923217 1.7163408994 -0.111166227 3.4458043289
4 0 -134.1229607193 1.7209210332 -0.112007726 3.7923555412
5 0 -134.1228364805 1.7219588214 -0.112178328 3.8174197804
6 0 -134.1228355256 1.7220081673 -0.112187764 3.8151642206

Analysis of Maximum Likelihood Estimates
Parameter DF Parameter
Estimate
Standard
Error
Chi-Square Pr > ChiSq Hazard
Ratio
95% Hazard Ratio Profile
Likelihood Confidence
Limits
LogBUN 1 1.72201 0.58379 8.7008 0.0032 5.596 1.761 17.231
HGB 1 -0.11219 0.06059 3.4279 0.0641 0.894 0.794 1.007
Contrived 1 3.81516 1.55812 5.9955 0.0143 45.384 5.406 6005.404