[R] Survival Analysis with an Historical Control

[Therneau, Terry M., Ph.D.]

You are asking for a one sample test.  Using your own data:

connection <- textConnection("
GD2  1   8 12  GD2  3 -12 10  GD2  6 -52  7
GD2  7  28 10  GD2  8  44  6  GD2 10  14  8
GD2 12   3  8  GD2 14 -52  9  GD2 15  35 11
GD2 18   6 13  GD2 20  12  7  GD2 23  -7 13
GD2 24 -52  9  GD2 26 -52 12  GD2 28  36 13
GD2 31 -52  8  GD2 33   9 10  GD2 34 -11 16
GD2 36 -52  6  GD2 39  15 14  GD2 40  13 13
GD2 42  21 13  GD2 44 -24 16  GD2 46 -52 13
GD2 48  28  9  GD2  2  15  9  GD2  4 -44 10
GD2  5  -2 12  GD2  9   8  7  GD2 11  12  7
GD2 13 -52  7  GD2 16  21  7  GD2 17  19 11
GD2 19   6 16  GD2 21  10 16  GD2 22 -15  6
GD2 25   4 15  GD2 27  -9  9  GD2 29  27 10
GD2 30   1 17  GD2 32  12  8  GD2 35  20  8
GD2 37 -32  8  GD2 38  15  8  GD2 41   5 14
GD2 43  35 13  GD2 45  28  9  GD2 47   6 15
")

hsv <- data.frame(scan(connection, list(vac="", pat=0, wks=0, x=0)))
hsv <- transform(hsv, status= (wks >0), wks = abs(wks))

fit1 <- survreg(Surv(wks, status) ~ 1, data=hsv, dist='exponential')
temp <- predict(fit1, type='quantile', p=.5, se=TRUE)

c(median= temp$fit[1], std= temp$se[1])
   median    std
24.32723  4.36930

--
The predict function gives the predicted median survival and standard deviation for each 
observation in the data set.  Since this was a mean only model all n of them are the same 
and I printed only the first.
For prediction they make the assumption that the std error for my future study will be the 
same as the std from this one, you want the future 95% CI to not include the value of 16, 
so the future mean will need to be at least 16 + 1.96* 4.369.

A nonparmetric version of the argument would be

> fit2 <- survfit(Surv(wks, status) ~ 1, data=hsv)
> print(fit2)
records   n.max n.start  events  median 0.95LCL 0.95UCL
      48      48      48      31      21      15      35

Then make the argument that in our future study, the 95% CI will stretch 6 units to the 
left of the median, just like it did here.  This argument is a bit more tenuous though. 
The exponential CI width depends on the total number of events and total follow-up time, 
and we can guess that the new study will be similar.  The Kaplan-Meier CI also depends on 
the spacing of the deaths, which is less likely to replicate.

Notes:
  1. Use summary(fit2)$table to extract the CI values.  In R the print functions don't 
allow you to "grab" what was printed, summary normally does.
  2. For the exponential we could work out the formula in closed form -- a good homework 
exercise for grad students perhaps but not an exciting way to spend my own afternoon.  An 
advantage of the above approach is that we can easily use a more realistic model like the 
weibull.
  3. I've never liked extracting out the "Surv(t,s)" part of a formula as a separate 
statement on another line.  If I ever need to read this code again, or even just the 
printout from the run, keeping it all together gives much better documentation.
  4. Future calculations for survival data, of any form, are always tenuous since they 
depend critically on the total number of events that will be in the future study.  We can 
legislate the total enrollment and follow-up time for that future study, but the number of 
events is never better than a guess.  Paraphrasing a motto found on the door of a well 
respected investigator I worked with 30 years ago (because I don't remember it exaclty):

   "The incidence of the condition under consideration and its subsequent death rate will 
both drop by 1/2 at the commencement of a study, and will not return to their former 
values until the study finishes or the PI retires."


Terry T.

---------------------------------------------------------------------------

On 07/10/2014 05:00 AM, r-help-request at r-project.org wrote:
> Hello All,
>
> I'm trying to figure out how to perform a survival analysis with an historical control. I've spent some time looking online and in my boooks but haven't found much showing how to do this. Was wondering if there is a R package that can do it, or if there are resources somewhere that show the actual steps one takes, or if some knowledgeable person might be willing to share some code.
>
> Here is a statement that describes the sort of analyis I'm being asked to do.
>
> A one-sample parametric test assuming an exponential form of survival was used to test the hypothesis that the treatment produces a median PFS no greater than the historical control PFS of 16 weeks.  A sample median PFS greater than 20.57 weeks would fall beyond the critical value associated with the null hypothesis, and would be considered statistically significant at alpha = .05, 1 tailed.
>
> My understanding is that the cutoff of 20.57 weeks was obtained using an online calculator that can be found at:
>
> http://www.swogstat.org/stat/public/one_survival.htm
>
> Thus far, I've been unable to determine what values were plugged into the calculator to get the cutoff.
>
> There's another calculator for a nonparamertric test that can be found at:
>
> http://www.swogstat.org/stat/public/one_nonparametric_survival.htm
>
> It would be nice to try doing this using both a parameteric and a non-parametric model.
>
> So my first question would be whether the approach outlined above is valid or if the analysis should be done some other way. If the basic idea is correct, is it relatively easy (for a Terry Therneau type genius) to implement the whole thing using R? The calculator is a great tool, but, if reasonable, it would be nice to be able to look at some code to see how the numbers actually get produced.