How to calculate the confidence interval around difference in median survival time between two (censored) groups? I need the confidence interval around the difference in median survival time in a two-sample test.  I have data for right-censored survival time for a treatment and a control groups, and I want to estimate the magnitude of the difference in median survival time. I have found ways to calculate confidence intervals around one population, or even confidence intervals around the difference in survival of 2 populations at a fixed time, but I'm looking for the confidence interval around difference in median survival time.
Example:
timeControl <- c(16, 17, 18, 31, 32, 33, 34, 35, 36, 37, 40, 45, 50, 81)
censorControl <- rep(1, 14)
timeTreatment <- c(10, 36, 38, 45, 47, 51, 63, 69, 102, 96, 105, 125)
censorTreatment <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0)

Desired output is in the form of a 95% CI, e.g. (x, y).
Are there packages in R or Python that do this calculation that I'm missing? If not, how would I implement the calculation?
Note, I've looked at the R packages: km.ci, controlTest, bpcp, and survival, and the python package lifelines, and they don't appear to do what I'm asking. 
 A: I ended up using survRM2::rmst2() to calculate the confidence interval around the difference in survival times. It doesn't calculate a confidence interval of the median, but it incorporates censoring so I can censor that have a large impact on the mean (e.g. censor events from the longest 20% of events). 
A: The problem of converting two single-sample confidence intervals (CI) into CI for their difference is difficult. The classic example is the Behrens-Fisher problem for inference on differences in means between two normal distributions whose means and variances have been estimated from the data and whose variances are not assumed to be the same.
With Kaplan-Meier (KM) survival curves the problem is compounded by its non-parametric form. The CI around the median survival time for a single KM curve are obtained by calculating the CI for survival fraction as a function of time and finding the points at which the lower and upper CI curves cross 50% survival. With two samples you get two separate medians and CI.
Fay et al. considered the problem of "Combining One-Sample Confidence Procedures for Inference in the Two-Sample Case" in Biometrics 71, 146–156 (2015). They expanded on prior methods that had required assumptions about shifts or continuity in survival curves, and developed a method for the median survival comparison if there is no censoring. That method is provided in the bpcp package, along with methods for comparisons at specific time points even with censoring. But to my knowledge there is yet no way to do what you want with KM curves and censoring.
One alternative is to fit a parametric model to the censored data, in which confidence intervals for parameter-value differences could be translated directly into confidence intervals for median-survival differences. Accelerated failure time models (as suggested in a comment on another answer) could be particularly straightforward, as the parameter-value differences model the relative squeezing of the time scale between the survival curves.
A second alternative would be bootstrapping. Even that might be troublesome, however, as the median-survival difference might not be close enough to pivotal to allow reliable application of simple empirical or percentile bootstrapping approaches.
