JuliaStats / Survival.jl

Survival analysis in Julia
MIT License
73 stars 23 forks source link

Possible bug in Nelson-Aalen #50

Open itsdfish opened 2 years ago

itsdfish commented 2 years ago

I noticed that NelsonAalen repeats the second to last value twice. I thought this was odd, so I compared it to R.

Julia

using Survival 

x = 1:5
v = fill(1, 5)

result = fit(NelsonAalen, x, v)
result.chaz

Result

5-element Vector{Float64}:
 0.2
 0.45
 0.7833333333333333
 1.2833333333333332
 1.2833333333333332

R

library("mice")
df = data.frame(time = c(1,2,3,4,5))
df$status = 1
nelsonaalen(df, time, status)

Result

[1] 0.2000000 0.4500000 0.7833333 1.2833333 2.2833333

Version

Julia 1.8.1
Survival 0.2.2
ararslan commented 2 years ago

I'm not familiar with the mice library in R. Does the same thing happen using the survival library?

itsdfish commented 2 years ago

Unfortunately, I had trouble figuring out how to adapt the example above to surival in R. Here is what I tried, but the output didn't look correct:

library("survival")
df = data.frame(time = c(1,2,3,4,5))
df$status = 1
fit <- survfit(Surv(time, status) ~ 1, ctype=1, data = df)

Do you have more experience with it?

ararslan commented 2 years ago

AFAIK there are two ways, one more manual than the other:

library("survival")
df <- data.frame(time=1:5, status=1)
na1 <- basehaz(coxph(Surv(time, status) ~ 1, data=df))$hazard
surv <- survfit(Surv(time, status) ~ 1, data=df)
na2 <- cumsum(surv$n.event / surv$n.risk)

They produce the same result, which matches your mice example.

However, when using other data (specifically t and s as defined in the Nelson-Aalen portion of test/runtests.jl), I get matching results for Julia and R.