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