Using SAS ODS Graphics with Example for Generating Kaplan-Meier Curves

Recently, I spent some time on exploring the feature of SAS ODS Graphics. ODS Graphics offers an easy way to generate high-quality statistical graphics without extensive programming using SAS/Graph.  ODS Graphics has been included in almost all of SAS/Stat Procedure.

In the discussions below, the example from SAS Proc Lifetest is used to illustrate the generation of Kaplan-Meier curve using ODS Graphics.

On the basis of SAS Example 49.2 Enhanced Survival Plot and Multiple-Comparison Adjustments, we can run the following program.

proc format;

      value risk 1='ALL' 2='AML-Low Risk' 3='AML-High Risk';

   

   data BMT;

      input Group T Status @@;

      format Group risk.;

      label T='Disease Free Time';

      datalines;

   1 2081 0 1 1602 0 1 1496 0 1 1462 0 1 1433 0

   1 1377 0 1 1330 0 1  996 0 1  226 0 1 1199 0

   1 1111 0 1  530 0 1 1182 0 1 1167 0 1  418 1

   1  383 1 1  276 1 1  104 1 1  609 1 1  172 1

   1  487 1 1  662 1 1  194 1 1  230 1 1  526 1

   1  122 1 1  129 1 1   74 1 1  122 1 1   86 1

   1  466 1 1  192 1 1  109 1 1   55 1 1    1 1

   1  107 1 1  110 1 1  332 1 2 2569 0 2 2506 0

   2 2409 0 2 2218 0 2 1857 0 2 1829 0 2 1562 0

   2 1470 0 2 1363 0 2 1030 0 2  860 0 2 1258 0

   2 2246 0 2 1870 0 2 1799 0 2 1709 0 2 1674 0

   2 1568 0 2 1527 0 2 1324 0 2  957 0 2  932 0

   2  847 0 2  848 0 2 1850 0 2 1843 0 2 1535 0

   2 1447 0 2 1384 0 2  414 1 2 2204 1 2 1063 1

   2  481 1 2  105 1 2  641 1 2  390 1 2  288 1

   2  421 1 2   79 1 2  748 1 2  486 1 2   48 1

   2  272 1 2 1074 1 2  381 1 2   10 1 2   53 1

   2   80 1 2   35 1 2  248 1 2  704 1 2  211 1

   2  219 1 2  606 1 3 2640 0 3 2430 0 3 2252 0

   3 2140 0 3 2133 0 3 1238 0 3 1631 0 3 2024 0

   3 1345 0 3 1136 0 3  845 0 3  422 1 3  162 1

   3   84 1 3  100 1 3    2 1 3   47 1 3  242 1

   3  456 1 3  268 1 3  318 1 3   32 1 3  467 1

   3   47 1 3  390 1 3  183 1 3  105 1 3  115 1

   3  164 1 3   93 1 3  120 1 3   80 1 3  677 1

   3   64 1 3  168 1 3   74 1 3   16 1 3  157 1

   3  625 1 3   48 1 3  273 1 3   63 1 3   76 1

   3  113 1 3  363 1

   ;

ods rtf file="c:\temp\test.doc" style=journal;

ods graphics on;

ods trace on;

proc lifetest data=BMT plots=survival(atrisk=0 to 2500 by 500);

      ods select SurvivalPlot;

      time T * Status(0);

      strata Group / test=logrank adjust=sidak;

      run;

ods trace off;

ods graphics off;

ods rtf close;

In the above statements,

ods rtf file="c:\temp\test.doc" specifies the ODS output as RTF file and the file location and file name. style=journal specifies the use of output in publication quality.

"ods graphics on" to invoke the ODS Graphics

"ods trace on" is not really necessary, but it is useful when we need to know the name/location of the SAS template. If we check the SAS Log Window, we will be able to see that the SAS template for the corresponding graph is Stat.Lifetest.Graphics.ProductLimitSurvival

Output Added:

-------------

Name:       SurvivalPlot

Label:      Survival Curves

Template:   Stat.Lifetest.Graphics.ProductLimitSurvival

Path:       Lifetest.SurvivalPlot

-------------

