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.
Below are some sample survival data and code in case this proves helpful.
Thanks,
Paul
###################################
#### Example Data: GD2 Vaccine ####
###################################
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, CENS=ifelse(WKS < 1, 1, 0), WKS=abs(WKS))
head(hsv)
require("survival")
survObj <- Surv(hsv$WKS, hsv$CENS==0) ~ 1
km <- survfit(survObj, type=c("kaplan-meier"))
print(km)
paraExp <- survreg(survObj, dist="exponential")
print(paraExp)