Multiple comparisons after Cox Model

[Roberto Vidri]

Hi all,

I'm seeking advice on the most appropriate way to perform multiple group comparison after running a Cox Model (code below). Using the pwcompare command provides p-values, but not HRs to go along. I also understand the Bonferroni corrections can be over conservative.

The main exposure (groups446) has 3 categories of patients (0, 1, 2). Some of the reviewers for this manuscript are asking for direct comparison of individual groups (ex. group 0 vs. group 1, 0 vs 2, 1 vs. 2).

What is the "recommended" or most acceptable approach?

Is it valid/acceptable to run the model as originally set up include two groups of interest at a time? For example adding "stock i.groups446 age......if groups446!=1"?

    stcox ib2.groups446 age i.sex i.race i.hispanic i.asaclass i.locheadneck i.radiographicstage i.natfl i.chemotherapyswitch i.neoadjuvantradiation i.adjuvantradiation i.veinresection i.arterialresection i.intraopbloodtransfusion lnr i.lymphovascularinvasion i.pathologiccompleteresponse i.year i.site, base

I would really appreciate some input, as I'm on a tight timeline with the reviewers.

Best,
Roberto

[Bruce Weaver]

Hello Roberto Vidri. I must start with a disclaimer: It is a long time since I last estimate a Cox regression model (other than some toy examples, maybe), and when I did, I think it was via SPSS rather than Stata. But I became interested in your post when I saw that your group variable has 3 levels. As you may know, Fisher's LSD properly controls the family-wise (FW) error rate when there are k = 3 groups. See Meier (2006), for example. Or see Statistical Methods for Psychology (by the late David Howell) or Serious Stats (by Thom Baguley).

Now before anyone points out that Fisher's LSD entails making all pair-wise comparisons among means (via modified t-tests) when the omnibus F-test (in an ANOVA model) is statistically significant, I know that. But I also know that I have read somewhere that the same logic and same two-step procedure can be applied in the context of other models. I.e., when one has k=3 groups, one must first examine the overall test (with df=2). Then, only if the overall test is statistically significant, one can carry out all pair-wise comparisons with no adjustment to the per contrast alpha. (If I remember where I read this, I will follow up, of course.)

So, if you buy that the logic of Fisher's LSD (with k=3) groups is applicable to your situation. then I think you could do something like the following. HTH.

sysuse cancer, clear
describe
* Generate a group variable with 3 levels.
generate group = mod(_n,3)
tabulate group

stcox i.group i.drug age, nolog
stcox, nohr

contrast group
* The overall test is significant.
* Now use -pwcompare- to get the pairwise comparisons.
pwcompare group, effects
* Store table in matrix pw
matrix pw = r(table_vs)'
* Save matrix to working dataset
svmat pw
generate hr = exp(pw1) // hazard ratio
generate lb = exp(pw5) // Lower bound of 95% CI
generate ub = exp(pw6) // Upper bound of 95% CI
list pw1 pw5 pw6 hr lb ub in 1/3
* The 3 HRs are 1 vs 0, 2 vs 0, and 2 vs 1.
* Notice that the first two results match exactly
* the coefficients from -stcox-.
stcox

Here is the last part of the output:

. list pw1 pw5 pw6 hr lb ub in 1/3

     +--------------------------------------------------------------------+
     |       pw1         pw5         pw6         hr         lb         ub |
     |--------------------------------------------------------------------|
  1. | -1.419171   -2.460257   -.3780854   .2419145    .085413    .685172 |
  2. | -.3685452   -1.209316    .4722261   .6917399   .2984012    1.60356 |
  3. |  1.050626    .0204999    2.080752    2.85944   1.020711   8.010488 |
     +--------------------------------------------------------------------+

. * The 3 HRs are 1 vs 0, 2 vs 0, and 2 vs 1.
. * Notice that the first two results match exactly
. * the coefficients from -stcox-.
. stcox

Cox regression with Breslow method for ties

No. of subjects =  48                                   Number of obs =     48
No. of failures =  31
Time at risk    = 744
                                                        LR chi2(5)    =  44.99
Log likelihood = -77.418701                             Prob > chi2   = 0.0000

------------------------------------------------------------------------------
          _t | Haz. ratio   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       group |
          1  |   .2419145   .1284991    -2.67   0.008      .085413     .685172
          2  |   .6917399   .2967376    -0.86   0.390     .2984012     1.60356
             |
        drug |
      Other  |   .1312363   .0709319    -3.76   0.000     .0454976    .3785469
         NA  |   .0315008   .0237897    -4.58   0.000     .0071695    .1384063
             |
         age |   1.151651   .0443722     3.66   0.000     1.067886    1.241987
------------------------------------------------------------------------------

[Roberto Vidri]

@Bruce Weaver - Thank you for your suggestion. I was not aware of LSD or other similar, less conservative methods, for postdoc comparisons - this is very helpful.

Best,
Roberto

[Simon Chang]

I found this very helpfull, but also discovered that there is an option to get the hazard ratios directly from pwcompare by adding the eform option which automatically gives the exp. estimates.

pwcompare group, effects eform