The PHREG procedure offers the HAZARDRATIO statement, which you can use to estimate hazard ratios, similar to the ODDSRATIO statement in the LOGISTIC procedure for estimating odds ratios. (See SAS Note 70221 for more information). 

While PHREG does not provide a PLOTS= option to create plots of hazard ratios, the following shows how you can produce plots of hazard ratios in several situations. For similar examples of plotting odds ratios from logistic models fit by PROC LOGISTIC, see SAS KB0044254, "Plot odds ratio change over a continuous variable in polynomial, spline, or interaction effect."

A hazard ratio for a continuous predictor is the change in hazard that results from increasing the predictor by some number of units. It is the ratio of the hazard at the interval maximum over the hazard at the interval minimum. You can specify the interval width in units of the predictor using the UNITS= option in the HAZARDRATIO statement. By default, UNITS=1 produces hazard ratios for a one-unit increase in the predictor. You can specify the starting point of the interval or intervals using the AT= option. For a categorical predictor, the hazard ratio compares the hazards at two levels of the predictor. You can also specify how levels are compared using the DIFF= option. By default, hazard ratios for all pairs of levels are estimated.

A continuous predictor or a binary categorical predictor that does not interact with another variable, is not involved in polynomial (higher-order) effects, and is not represented in the model using a spline has only a single hazard ratio. A multi-level categorical predictor in this scenario has a single set of hazard ratios that compare its levels. However, if the predictor is involved in interactions, splines, or polynomial effects (such as a quadratic effect), then the hazard ratio for the predictor depends on the value of the interacting variable or on its own value if the variable appears in a spline or a higher order effect.

For more information about hazard ratios, see "Hazard Ratios" in the Details section of the PROC PHREG documentation. The following examples use the myeloma data shown in the example titled "Stepwise Regression" of the PHREG procedure documentation to explain how you can produce plots of changing hazard ratios in a variety of models.

This article contains the following sections:

Hazard Ratios for a Binary Variable Interacting with a Continuous Variable

The following statements fit a model in which the binary Platelet variable (with levels 0 and 1) interacts with the continuous Protein variable.

The HAZARDRATIO statement requests hazard ratios comparing the two Platelet levels at each integer value of Protein from 0 to 5. The DIFF=REF option specifies that each hazard ratio has the level 0 hazard in the denominator because 0 is specified as the reference level for Platelet in the CLASS statement. The ODS OUTPUT statement saves the table of hazard ratio estimates in a data set. You can use PROC PRINT to examine the variable names and contents of the saved data set. PROC SGPLOT then plots the hazard ratios. 

proc phreg data=Myeloma;
   class platelet(ref='0');
   model Time*VStatus(0)=protein|platelet;
   hazardratio platelet / at (protein=0 to 5) diff=ref;
   ods output hazardratios=hr;
   run;
proc sgplot noautolegend;
   highlow high=waldupper low=waldlower x=description / lowcap=serif highcap=serif;
   scatter y=hazardratio x=description;
   xaxis label='Protein' grid; 
   yaxis label='Hazard Ratio' grid;
   refline 1 / axis=y;
   run;

The resulting plot shows an increasing hazard ratio for level 1 versus level 0 of Platelet over increasing Protein values. The hazard ratios are significantly less than one for Platelet values below three, suggesting that level 1 significantly reduces the hazard.

Hazard Ratios for a Continuous Variable Interacting with a Continuous Variable

The following model includes two continuous variables, HGB and Age, and their interaction.

In the HAZARDRATIO statement, hazard ratios are requested for HGB at Age values 40, 50, 60, 70, and 80. Since the UNITS= option is not specified, each estimated hazard ratio reflects a one-unit increase in HGB. 

proc phreg data=Myeloma;
      model Time*VStatus(0)=age|hgb ;
      hazardratio hgb / at (age=40 to 80 by 10);
      ods output hazardratios=hr;
      run;
