1
$\begingroup$

Hello I am learning about survival analysis and I noticed that the variance estimate from R's survival package is different from when I calculate variance by hand using green wood's formula

$$\widehat{Var}[\hat{S}(t)] = \hat{S}(t)^2 \sum_{t_i \leq t}\frac{d_i}{r_i(r_i - d_i)} $$

Where $r_i$ is the number at risk right before $t$ and $d_i$ is the number of deaths occurred at time $t_i$

library(survival)
library(dplyr)
library(broom)
size = 10
deathtime <- seq(1, size)
death = rep(1, size)
df <- data.frame(deathtime, death)

surv.obj<- Surv(df$deathtime, df$death)

surv.fit <- survfit(surv.obj ~ 1)

df.fit <- tidy(surv.fit)


df.fit %>% 
  mutate(var.hand = estimate^2 * cumsum(n.event/(n.risk * (n.risk - n.event))), # Green Wood
         var.survival = std.error^2) #Survival Package

enter image description here

The var.hand calculated varies very differently then the variance calculated from the survival package I am curious why?

$\endgroup$

1 Answer 1

1
$\begingroup$

The std.err reported by survfit() is the standard error of the cumulative hazard. From the help page for survfit.object:

for a survival curve this contains standard error of the cumulative hazard or -log(survival)

That's just this part of the Greenwood formula:

$$ \sum_{t_i \leq t}\frac{d_i}{r_i(r_i - d_i)}.$$

You must multiply the corresponding variance of cumulative hazard by $\hat S(t)^2$ to get the variance in survival at time $t$.

survSquared <- seq(.9,0,-.1)^2
survSquared
 [1] 0.81 0.64 0.49 0.36 0.25 0.16 0.09 0.04 0.01 0.00
surv.fit <- survfit(surv.obj~1)
surv.fit$std.err
 [1] 0.1054 0.1581 0.2070 0.2582 0.3162 0.3873 0.4830 0.6325 0.9487    Inf
varCumSurv <- surv.fit$std.err^2
survSquared * varCumSurv
 [1] 0.009 0.016 0.021 0.024 0.025 0.024 0.021 0.016 0.009   NaN

Those values agree with yours.

$\endgroup$
1

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.