On Apr 15, 2014, at 6:32 AM, Therneau, Terry M., Ph.D. wrote:
> You can do statistical tests within a single model, for whether portions of
> it fit or do not fit. But one cannot take three separate fits and compare
> them. The program needs context to know how the three relate to one another.
> Say that "group" is your strata variable, trt the variable of interest, and
> x1, x2 are adjusters.
>
> fit <- coxph(Surv(time,status) ~ trt * strata(group) + x1 + x2, data=mydata)
>
> Will fit a model with a separate treatment coefficient for each of the
> groups, and a separate baseline hazard for each. One can now create a
> contrast that corresponds to your trend test, using vcov(fit) for the
> variance matrix and coef(fit) to retrieve the coefficients.
>
I have at the moment on my workspace a dataset of breast cancer cases extracted
from SEER that has a factor representing three grades of histology:
$Grade.abb. $ Grade.abb : Factor w/ 3 levels "Well","Moderate", "Poorly"
It would be reasonable to test whether the grading has a "trend" in its effect
when controlling for other factors (and it would be surprising to a medical
audience if there were no effect.). So I compare across models using
AgeStgSiz.NG.Rad as my "linear" trend" model (with one df for an `as.numeric`
version, AgeStgSizRad as my no-grade-included baseline, and AgeStgSiz.FG.Rad as
my full factor version:
> AgeStgSiz.NG.Rad <- coxph( Surv(Survmon,
> Vital.status.recode..study.cutoff.used.=="Dead") ~ AgeDx+
> SEER.historic.stage.A+ as.numeric(Size.cat) + as.numeric(Grade.abb) + Rgrp ,
> data=BrILR)
> AgeStgSizRad <- coxph( Surv(Survmon,
> Vital.status.recode..study.cutoff.used.=="Dead") ~ AgeDx+
> SEER.historic.stage.A+ as.numeric(Size.cat) + Rgrp , data=BrILR)
> AgeStgSiz.FG.Rad <- coxph( Surv(Survmon,
> Vital.status.recode..study.cutoff.used.=="Dead") ~ AgeDx+
> SEER.historic.stage.A+ as.numeric(Size.cat) + Grade.abb + Rgrp , data=BrILR)
> 2*diff(summary(AgeStgSizRad)[['loglik']])
[1] 5166.291
> 2*diff(summary(AgeStgSiz.NG.Rad)[['loglik']])
[1] 5888.492
> 2*diff(summary(AgeStgSiz.FG.Rad)[['loglik']])
[1] 5889.379
So there is strong evidence that adding grade to the existing covariates
improves the fit but that representing as separate factor values with one extra
degree of freedom may not be needed. When I add grade as a stratum I get a very
different 2*loglikelihood:
> AgeStgSiz.SG.Rad <- coxph( Surv(Survmon,
> Vital.status.recode..study.cutoff.used.=="Dead") ~ AgeDx+
> SEER.historic.stage.A+ as.numeric(Size.cat) + strata(Grade.abb) + Rgrp ,
> data=BrILR)
> 2*diff(summary(AgeStgSiz.SG.Rad)[['loglik']])
[1] 3980.908
> dput( vcov(AgeStgSiz.SG.Rad) )
structure(c(9.00241385282728e-06, -4.45446264398645e-07, 5.18927440846587e-07,
2.62020260612094e-07, 7.47434378232446e-07, -4.45446264398645e-07,
0.0046168537719431, 0.00445692601518848, -8.67833275051278e-07,
-3.68093395861629e-05, 5.18927440846587e-07, 0.00445692601518848,
0.00464685164887969, -1.61616621634903e-05, -3.77256742079467e-05,
2.62020260612094e-07, -8.67833275051278e-07, -1.61616621634903e-05,
7.84049821976807e-06, 1.8221575745622e-06, 7.47434378232446e-07,
-3.68093395861629e-05, -3.77256742079467e-05, 1.8221575745622e-06,
0.000313989310316303), .Dim = c(5L, 5L), .Dimnames = list(c("AgeDx",
"SEER.historic.stage.ALocalized", "SEER.historic.stage.ARegional",
"as.numeric(Size.cat)", "RgrpTRUE"), c("AgeDx",
"SEER.historic.stage.ALocalized",
"SEER.historic.stage.ARegional", "as.numeric(Size.cat)", "RgrpTRUE"
)))
> dput(coef(AgeStgSiz.SG.Rad))
structure(c(0.0393472050734995, 0.901971276489615, 1.46695753267995,
0.108860100649677, -0.273688779502084), .Names = c("AgeDx",
"SEER.historic.stage.ALocalized",
"SEER.historic.stage.ARegional", "as.numeric(Size.cat)", "RgrpTRUE"
))
I'm not particularly facile with contrast construction with var-covar matrices,
so hoping for a worked example. Also wondering of the cross-model comparisons
are invalid or less powerful?
> Terry T.
>
>
>
> On 04/15/2014 05:00 AM, r-help-requ...@r-project.org wrote:
>> Hello,
>>
>> I have the following problem. I stratified my patient cohort into three
>> ordered groups and performed multivariate adjusted Cox regression analysis
>> on each group separately. Now I would like to calculate a p for trend across
>> the hazard ratios that I got for the three groups. How can I do that if I
>> only have the HR and the confidence interval? For example I got the
>> following HRs for one endpoint:
>>
>> 1.09(0.68-1.74), 1.29(0.94-1.76) and 1.64(1.01-2.68).
>>
>> There is a trend but how do I calculate if it is significant?
>>
>> Best regards
>>
>> Marcus Kleber
>>
>
> ______________________________________________
> 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.
David Winsemius
Alameda, CA, USA
______________________________________________
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.
