With the help of ChatGPT, I review the code in survival package (survival source: R/quantile.survfit.R). A portion of the description matches what the code demonstrates.
# quantile = where a horzontal line at p intercects the curve. At each
# x the curve of 1-y jumps up to a new level
# The most work is to check for horizontal lines in the survival
# curve that match one of our quantiles within tolerance. If any
# p matches, then our quantile is the average of the given x and
# the x value of the next jump point, i.e., the usual midpoint rule
# used for medians.
# A flat at the end of the curve is a special case, as is the quantile
# of 0.
indx1 <- approx(y+tol, 1:n, p, method="constant", f=1)$y
indx2 <- approx(y-tol, 1:n, p, method="constant", f=1)$y
quant <- (x[indx1] + x[indx2])/2
quant[p==0] <- x[1]
if (!is.na(y[n])) {
lastpt <- (abs(p- y[n]) < tol) # end of the curve
if (any(lastpt)) quant[lastpt] <- (x[indx1[lastpt]] + xmax)/2
}
The results of your code were
> quantile(kmfit, probs = c(0.25, 0.50, 0.75), conf.int = FALSE)
25 50 75
6.93 18.86 NA
Also I run the summary command
> summary(kmfit)
Call: survfit(formula = Surv(Duration, Response == 0) ~ 1, data = dt,
conf.int = 0.95, conf.type = "log-log")
time n.risk n.event survival std.err lower 95% CI upper 95% CI
4.14 8 1 0.875 0.117 0.387 0.981
5.00 7 1 0.750 0.153 0.315 0.931
8.86 6 1 0.625 0.171 0.229 0.861
12.29 5 1 0.500 0.177 0.152 0.775
And you could see that 18.86 equal to
12.29/2+25.43/2
I could not run the SAS code that you provide, so I rewrite the code here. With the code, you could see the plot with 95% confidence interval.
data dt;
input Duration Response trt;
datalines;
4.14 0 1
12.29 0 1
16.57 1 1
15.43 1 1
5.00 0 1
15.71 1 1
8.86 0 1
25.43 1 1
;
run;
proc lifetest data=dt conftype=loglog plots=survival(cl);
time Duration * Response(1);
strata trt;
run;
The tables would be
As why the table did not show the median value, I read the SAS document (https://support.sas.com/documentation/onlinedoc/stat/123/lifetest.pdf) and it said
Then I redefined the outcome event .
The R code would be
kmfit <- survfit(
Surv(Duration, Response == 1) ~ 1,
data = dt,
conf.int = 0.95,
conf.type = "log-log"
)
quantile(kmfit, probs = c(0.25, 0.50, 0.75), conf.int = FALSE)
## the output
25 50 75
15.57 16.14 21.00
The SAS code would be
proc lifetest data=dt conftype=loglog plots=survival(cl);
time Duration * Response(0);
strata trt;
run;
Please see
75 21.0000 LOGLOG 15.4300 .
50 16.1400 LOGLOG 15.4300 .
25 15.5700 LOGLOG 15.4300 16.5700
Not surprisingly, they came out the same.
