Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
A routine showing how to calculate quantiles.
library(survival)
library(Hmisc)
# 1. a weibull random number generator, see http://statwonk.com/weibull.html
rweibull_cens <- function(n, shape, scale) {
# will happen to see the death time first or censoring?
rweibull(n, shape = shape, scale = scale) -> a_random_death_time
rweibull(n, shape = shape, scale = scale) -> a_random_censor_time
pmin(a_random_censor_time, a_random_death_time) -> observed_time
(observed_time == a_random_death_time) -> censor
data.frame(time = observed_time, censor = censor, stringsAsFactors = FALSE)
}
# 2. generate sample data
rweibull_cens(1e3, 1, 500) -> d
# 3. calculate a single median 0.5 quantile (percentile/100)
(survival:::quantile.survfit(
survfit(Surv(time, censor) ~ 1, data = d),
probs = 0.5
) -> single_estimate)
# 4. calculate 1,000 of them and visualize a distribution
hist(Hmisc::bootkm(Surv(d$time, d$censor), B = 1e3), breaks = 25, col = "orange",
xlab = "Time until half convert", main = "Time until half convert (median failure time)")
abline(v = single_estimate, col = "black", lwd = 2, lty = 2)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.