0
$\begingroup$

My data have non-proportional hazards with clear separation. Should I handle it via stratified Cox regression or using separate Log-rank test within subsets?

I will use R only to illustrate.

I want to split the follow up into two periods, T<100 days and T>= 100 days. For the Cox analysis I can use the survSplit(), but survdiff doesn't work with it.

lung$ph.karno_cat <- ifelse(lung$ph.karno < 100, "A", "B")
lung.split <- survSplit(Surv(time, status) ~ ., data= lung, cut=c(180, 350), episode= "tgroup", id="id")

> survdiff(Surv(tstart, time, status) ~ ph.karno_cat:strata(tgroup), data=lung.split)
Error in survdiff(Surv(tstart, time, status) ~ ph.karno_cat:strata(tgroup),  : 
  Right censored data only

Having just one categorical covariate, can I just use the p-values from the coxph() per each stratum? They test whether the coefficient is non-zero and should correspond to the log-rank, isn't it?

> coef(summary(coxph(Surv(tstart, time, status) 
~ ph.karno_cat:strata(tgroup), data=lung.split) ))

                                            coef exp(coef)  se(coef)           z  Pr(>|z|)
ph.karno_catA:strata(tgroup)tgroup=1  0.82367716  2.278864 0.5170404  1.59306146 0.1111464
ph.karno_catB:strata(tgroup)tgroup=1          NA        NA 0.0000000          NA        NA
ph.karno_catA:strata(tgroup)tgroup=2  0.65219820  1.919756 0.4709701  1.38479747 0.1661144
ph.karno_catB:strata(tgroup)tgroup=2          NA        NA 0.0000000          NA        NA
ph.karno_catA:strata(tgroup)tgroup=3 -0.03294059  0.967596 0.3683481 -0.08942787 0.9287419
ph.karno_catB:strata(tgroup)tgroup=3          NA        NA 0.0000000          NA        NA

Or should I use the survdiff on filtered data?

> survdiff(Surv(time, status) ~ ph.karno_cat, data=lung, subset = time < 180)
Call:
survdiff(formula = Surv(time, status) ~ ph.karno_cat, data = lung, 
    subset = time < 180)

n=67, 1 observation deleted due to missingness.

                N Observed Expected (O-E)^2/E (O-E)^2/V
ph.karno_cat=A 63       57    59.17    0.0795      2.72
ph.karno_cat=B  4        4     1.83    2.5704      2.72

 Chisq= 2.7  on 1 degrees of freedom, p= 0.1 


> survdiff(Surv(time, status) ~ ph.karno_cat, data=lung, subset = time >= 180 & time < 350)
Call:
survdiff(formula = Surv(time, status) ~ ph.karno_cat, data = lung, 
    subset = time >= 180 & time < 350)

                N Observed Expected (O-E)^2/E (O-E)^2/V
ph.karno_cat=A 72       46     40.3      0.82      4.07
ph.karno_cat=B 13        5     10.7      3.07      4.07

 Chisq= 4.1  on 1 degrees of freedom, p= 0.04 

and

> survdiff(Surv(time, status) ~ ph.karno_cat, data=lung, subset = time >= 350)
Call:
survdiff(formula = Surv(time, status) ~ ph.karno_cat, data = lung, 
    subset = time >= 350)

                N Observed Expected (O-E)^2/E (O-E)^2/V
ph.karno_cat=A 63       43    43.08  0.000143  0.000839
ph.karno_cat=B 12        9     8.92  0.000689  0.000839

 Chisq= 0  on 1 degrees of freedom, p= 1 

The p-values are different, but quite close to the p-values for the Cox model. Which option would you prefer and why?

EDIT: I used the statistical package only to illustrate the problem. I am not asking for anyone to write code for me and I don't ask for datasets or debugging.

$\endgroup$

1 Answer 1

0
$\begingroup$

The way you did the separate survdiff() evaluations, based on subsetting from the original lung data set, means that each evaluation was only done for those with time values within each time-period subset. It ignores all those who might have been at risk during that time period but only have event/censoring times beyond that period. That's incorrect, as it's ignoring individuals with informative right-censored survival times (from the perspective of your first 2 time subsets).

That's why survSplit() defines each tgroup to include all at risk during the corresponding time period. (Also, survSplit() uses left-open/right-closed time limits rather than the left-closed/right-open time limits you used for subsetting the data.)

I suppose you could do survdiff() on each tgroup, but what does that gain you over the Cox model? With the Cox model you can incorporate the additional covariates that typically are needed to describe a data set adequately.

Finally, although you say "My data have non-proportional hazards with clear separation," that might not necessarily be the case. It's possible for a mis-specification of the functional form of a Cox model to masquerade as non-proportionality. Omission of critical predictors or interaction terms, or improper transformations of continuous predictors, can lead to apparent violations of proportional hazards that can be fixed by fixing the model itself. Unless there's a compelling reason to believe that there is a clear separation based on your knowledge of the subject matter rather than inspection of your results, consider that possibility before you go on to time stratification.

$\endgroup$
7
  • $\begingroup$ Part #1: Thank you. I have only one covariate, nothing else. In my real example (I couldn't share, it's secret, I have NDA), for 90% of the follow-up period the curves are identical, no difference (maybe just 1-2 events). In the remaining 10%, where I clearly see the time point, they rapidly diverge, with HR=12 (twelve). When I use logrank for the entire period, it finds nothing (p=0.7). Fleming-Harrington (0,1) detects it perfectly, but I wanted the HRs to be calculated too. So I decided to split the whole follow up into 2 pieces and calculate HRs and their log-ranks per time period. $\endgroup$
    – Aronovsky
    Commented May 7, 2022 at 16:46
  • $\begingroup$ Part #2. So I think the first method will be fine, and the p-values for the levels of the categorical covariate (the only I have) will correspond to the log-ranks per period, as I wanted, the same way it works for the ordinary Cox with single categorical covariate (the beta coefficient p-value = overall logrank). Am I right? Regarding the expectation - I did know they will diverge at a very late time point, but I couldn't anticipate where exactly. It's typical to have them equal during about 70-90% of the follow up. I could use the FH(0,1), but I want the piecewise HR and corresponding p-value $\endgroup$
    – Aronovsky
    Commented May 7, 2022 at 16:49
  • $\begingroup$ @Aronovsky if you use the survdiff method, then for the first time period you need to include all cases, with those extending into the second time period coded as censored at the end of the first time period. That's the way the survSplit correctly handles data. The last time period will be limited to those who got as far as the beginning of that time period; that's OK. $\endgroup$
    – EdM
    Commented May 7, 2022 at 17:23
  • $\begingroup$ @Aronovsky As this vignette says right after discussing the choice of p and q in FH: "Remember that test selection should be performed at the research design level! Not after looking in the dataset." Your p-values will not be reliable, as you looked at the results to decide where to place the break-time point and only then did the Cox model to get a coefficient. And, maybe, even to make your choice of (0,1) for FH. $\endgroup$
    – EdM
    Commented May 7, 2022 at 17:30
  • $\begingroup$ Then I have no option left. I have to do something. I decided to use the MaxCombo test, which is what the pharmaceutical industry pushed for, which uses FH(0,0), FH(0,1), FH(1,0) and FH(1,1) able to detect any kind of change - and it did well. But still I need the piecewise HR and logrank tests or at least CIs. I did foresee that there'll be late effects (so I'm allowed to use FH(0,1)), but I could NOT anticipate the exact day to split at. It's impossible in my case. Since I must complete the analysis, I will use the piecewise Cox with strata and take the p-values as log-rank ones. Thank you $\endgroup$
    – Aronovsky
    Commented May 7, 2022 at 17:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.