PROC LIFETEST is invoked to compute the product-limit estimate of the survivor function for each of three risk categories (1='ALL' 2='AML-Low Risk' 3='AML-High Risk'). Using ODS Graphics, you can display the number of subjects at risk in the survival plot. The PLOTS= option requests that the survival curves be plotted, and the ATRISK= suboption specifies the time points at which the at-risk numbers are displayed. In the STRATA statement, the ADJUST=SIDAK option requests the idák multiple-comparison adjustment, and by default, all paired comparisons are carried out.

After the run, we will need to use the following statements to close ODS.

ods trace off;

ods graphics off;

ods rtf close;

The above program will create a high quality figure like below:

Suppose we would like to change the title and labels of the figure above, we can use PROC Template. First thing we need to do is to identify which Template we need to modify. The SAS Log window from ods trace on statement indicates the following template.

Stat.Lifetest.Graphics.ProductLimitSurvival

To open this template, we would need to the pull down manu in SAS window, then choose View -> Results -> View -> Template

From Template window, we will locate the following folder:

SAShelp.tmplmst -> stat -> lifetest -> graphics

We will then see the template called ProductLimitSurvival. Double click this file, we will be able to see the contents in this template file.

The template file can be copied and pasted into the program window. We can then edit the template file to modify the title, labels, legends,…

