Example 7.40: Nelson-Aalen plotting

In our previous entry, we described how to calculate the Nelson-Aalen estimate of cumulative hazard.

In this entry, we display the estimates for the time to linkage to primary care for both the treatment and control groups in the HELP study.

R

We use the previously defined function, after removing missing values and sorting by the time to event or censoring.

ds = read.csv("http://www.math.smith.edu/sasr/datasets/help.csv")
smallds = data.frame(dayslink=ds$dayslink, 
   linkstatus=ds$linkstatus, treat=ds$treat)

# drop subjects with missing data
smallds = na.omit(smallds)

# order by dayslink
smallds = smallds[order(smallds$dayslink),]

rm(ds)  # clean up
attach(smallds)
library(survival)
calcna = function(time, event) {
   na.fit = survfit(coxph(Surv(time,event)~1),
     type="aalen")
   jumps = c(0, na.fit$time, max(time))
   # need to be careful at the beginning and end
   surv = c(1, na.fit$surv, na.fit$surv[length(na.fit$surv)])

   # apply appropriate transformation
   neglogsurv = -log(surv)   
   
   # create placeholder of correct length
   naest = numeric(length(time))  
   for (i in 2:length(jumps)) {
      naest[which(time>=jumps[i-1] & time<=jumps[i])] = 
         neglogsurv[i-1]   # snag the appropriate value
   }
   return(naest)
}

nacontrol = calcna(dayslink[treat==0], linkstatus[treat==0])
natreat = calcna(dayslink[treat==1], linkstatus[treat==1])

plot(dayslink[treat==1], natreat, type="s", col="blue",
   ylab="Nelson-Aalen estimate",
   xlab="Number of days", lwd=2)
lines(dayslink[treat==0], nacontrol, lty=2, lwd=2,
  col="red", type="s")
legend(250, 0.55, legend=c("Intervention", "Control"), 
   lty=1:2, lwd=2, col=c("blue", "red"))

The time to linkage was much lower for the treatment group than for the control group.

SAS

In SAS we'll use proc lifetest with the strata statement, as in section 5.6.3, removing subjects with missing time, censoring, or treatment indicators using the nmiss function (section 1.4.14) and subsetting if statement (section 1.5.1). We supress all printed output and save the desired data set using ODS commands.

proc import file="c:\book\help.csv"
   out=help dbms=dlm;
   delimiter=',';
   getnames=yes;
run;

data h2;
set help;
  if nmiss(dayslink, linkstatus,treat) eq 0;
run;

ods select none;
ods output productlimitestimates=naout1;
proc lifetest data=h2 nelson;
  time dayslink*linkstatus(0);
  strata treat;
run;
ods select all;

Then we just have to plot the data. We make it a little more self-explanatory by defining a format for the treatment group. Separate curves are requested with the y*x=z syntax, which also produces a legend by default. We use the symbol statement with the stepj interpolation (section 5.1.19) to ask for a step function plot of each curve. The axis statement is used to rotate the title of the y-axis. The label for the y-axis is inherited from proc lifetest.

proc format;
  value treat 0="Control" 1="Intervention";
run;

axis1 label = (angle = 90) minor = none order = 0 to 1.1 by .1;
symbol1 i = stepj l = 3 w = 5 c = red;
symbol2 i = stepj l = 1 w = 5 c = blue;
proc gplot data = naout1;
  plot cumhaz * dayslink = treat / vaxis = axis1;
  label treat = "Treatment group";
  format treat treat. cumhaz dayslink 4.1;
run;
quit;