proc sgplot data=hr noautolegend;
      band upper=waldupper lower=waldlower x=description / transparency=.5;
      series y=hazardratio x=description;
      xaxis grid; 
      yaxis label='Hazard Ratio' grid;
      refline 1 / axis=y;
      run;

The resulting plot shows that the one-unit hazard ratios for HGB increase with Age but are only significant at Age=60.

Hazard Ratios for a Polynomial or Splined Variable

Adding a spline or polynomial effect in the model for a continuous predictor allows for a nonlinear fit to the response function. Although polynomial effects (such as quadratic or cubic effects) allow for a variety of shapes, splines are even more flexible. When the predictor is involved in a nonlinear effect, its hazard ratio changes over the range of the predictor. In the following statements, the EFFECT statement produces a natural cubic spline named SPL for the Age predictor. The spline then represents the effect of Age in the model.

The HAZARDRATIO statement requests hazard ratio estimates for increases over one unit intervals of AGE (since UNITS= is not specified). Since the hazard ratio for AGE depends on where those intervals are located in the range of AGE, the AT= option is used to specify the starting points of several one-unit intervals: 45 to 46, 46 to 47, and so on. Rather than using the character Description variable in the saved data set, the SCAN function in the DATA step below captures the values of Age that appear at the end of the Description values and uses the resulting numeric Age variable for the horizontal axis. 

proc phreg data=Myeloma;
      effect spl=spline(age/naturalcubic basis=tpf(noint));
      model Time*VStatus(0)=spl;
      hazardratio age / at (age=45 to 75);
      ods output hazardratios=hr;
      run;
data hr; set hr;
      length Age 8;
      Age=scan(description,-1,'=');
      run;
proc sgplot data=hr noautolegend;
      band upper=waldupper lower=waldlower x=age / transparency=.5;
      series y=hazardratio x=age;
      xaxis grid; 
      yaxis label='Hazard Ratio' grid;
      refline 1 / axis=y;
      title 'Hazard ratios for 1 unit increase in Age';
      run; 

While never significantly different from 1, the plot shows increasing one-unit hazard ratios in the middle of the Age range and plateaus at more extreme age values.

Hazard Ratios for a Splined Variable Interacting with a Continuous Variable

In a model that contains two interacting continuous variables, the hazard ratio for one variable depends on the value of the other. In the case of a splined continuous variable that interacts with another continuous variable, the hazard ratio for the splined variable depends on the values of both variables. In this latter scenario, a plot of the changing hazard ratios becomes three-dimensional, varying both variables with the hazard ratio estimate at each combination. You can present this scenario as a two-dimensional contour plot with contour lines that indicate the hazard ratio values.

As in the spline example above, the EFFECT statement defines a spline on the LOGBUN predictor. The MODEL statement specifies the model that includes the effects of HGB, the spline on LOGBUN, and their interaction. In the HAZARDRATIO statement, the AT= option specifies sets of values for HGB and LOGBUN at which to estimate the hazard ratio for one unit increases in LOGBUN. As mentioned before, the ODS OUTPUT statement saves the estimated hazard ratios in a data set.

proc phreg data=Myeloma;
      effect spl=spline(logbun/naturalcubic basis=tpf(noint));
      model Time*VStatus(0)=spl|hgb ;
      hazardratio logbun / at (hgb=5 to 12 logbun=1 to 2 by .1);
      ods output hazardratios=hr;
      run; 

