Dear Terry Therneau:
Thank you for replying. Please forgive me for carrying on.
Here is the example I gave, now with some output shown:
require(survival)
data(pbc)
coxph(Surv(time,status)~edtrt,data=pbc)
--- OUTPUT
Call:
coxph(formula = Surv(time, status) ~ edtrt, data = pbc)
coef exp(coef) se(coef) z p
edtrt 2.22 9.22 0.241 9.22 0
Likelihood ratio test=61.1 on 1 df, p=5.44e-15 n= 418
--- END OUTPUT
m<-coxph(Surv(time,status)~edtrt+log(bili),data=pbc)
m # log(bili) is a strong confounder
--- OUTPUT
Call:
coxph(formula = Surv(time, status) ~ edtrt + log(bili), data = pbc)
coef exp(coef) se(coef) z p
edtrt 1.435 4.20 0.2495 5.75 8.9e-09
log(bili) 0.895 2.45 0.0807 11.10 0.0e+00
Likelihood ratio test=181 on 2 df, p=0 n= 418
--- END OUTPUT
plot(survfit(Surv(time,status)~edtrt,data=pbc))
lines(survexp(~edtrt+ratetable(edtrt=edtrt,bili=bili),data=pbc,ratetable
=m,cohort=TRUE),col="purple")
If I understand the PBC-dataset correctly, the 'edtrt'-variable is not
treatment, but edema score. Unlike in the survexp-example, I do not use
the 'trt' variable in my example at all. What I meant by "strong
confounding" was the change in the HR for edema, from 9.22 without
adjustment for log(bili) to 4.20 with adjustment. I have read your book
several times, but I left it at work, and Google Print doesn't show page
281, where I believe that "direct-adjusted survival" is mentioned. Also,
if I remember correctly, the entire chapter 10 concernes using rates
from one population to predict rates in another population. If a revised
version of your book is forthcoming, as I strongly hope, I would love to
see a discussion of how the Ederer method can be used to adjust for
confounding, as in the above example where I (am trying to) plot the
survival probabilities for the three edema groups, as they would have
been if they all had the same bilirubin distribution. I'm in over my
head here, so please forgive me if I'm overlooking something obvious.
Best regards,
Peter.
-----Oprindelig meddelelse-----
Fra: Terry Therneau [mailto:[EMAIL PROTECTED]
Sendt: 31. januar 2008 15:52
Til: Peter Jepsen
Cc: r-help@r-project.org
Emne: Re: Direct adjusted survival?
> The lines that I hoped to be the survival probabilities for each
edtrt-group
> adjusted for confounding by log(bili) are nearly identical to the
KM-lines,
> and they certainly don't appear adjusted for the very strong
confounding by
> log(bili). I'm not quite sure what they are, though.
Yes, survexp will fit direct adjusted curves (and also the Hakulinen
and
conditional methods). For your example, I would expect that the
ordinary
Kaplan-Meier curves for treatment 1 vs 2 should be almost identical to
the
adjusted curves for treatment 1 vs 2. The PBC data is from a randomized
trial,
the two treatment arms are (not surprisingly) very well balanced with
respect
to bilirubin values, and so adjusting for imbalance makes no real
change. This
is exactly what the survexp code that you gave does.
If you are expected the curves to change, then I guess I'm not sure
what you
mean by "strong confounding". Bilirubin is perhaps the most important
clinical
measure of severity for any of the cholestatic liver diseases, of which
PBC is
one; but being a strong predictor of mortality does not necessarily
imply
confounding.
Standard errors for the direct curve are daunting -- it is several
pages of
code in a Gail and Benichou (?) paper. I need to create an example for
doing
this with the bootstrap. One problem is the two sources of variation.
The
original Cox model's curves have variance of course, but do we consider
the
population of subjects being averaged over (for the DA curve) to be
fixed or
random?
For a long explanation of expected survival I would refer you to
chapter 10
of Therneau and Grambsch, "Modeling Survival Data". One of the more
confusing
aspects is that things get re-discovered and renamed, the "direct
adjusted
survival" curve for instance is Ederer's method (1961) brought forward
to a
Cox model. The ideas are not hard, but it does take a whole chapter.
Terry Therneau
______________________________________________
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.
