Survival version 2.40 has been relased to CRAN. This is a warning that some users may see
changes in results, however.
The heart of the issue can be shown with a simple example. Calculate the following simple
set of intervals:
<<interval1>>=
birth <- as.Date("1973/03/10")
start <- as.Date("1998/09/13") + 1:40
end <- as.Date("1998/12/03") + rep(1:10, 4)
interval <- (end-start)
table(interval)
51 61 71 81
10 10 10 10
@
Each interval has a different start and end date, but there are only 4 unique intervals,
each of which appears 10 times.
Now convert this to an age scale.
<<interval2>>=
start.age <- as.numeric(start-birth)/365.25
end.age <- as.numeric(end -birth)/365.25
age.interval <- end.age - start.age
table(match(age.interval, unique(age.interval)))
1 2 3 4 5 6 7 8
9 1 5 5 1 9 7 3
@
There are now eight different age intervals instead of 4, and the 8 unique values appear
between 1 and 9 times each. Exact results likely will depend on your computer system. We
have become a victim of round off error.
Some users prefer to use time in days and some prefer time in years, and those latter
users expect, I am sure, survival analysis results to be identical on the two scales.
Both the coxph and survfit routines treat tied event times in a special way, and this
roundoff can make actual ties appear as non-tied values, however. Parametric survival such
as \code{survreg} is not affected by this issue.
In survival version 2.40 this issue has been addressed for the coxph and survfit routines;
input times are subjected to the same logic found in the all.equal routine in order to
determine actual ties. The upshot is that some users may experience a changed results.
For the following test case cox1 and cox2 are identical coefficients in version 2.40, but
different in prior versions.
<<>>=
ndata <- data.frame(id=1:30,
birth.dt = rep(as.Date("1953/03/10"), 30),
enroll.dt= as.Date("1993/03/10") + 1:30,
end.dt = as.Date("1996/10/21") + 1:30 +
rep(1:10, 3),
status= rep(0:1, length=30),
x = 1:30)
ndata$enroll.age <- with(ndata, as.numeric(enroll.dt - birth.dt))/365.25
ndata$end.age <- with(ndata, as.numeric(end.dt - birth.dt))/365.25
fudays <- with(ndata, as.numeric(end.dt - enroll.dt))
fuyrs <- with(ndata, as.numeric(end.age- enroll.age))
cox1 <- coxph(Surv(fudays, status) ~ x, data=ndata)
cox2 <- coxph(Surv(fuyrs, status) ~ x, data=ndata)
@
This general issue of floating point precision arises often enough in R that is part of
the frequently asked questions, see FAQ 7.31 on CRAN. The author of the survival routines
(me) has always used days as the scale
for analysis -- just by habit, not for any particluarly good reason -- so the issue had
never appeared in my work nor in the survival package's test suite. Due to user input,
this issue had been addressed earlier in the survfit routine, but only when the status
variable was 0/1, not when it is a factor.
As a final footnote, the simple data set above also gives different results when coded in
SAS: I am not alone in overlooking it. As a consequence, the maintainer expects to get
new emails that ``we have found a bug in your code: it gives a different answer than
SAS''. (This is an actual quote.)
Terry Therneau
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