Suppose we run the following program to generate the Kaplan-Meier curve (upward curve instead of the original downward curve). The upward curve is sometimes easier to be understood. For example, it is useful in the time to event variable where the event is a better outcome (for example, time to recovery, time to release from the hospital, time to immune tolerance (in Hemophilia). By changing the layout of the figure, we would also need to modify the label.

proc lifetest data=BMT plots=survival(f test);

      ods select failurePlot;

      time T * Status(0);

      strata Group / test=logrank adjust=sidak;

run;

plots=survival(f test) requests the failure Plot (upward) instead of the survival plot (downward). The ‘test’ indicates the display of the logrank test results.

From the LOG window trigged by ODS Trace On statement, we will be able to know that the template is Stat.Lifetest.Graphics.ProductLimitFailure which is located in the same folder as Stat.Lifetest.Graphics.ProductLimitSurvival discussed before.

Output Added:

-------------

Name:       FailurePlot

Label:      Failure Curves

Template:   Stat.Lifetest.Graphics.ProductLimitFailure

Path:       Lifetest.FailurePlot 

We can copy and paste the template statements and run as part of the program.

ods rtf file="c:\temp\test.doc" style=journal;

ods graphics on;

ods trace on;

proc template;                                                               

   define statgraph Stat.Lifetest.Graphics.ProductLimitFailure;              

      dynamic NStrata xName maxTime plotAtRisk plotCensored plotCL plotHW    

         plotEP labelCL labelHW labelEP yMin StratumID classAtRisk plotTest  

         GroupName Transparency SecondTitle TestName pValue;                 

      BeginGraph;                                                             

         if (NSTRATA=1)                                                      

            if (EXISTS(STRATUMID))                                           

            entrytitle "Kaplan-Meier Curve" " for " STRATUMID;       /*Revised the Figure title*/

         else                                                                

            entrytitle "Kaplan-Meier Curve";                        

         endif;                                                               

         if (PLOTATRISK)                                                     

            entrytitle "with Number of Subjects at Risk" / textattrs=        

            GRAPHVALUETEXT;                                                   

         endif;                                                              

         layout overlay / xaxisopts=(label="Time to Recovery" offsetmin=.05 linearopts /*xaxisopts=(label=" ") to change the x-axis label*/

            =(viewmax=MAXTIME)) yaxisopts=(label="Probability of Achieving Tolerance"       

            shortlabel="Failure" linearopts=(viewmin=0 viewmax=1 tickvaluelist

            =(0 .2 .4 .6 .8 1.0)));                                          

            if (PLOTHW=1 AND PLOTEP=0)                                       

               bandplot LimitUpper=eval (1-HW_LCL) LimitLower=eval (1-HW_UCL)

               x=TIME / modelname="Failure" fillattrs=GRAPHCONFIDENCE name=  

               "HW" legendlabel=LABELHW;                                      

            endif;                                                           

            if (PLOTHW=0 AND PLOTEP=1)                                       

               bandplot LimitUpper=eval (1-EP_LCL) LimitLower=eval (1-EP_UCL)

               x=TIME / modelname="Failure" fillattrs=GRAPHCONFIDENCE name=  

               "EP" legendlabel=LABELEP;                                     

            endif;                                                           

            if (PLOTHW=1 AND PLOTEP=1)                                       

               bandplot LimitUpper=eval (1-HW_LCL) LimitLower=eval (1-HW_UCL)

               x=TIME / modelname="Failure" fillattrs=GRAPHDATA1             

               datatransparency=.55 name="HW" legendlabel=LABELHW;           

            bandplot LimitUpper=eval (1-EP_LCL) LimitLower=eval (1-EP_UCL) x=

               TIME / modelname="Failure" fillattrs=GRAPHDATA2               

               datatransparency=.55 name="EP" legendlabel=LABELEP;           

            endif;                                                           

            if (PLOTCL=1)                                                    

               if (PLOTHW=1 OR PLOTEP=1)                                      

               bandplot LimitUpper=eval (1-SDF_LCL) LimitLower=eval (1-SDF_UCL

               ) x=TIME / modelname="Failure" display=(outline) outlineattrs=

               GRAPHPREDICTIONLIMITS name="CL" legendlabel=LABELCL;           

            else                                                             

               bandplot LimitUpper=eval (1-SDF_LCL) LimitLower=eval (1-SDF_UCL

               ) x=TIME / modelname="Failure" fillattrs=GRAPHCONFIDENCE name=

               "CL" legendlabel=LABELCL;                                     

            endif;                                                           

            endif;                                                           

            stepplot y=eval (1-SURVIVAL) x=TIME / name="Failure" rolename=(  

               _tip1=ATRISK _tip2=EVENT) tip=(y x Time _tip1 _tip2)          

               legendlabel="Failure";                                        

            if (PLOTCENSORED)                                                 

               scatterplot y=eval (1-CENSORED) x=TIME / markerattrs=(symbol= 

               plus) name="Censored" legendlabel="Censored";                 

            endif;                                                           

            if (PLOTCL=1 OR PLOTHW=1 OR PLOTEP=1)                            

               discretelegend "Censored" "CL" "HW" "EP" / location=outside   

               halign=center;                                                 

            else                                                             

               if (PLOTCENSORED=1)                                           

               discretelegend "Censored" / location=inside autoalign=(topleft

               bottomright);                                                 

            endif;                                                           

            endif;                                                           

            if (PLOTATRISK=1)                                                

               innermargin / align=bottom;                                   

               blockplot x=TATRISK block=ATRISK / repeatedvalues=true display=

                  (values) valuehalign=start valuefitpolicy=truncate         

                  labelposition=left labelattrs=GRAPHVALUETEXT valueattrs=   

                  GRAPHDATATEXT (size=7pt) includemissingclass=false;        

            endinnermargin;                                                  

            endif;                                                           

         endlayout;                                                           

         else                                                                

            entrytitle "Kaplan-Meier Curve";                       

         if (EXISTS(SECONDTITLE))                                            

            entrytitle SECONDTITLE / textattrs=GRAPHVALUETEXT;               

         endif;                                                              

         layout overlay / xaxisopts=(label="Time to Recovery" offsetmin=.05 linearopts

            =(viewmax=MAXTIME)) yaxisopts=(label="Probability of Recovery"     /*y-axis label*/  

            shortlabel="Failure" linearopts=(viewmin=0 viewmax=1 tickvaluelist

            =(0 .2 .4 .6 .8 1.0)));                                          

            if (PLOTHW=1)                                                    

               bandplot LimitUpper=eval (1-HW_LCL) LimitLower=eval (1-HW_UCL)

               x=TIME / group=STRATUM index=STRATUMNUM modelname="Failure"   

               datatransparency=Transparency;                                

            endif;                                                           

            if (PLOTEP=1)                                                    

               bandplot LimitUpper=eval (1-EP_LCL) LimitLower=eval (1-EP_UCL)

               x=TIME / group=STRATUM index=STRATUMNUM modelname="Failure"   

               datatransparency=Transparency;                                

            endif;                                                           

            if (PLOTCL=1)                                                    

               if (PLOTHW=1 OR PLOTEP=1)                                     

               bandplot LimitUpper=eval (1-SDF_LCL) LimitLower=eval (1-SDF_UCL

               ) x=TIME / group=STRATUM index=STRATUMNUM modelname="Failure" 

               display=(outline);                                            

            else                                                             

               bandplot LimitUpper=eval (1-SDF_UCL) LimitLower=eval (1-SDF_LCL

               ) x=TIME / group=STRATUM index=STRATUMNUM modelname="Failure" 

               datatransparency=Transparency;                                

            endif;                                                            

            endif;                                                           

            stepplot y=eval (1-SURVIVAL) x=TIME / group=STRATUM index=       

               STRATUMNUM name="Failure" rolename=(_tip1=ATRISK _tip2=EVENT) 

               tip=(y x Time _tip1 _tip2);                                   

            if (PLOTCENSORED)                                                

               scatterplot y=eval (1-CENSORED) x=TIME / group=STRATUM index= 

               STRATUMNUM markerattrs=(symbol=plus);                         

            endif;                                                           

            if (PLOTATRISK)                                                  

               innermargin / align=bottom;                                   

               blockplot x=TATRISK block=ATRISK / class=CLASSATRISK          

                  repeatedvalues=true display=(label values) valuehalign=start

                  valuefitpolicy=truncate labelposition=left labelattrs=     

                  GRAPHVALUETEXT valueattrs=GRAPHDATATEXT (size=7pt)         

                  includemissingclass=false;                                 

            endinnermargin;                                                   

            endif;                                                           

            DiscreteLegend "Failure" / title="Risk Group" location=inside autoalign=(bottomright); /*Revise the legend*/   

            if (PLOTCENSORED)                                                

               if (PLOTTEST)                                                 

               layout gridded / rows=2 autoalign=(TOPLEFT BOTTOMRIGHT BOTTOM 

               TOP) border=true BackgroundColor=GraphWalls:Color Opaque=true;

               entry "+ Censored";                                           

               if (PVALUE < .0001)                                           

                  entry TESTNAME " p " eval (PUT(PVALUE, PVALUE6.4));        

               else                                                          

                  entry TESTNAME " p=" eval (PUT(PVALUE, PVALUE6.4));        

               endif;                                                        

            endlayout;                                                       

            else                                                             

               layout gridded / rows=1 autoalign=(TOPLEFT BOTTOMRIGHT BOTTOM 

               TOP) border=true BackgroundColor=GraphWalls:Color Opaque=true;

               entry "+ Censored";                                           

            endlayout;                                                        

            endif;                                                           

            else                                                             

               if (PLOTTEST)                                                  

               layout gridded / rows=1 autoalign=(TOPLEFT BOTTOMRIGHT BOTTOM 

               TOP) border=true BackgroundColor=GraphWalls:Color Opaque=true;

               if (PVALUE < .0001)                                           

                  entry TESTNAME " p " eval (PUT(PVALUE, PVALUE6.4));        

               else                                                          

                  entry TESTNAME " p=" eval (PUT(PVALUE, PVALUE6.4));        

               endif;                                                        

            endlayout;                                                       

            endif;                                                           

            endif;                                                            

         endlayout;                                                          

         endif;                                                              

      EndGraph;                                                              

   end;                                                                      

run; 

proc lifetest data=BMT plots=survival(f test);

      ods select quartiles FailurePlot;

      time T * Status(0);

      strata Group / test=logrank adjust=sidak;

      run;

ods trace off;

ods graphics off;

ods rtf close;

At the end of the program, we use the following program to delete the user defined Template.

proc template;

   delete Stat.Lifetest.Graphics.ProductLimitFailure;

run;

Further reading:

• Resource for ODS Statistics Graphics
• Getting Started with ODS Statistical Graphics in SAS® 9.2 presents the fundamental information you need to get started with ODS Statistical Graphics in SAS 9.2. (Revised 2009)
• An Overview of ODS Statistical Graphics in SAS 9.3
• Summary of ODS Graphics Capabilities provides a concise overview of the capabilities of ODS Graphics.
• The Graph Template Language and the Statistical Graphics Procedures:
• An Example-Driven Introduction
• Customizing ODS Statistical Graphs
• Customizing ODS Statistical Graphs Step-by-Step
• Modifying Graph Titles and Axis Labels
• Graph Templates

          

Using Area Under the Curve (AUC) as Clinical Endpoints

Area Under Curve (AUC) has been frequently used as the endpoint measure in clinical trials. We use AUC commonly in clinical pharmacology - Area under the time concentration curve or in diagnostic research – Area Under the ROC curve. The use of AUC is much more broader than what we think. Many clinical endpoints can utilize the AUC as a measure for the aggregate effect over a period of time. Below are some of the examples that I have experienced where AUC is used in clinical trials not for the purpose of pharmacokinetics measure or ROC measure.  

AUC Used in Pain Assessment

In the study of pain medications (usually acute pain medications), the pain intensity or pain relief in scales are measured at pre and serial time points post analgesic drug administration. The Summed Pain Intensity Difference (SPID) and total pain relief (TOTPAR) are usually calculated and used as the efficacy endpoints. TOTPAR is a time-weighted measure of AUC or total area under the pain relief curve and is a summary measure that integrates serial assessments of a subject’s pain over the duration of the study. The area under the pain relief vs. time curve can be used to derive the proportion of patients experiencing typically 50% pain relief over a specified time frame. This can be calculated as the ratio of two AUCs: TOTPAR vs. maxTOTPAR (maximum potential value for TOTPAR) as illustrated in the “Analysis of scale results – summary measures” of pain.

FEV1 AUC

In Asthma and COPD studies, FEV1 can be measured at pre-dose and at several serial time points post the treatment. The AUC will then be calculated from the time-FEV1 curve.

In  DULERA drug label, FEV1 AUC(0-12 hr) was mentioned as the efficacy measure:

“FEV1 AUC (0-12hr) was assessed as a co-primary efficacy endpoint to evaluate the contribution of the formoterol component to DULERA. Patients receiving DULERA 100 mcg/5 mcg had significantly higher increases from baseline at Week 12 in mean FEV1 AUC (0-12 hr) compared to mometasone furoate 100 mcg (the primary treatment comparison) and vs. placebo ......"

In a recent news release “Results of Phase II Study of Boehringer Ingelheim's Investigational Bronchodilator for COPD Presented at 2012 ATS International Conference”, FEV1 AUC was used to measure the treatment effect in COPD.

“Results of the study found olodaterol 5 microgram QD provided significant improvement in lung function as measured by FEV1 AUC(0-12) versus twice-daily olodaterol 2 microgram, while twice-daily dosing of olodaterol 5 microgram had a better FEV1 AUC(0-12) profile versus once-daily olodaterol 10 microgram”

AUC in Type 1 Diabetes

In type-1 diabetes research, the main purpose of the treatment is to preserve the beta-cell function. The assessment of beta-cell function is through the measurement of the C-peptide concentration after simulated Mixed Meal Tolerance Test (MMTT) - the gold standard measure of endogenous insulin secretion

In the mixed-meal tolerance test (MMTT), commonly used in the U.S., a liquid meal (Sustacal/Boost) is ingested in the fasting state with timed measurements of C-peptide over the subsequent 2–4 h. The AUC is then calculated for the area under time-C-Peptide curve over 2 hour (AUC_0-2hr) or 4 hours (AUC_0-4hr) (see Greenbaum at al “Mixed-Meal Tolerance Test Versus Glucagon Stimulation Test for the Assessment of β-Cell Function in Therapeutic Trials in Type 1 Diabetes”).

In type 1 diabetes research, a concept of mean AUC is also used. Mean AUC is calculated by the AUC divided by the time duration (i.e., AUC_{0-2 hr} / 120 minutes or AUC_{0-4 hr} / 240 minutes)

AUCs to Assess the Responsiveness

We recently published a paper “Vigorimeter grip strength in CIDP: a responsive tool that rapidly measures the effect of IVIG – the ICE study” where we used AUCs to compare the responsiveness of two difference measures. Since two different measures used different scales, we had to calculate the SRM (standardized response mean) before we calculated the AUCs for INCAT scale and for Grip Strength . The larger the AUC, the higher the responsiveness to the treatment. The results indicated that the Vigorimeter grip strength could be more sensitive measure comparing to INCAT scale to evaluate the treatment effect of IVIG in CIDP patients.

AUC for Visual Analog Scale (VAS) for Dyspnea in Acute Heart Failure

In FDA's Cardiovascular and Renal Drugs Advisory Committee Meeting in March 27, 2014 for Serelaxin for Acute Heart Failure, one of the statistical issues discussed was the use of VAS AUC to assess the dyspnea in Acute Heart Failure. The FDA presentation included the detail calculation of the VAS-AUC and the results from this endpoint.

Missingness Mechanism (MCAR, MAR, and MNAR) - A Great Explanation of These Terms

For statisticians working in clinical trial field, the best challenge may not be in the statistical methodologies. The best challenge may be in communication with non-statisticians (such as physicians, clinical team members, corporate executives) about the statistical concepts and the statistical terminologies in plain languages.  

The missing data is very common in clinical trials and the concept of the missing data is very easy to understand. However, the categories for missing data mechanisms (or taxonomy of missingness) are not so easy to understand. A formal taxonomy exists for classifying missing data mechanisms, including for longitudinal and event history data. The mechanisms can be classified as MCAR (missing completely at random), MAR (missing at random), and MNAR (missing not at random). Take a look at the definition of MCAR, MAR, and MNAR below, you will see that these definitions are not easy to be understood by non-statisticians.

	EMA	Academy (US)
MCAR	For the dependent variable (conditional on the covariates in the model), if the probability of an observation being missing does not depend on observed or unobserved measurements then the observation is Missing Completely At Random (MCAR).	In the case of MCAR, the missing data are unrelated to the study variables: thus, the participants with completely observed data are in effect a random sample of all the participants assigned a particular intervention. With MCAR, the random assignment of treatments is assumed to be preserved, but that is usually an unrealistically strong assumption in practice.
MAR	Conditional on the covariates in the model, if the probability of an observation being missing depends only on observed measurements then the observation is Missing At Random (MAR).	In the case of MAR, whether or not data are missing may depend on the values of the observed study variables. However, after conditioning on this information, whether or not data are missing does not depend on the values of the missing data.
MNAR	When observations are neither MCAR nor MAR, they are classified as Missing Not At Random (MNAR), i.e. the probability of an observation being missing depends on unobserved measurements. In this scenario, the value of the unobserved responses depends on information not available for the analysis (i.e. not the values observed previously on the analysis variable or the covariates being used), and thus, future observations cannot be predicted without bias by the model.	In the case of MNAR, whether or not data are missing depends on the values of the missing data.

Thanks to Ziad Taib, the following example for three different missingness mechanisms were explained very well and were easy to be understood by the non-statisticians.  

Suppose you are modelling weight (Y) as a function of sex (X). Some respondents wouldn't disclose their weight, so you are missing some values for Y. There are three possible mechanisms for the nondisclosure:

• There may be no particular reason why some respondents told you their weights and others didn't. That is, the probability that Y is missing may has no relationship to X or Y. In this case our data is missing completely at random (MCAR)
• One sex may be less likely to disclose its weight. That is, the probability that Y is missing depends only on the value of X. Such data are missing at random (MAR)
• Heavy (or light) people may be less likely to disclose their weight. That is, the probability that Y is missing depends on the unobserved value of Y itself. Such data are not missing at random or missing not at random (MNAR)

 

Understanding the concept of missing mechanism is one thing, fully understanding missing mechanism in practice is another story. The reason for missing data is often not collected or incompletely collected in the clinical trials. Patients may not tell the real reason for them to withdraw from the study (discontinue from the study earlier). Academy’s suggestions below are reasonable, however, ‘full and detailed documentation for each individual of the reasons for missing records or missing observations’ is not the reality in the current clinical trial practice.

"Reasons for missing data must be documented as much as possible. This includes full and detailed documentation for each individual of the reasons for missing records or missing observations. Knowing the reason for missingness permits formulation of sensible assumptions about observations that are missing, including whether those observations are well defined.

Missing data in clinical trials can seriously undermine the benefits provided by randomization into control and treatment groups. Two approaches to the problem are to reduce the frequency of missing data in the first place and to use appropriate statistical techniques that account for the missing data. The former approach is preferred, since the choice of statistical method requires unverifiable assumptions concerning the mechanism that causes the missing data, and so always involves some degree of subjectivity.”