Hi Chris, Thanks for your reply. Certainly no need to apologize for a delayed response. Appreciate your taking the time to answer my question.
My concern was about the value of "16". My understanding is it is based on one expert's opinion. I wondered if this is in keeping with whatever the standard practice might be. I feel uncomfortable with the idea of using this as the basis but perhaps I am alone in feeling this way. My preference would be to see citations from previous studies to justify the number. A meta analysis of results from previous studies might be even better. Your code is very interesting. Especially the parametric part. Nice to get an actual p-value for the test. Thanks, Paul -------------------------------------------- On Mon, 7/21/14, Andrews, Chris <chri...@med.umich.edu> wrote: Subject: RE: [R] Survival Analysis with an Historical Control Cc: "r-help@r-project.org" <r-help@r-project.org> Received: Monday, July 21, 2014, 10:33 AM Hi Paul, Sorry for the delayed reply. I was away last week. I'm not clear what you want confirmed about your approach. (a) "20.57"- computing the rejection region of the analysis. The formulas implemented at the addresses you gave in your original post are from a reputable source - Lawless (1982) - at least that is claimed in the help file (http://www.swogstat.org/stat/Public/Help/survival1.html). It seems that another person derived 20.57 from some combination of input. I don't see how to back calculate what the input was from the information provided. Perhaps elsewhere in your study protocol there is discussion of accrual time, et al. that can help you. (b) "16" - the value of the parameter in the null hypothesis is very important as you noted in your response to Terry. But this is not really a statistical question. It may be derived from historical data, expert opinion, regulatory mandate, or some combination of these and other factors. Presumably this study was undertaken because 16 was an important number to somebody. (c) performing the one-sample survival analysis itself. This is what I did with the data you provided. # Non-parametric (km <- survfit(Surv(WKS, 1-CENS) ~ 1, data=hsv, type="kaplan-meier", conf.type="log", conf.int=0.9)) summary(km) # Compare to median survival = 16 # (Used 90% CI above to get 0.05 one sided test here) quantile(km, prob=0.5)$lower > 16 # Parametric (paraExp <- survreg(Surv(WKS, 1-CENS) ~ 1, data=hsv, dist="exponential")) summary(paraExp) # Compare to median survival = 16 # That is, compare to beta0 = log(16/log(2)) = 3.1391... # one sided test pnorm(c((coef(paraExp) - log(16/log(2))) / sqrt(vcov(paraExp))), lower.tail=FALSE) Chris ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.