Since the saved data set only shows the values of HGB and LOGBUN within the character Description variable, the following DATA step adds separate HGB and LOGBUN variables, taking care that they match the values as shown in the Description. Following the DATA step, the TEMPLATE procedure defines a contour plot template named HRcontourplot to produce the desired contour plot. In the CONTOURPLOTPARM statement, LOGBUN is specified as the horizontal axis variable, HGB as the vertical axis variable, and HAZARDRATIO (created by the HAZARDRATIO statement) as the variable that defines the contours. The NHINT=12 option requests approximately twelve contour lines. The CONTOURTYPE=LABELEDLINEFILL option requests labels and colored contour intervals on the contour lines. The CONTINUOUSLEGEND statement provides a scale of colors and the associated hazard ratio values. For more information about the CONTOURPLOTPARM and other statements in PROC TEMPLATE, see SAS^® Graph Template Language: Reference in the SAS documentation.

   data hr2; 
      do hgb=5 to 12;
      do logbun=1 to 2 by .1;
      set hr; output;
      end; end;
      run;
   proc template;
      define statgraph HRcontourplot;
      begingraph;
        entrytitle "Hazard Ratio Plot";
        layout overlay /
            xaxisopts=(offsetmin=0 offsetmax=0
              linearopts=(thresholdmin=0 thresholdmax=0))
            yaxisopts=(offsetmin=0 offsetmax=0
              linearopts=(thresholdmin=0 thresholdmax=0));
          contourplotparm x=logbun y=hgb z=hazardratio /
            contourtype=labeledlinefill nhint=12 name="Contour";
          continuouslegend "Contour" / title="Hazard Ratio";
        endlayout;
      endgraph;
      end;
      run;
   proc sgrender data=hr2 template=HRcontourplot;
      run; 

The resulting contour plot shows that the hazard ratio increases rapidly as LOGBUN increases and HGB decreases.

Hazard Ratios for Splined Variable with Changing Starting Points and Units of Change

Returning to a model on a single splined predictor, this example produces a plot that enables you to see how the hazard ratio changes as you change both the size and the location of the interval over which the hazard ratio is computed. Recall that the hazard ratio on a continuous predictor is the ratio of the hazard at the interval endpoint divided by the hazard at the interval starting point.

The statements below fit a model on a spline of the continuous HGB variable. In the HAZARDRATIO statement, the AT= option varies the interval starting point from 5 to 12, and the UNITS= option varies the width of the intervals from 1 to 3 units of HGB. As mentioned before, the table of hazard ratio estimates is saved in a data set. The DATA step then adds separate HGB and UNITS variables (repeating the values as seen embedded in the Description variable) in the saved data set, making them conveniently available for plotting.

Similar to the previous example, PROC TEMPLATE defines a template named HRcontourplot for producing the contour plot with the starting point of the HGB interval on the horizontal axis, the width of the interval (UNITS) on the vertical axis, and the hazard ratio estimate (HAZARDRATIO) as the variable defining the contours.

   proc phreg data=Myeloma;
      effect spl=spline(hgb/naturalcubic);
      model Time*VStatus(0)=spl ;
      hazardratio hgb / at (hgb=5 to 12) units=1 2 3;
      ods output hazardratios=hr;
      run;
   data hr2; 
      do hgb=5 to 12;
      do units=1 to 3;
      set hr; output;
      end; end;
      run;
   proc template;
      define statgraph HRcontourplot;
      begingraph;
        entrytitle "Hazard Ratio Plot";
        layout overlay /
            xaxisopts=(offsetmin=0 offsetmax=0
              linearopts=(thresholdmin=0 thresholdmax=0))
            yaxisopts=(offsetmin=0 offsetmax=0
              linearopts=(thresholdmin=0 thresholdmax=0));
          contourplotparm x=hgb y=units z=hazardratio /
            contourtype=labeledlinefill nhint=12
              name="Contour";
          continuouslegend "Contour" / title="Hazard Ratio";
        endlayout;
      endgraph;
      end;
      run;
   proc sgrender data=hr2 template=HRcontourplot;
      run; 

The resulting plot shows that the hazard ratios are generally below 1, which indicates a decreasing hazard ratio for increases in HGB. It is perhaps most useful to consider how the hazard ratio changes with HGB for a fixed interval width. For example, to see how the hazard ratio changes when HGB increases by two units, scan along a horizontal line at UNITS=2. Along that line, the hazard ratio is about 1 for a two-unit increase starting at HGB=5, is about 0.75 for a two-unit increase starting at HGB=8, and is between 0.6 and 0.65 for a two-unit HGB intervals starting above ten.