PH violation occurs after changing vce() option in stcox

[Nazli Uludere Aragon]

Dear Statalisters working on survival analyses,

Hope you are all well.

I would appreciate your suggestions to help me understand what might be underlying the following:

1. I estimate stcox x1 x2 x3 i.x4, vce(robust) followed by estat phtest, detail. I see no violation of PH.
2. However, stcox x1 x2 x3 i.x4, vce(cluster clusterID) results in both x1 and x2 violating PH.
3. Then I run stcox x1 x2 x3 i.x4, tvc(x1 x2) texp(ln(_t)) vce(cluster clusterID). Coefficients under tvc (time interactions) are not significant.

My sample is approx. 17000 (farms) spread over 1600 clusters (counties).

I have also read the intuition behind frailty. There could be county level frailty; it is not implausible. When I estimate a frailty model Stata takes a long time to estimate. It has succeeded once, and had a significant theta.

I am also wondering if some of my PH violations are due to use of a large dataset.

I use Stata Now 19.5 / SE for Windows.

Thank you.

[Carlo Lazzaro]

Nazli:
as per FAQ, sharing what Stata gave you back can help interested readers enormously. Thanks.

[Nazli Uludere Aragon]

Thank you, Carlo. Of course. I appreciate it. This is what I describe in step 1:

stcox cov5 cov4 cov3 cov2 cov1 group2 nicc_mean_inv log_SRR_2017  c.log_SRR_2017#c.nicc_mean_inv avg_sh_cty_broad1or2 avgLR
> abanrate , vce(robust) nolog

        Failure _d: fail1_inst==1
  Analysis time _t: stime1_months
note: cov1 omitted because of collinearity.

Cox regression with Breslow method for ties

No. of subjects =  17,393                              Number of obs =  17,393
No. of failures =  11,363
Time at risk    = 870,582
                                                       Wald chi2(11) = 2605.18
Log pseudolikelihood = -108651.49                      Prob > chi2   =  0.0000

------------------------------------------------------------------------------------------------
                               |               Robust
                            _t | Haz. ratio   std. err.      z    P>|z|     [95% conf. interval]
-------------------------------+----------------------------------------------------------------
                          cov5 |    .847145   .0209392    -6.71   0.000     .8070832    .8891954
                          cov4 |   .8833787   .0402164    -2.72   0.006     .8079702     .965825
                          cov3 |   .8876805   .0278701    -3.79   0.000     .8347028    .9440206
                          cov2 |   .3373182   .0167405   -21.90   0.000     .3060527    .3717777
                          cov1 |          1  (omitted)
                        group2 |   .8878377   .0219666    -4.81   0.000     .8458111    .9319525
                 nicc_mean_inv |   1.088177   .0965977     0.95   0.341     .9144044    1.294974
                  log_SRR_2017 |   1.115361   .1252484     0.97   0.331     .8950153    1.389953
                               |
c.log_SRR_2017#c.nicc_mean_inv |   .9886333   .0192872    -0.59   0.558     .9515446    1.027168
                               |
          avg_sh_cty_broad1or2 |   1.613051   .0693366    11.12   0.000     1.482721    1.754837
                         avgLR |   .6711678   .0222421   -12.03   0.000     .6289596    .7162085
                      abanrate |   .9782676   .0019573   -10.98   0.000     .9744388    .9821114
------------------------------------------------------------------------------------------------

. estat phtest , detail

Test of proportional-hazards assumption

Time function: Analysis time
--------------------------------------------------------
             |        rho     chi2       df    Prob>chi2
-------------+------------------------------------------
        cov5 |    0.00483     0.12        1       0.7244
        cov4 |   -0.00051     0.00        1       0.9718
        cov3 |    0.00552     0.16        1       0.6903
        cov2 |    0.01719     2.84        1       0.0919
      o.cov1 |          .        .        1           .
      group2 |   -0.00203     0.02        1       0.8786
nicc_mean_~v |    0.02040     2.14        1       0.1435
log_SRR_2017 |    0.02488     3.23        1       0.0725
c.log_SRR_~v |   -0.02074     2.16        1       0.1413
avg_sh_cty~2 |   -0.01355     0.95        1       0.3292
       avgLR |   -0.01827     1.85        1       0.1737
    abanrate |    0.01502     1.28        1       0.2575
-------------+------------------------------------------
 Global test |               14.45       11       0.2091
--------------------------------------------------------
Note: Robust variance–covariance matrix used.

And this is what I get in step 2:

stcox cov5 cov4 cov3 cov2 cov1 group2 nicc_mean_inv log_SRR_2017  c.log_SRR_2017#c.nicc_mean_inv avg_sh_cty_broad1or2 avgLR 
> abanrate , vce(cluster geoid_set) nolog

        Failure _d: fail1_inst==1
  Analysis time _t: stime1_months
note: cov1 omitted because of collinearity.

Cox regression with Breslow method for ties

No. of subjects =  17,393                               Number of obs = 17,393
No. of failures =  11,363
Time at risk    = 870,582
                                                        Wald chi2(11) = 819.65
Log pseudolikelihood = -108651.49                       Prob > chi2   = 0.0000

                                            (Std. err. adjusted for 1,643 clusters in geoid_set)
------------------------------------------------------------------------------------------------
                               |               Robust
                            _t | Haz. ratio   std. err.      z    P>|z|     [95% conf. interval]
-------------------------------+----------------------------------------------------------------
                          cov5 |    .847145   .0358974    -3.91   0.000     .7796298     .920507
                          cov4 |   .8833787   .0520292    -2.11   0.035     .7870691    .9914732
                          cov3 |   .8876805   .0386582    -2.74   0.006     .8150554    .9667768
                          cov2 |   .3373182   .0243004   -15.09   0.000     .2928999    .3884726
                          cov1 |          1  (omitted)
                        group2 |   .8878377   .0332658    -3.18   0.001     .8249744    .9554912
                 nicc_mean_inv |   1.088177    .155555     0.59   0.554     .8222811    1.440055
                  log_SRR_2017 |   1.115361   .2016798     0.60   0.546     .7825298    1.589754
                               |
c.log_SRR_2017#c.nicc_mean_inv |   .9886333   .0307962    -0.37   0.714     .9300795    1.050873
                               |
          avg_sh_cty_broad1or2 |   1.613051   .1393731     5.53   0.000     1.361763     1.91071
                         avgLR |   .6711678   .0470844    -5.68   0.000     .5849474     .770097
                      abanrate |   .9782676   .0042931    -5.01   0.000     .9698894    .9867181
------------------------------------------------------------------------------------------------

. estat phtest , detail

Test of proportional-hazards assumption

Time function: Analysis time
--------------------------------------------------------
             |        rho     chi2       df    Prob>chi2
-------------+------------------------------------------
        cov5 |    0.00475     0.36        1       0.5470
        cov4 |    0.00443     0.17        1       0.6793
        cov3 |    0.00961     0.96        1       0.3284
        cov2 |    0.01854     7.13        1       0.0076
      o.cov1 |          .        .        1           .
      group2 |   -0.00915     1.14        1       0.2846
nicc_mean_~v |    0.01258     2.50        1       0.1136
log_SRR_2017 |    0.01876     5.76        1       0.0164
c.log_SRR_~v |   -0.01259     2.42        1       0.1198
avg_sh_cty~2 |   -0.01488     4.79        1       0.0286
       avgLR |   -0.01803     8.25        1       0.0041
    abanrate |    0.01403     5.41        1       0.0200
-------------+------------------------------------------
 Global test |               48.83       11       0.0000
--------------------------------------------------------
Note: Robust variance–covariance matrix used.

This is what happens in step 3. Not all time interactions are significant.

. stcox cov5 cov4 cov3 cov2 cov1 group2 nicc_mean_inv log_SRR_2017  c.log_SRR_2017#c.nicc_mean_inv avg_sh_cty_broad1or2 avgLR 
> abanrate , tvc( avgLR abanrate log_SRR_2017 avg_sh_cty_broad1or2  ) vce(cluster geoid_set) nolog

        Failure _d: fail1_inst==1
  Analysis time _t: stime1_months

note: cov1 omitted because of collinearity.

Cox regression with Breslow method for ties

No. of subjects =  17,393                               Number of obs = 17,393
No. of failures =  11,363
Time at risk    = 870,582
                                                        Wald chi2(15) = 843.71
Log pseudolikelihood = -108639.8                        Prob > chi2   = 0.0000

                                            (Std. err. adjusted for 1,643 clusters in geoid_set)
------------------------------------------------------------------------------------------------
                               |               Robust
                            _t | Haz. ratio   std. err.      z    P>|z|     [95% conf. interval]
-------------------------------+----------------------------------------------------------------
main                           |
                          cov5 |   .8450982   .0357445    -3.98   0.000     .7778655    .9181419
                          cov4 |   .8787719   .0515008    -2.21   0.027     .7834137    .9857373
                          cov3 |   .8842318   .0384633    -2.83   0.005     .8119693    .9629253
                          cov2 |   .3362626   .0242239   -15.13   0.000     .2919841     .387256
                          cov1 |          1  (omitted)
                        group2 |   .8886637   .0332612    -3.15   0.002     .8258067    .9563052
                 nicc_mean_inv |   1.079638   .1537484     0.54   0.591     .8166963    1.427235
                  log_SRR_2017 |   1.032511    .188727     0.18   0.861     .7216193    1.477342
                               |
c.log_SRR_2017#c.nicc_mean_inv |   .9902887   .0307763    -0.31   0.754     .9317686    1.052484
                               |
          avg_sh_cty_broad1or2 |   1.721229   .1466454     6.37   0.000     1.456525    2.034039
                         avgLR |   .6898414   .0529716    -4.84   0.000     .5934541    .8018837
                      abanrate |   .9762529   .0044416    -5.28   0.000     .9675862    .9849972
-------------------------------+----------------------------------------------------------------
tvc                            |
                         avgLR |   .9987572   .0019711    -0.63   0.529     .9949014    1.002628
                      abanrate |   1.000134   .0001003     1.33   0.183      .999937     1.00033
                  log_SRR_2017 |   1.007067   .0022995     3.08   0.002     1.002571    1.011584
          avg_sh_cty_broad1or2 |   .9943387   .0025936    -2.18   0.030     .9892682    .9994351
------------------------------------------------------------------------------------------------
Note: Variables in tvc equation interacted with _t.

Lastly, this is the shared frailty model where I specify option shared(geoid_set). Stata worked overnight to provide estimates for this version.

 stcox cov5 cov4 cov3 cov2 cov1 group2 nicc_mean_inv log_SRR_2017  c.log_SRR_2017#c.nicc_mean_inv avg_sh_cty_broad1or2 avgLR 
abanrate , tvc( avgLR abanrate log_SRR_2017 avg_sh_cty_broad1or2) shared(geoid_set) nolog

note: cov1 omitted because of collinearity.

        Failure _d: fail1_inst==1
  Analysis time _t: stime1_months

Cox regression with Breslow method for ties
Gamma shared frailty                                Number of obs     = 17,393
Group variable: geoid_set                           Number of groups  =  1,643
                                                    Obs per group:   
No. of subjects =  17,393                                         min =      1
No. of failures =  11,363                                         avg =     11
Time at risk    = 870,582                                         max =     96
                                                    Wald chi2(15)     = 918.62
Log likelihood = -108547.32                         Prob > chi2       = 0.0000

------------------------------------------------------------------------------------------------
                            _t | Haz. ratio   Std. err.      z    P>|z|     [95% conf. interval]
-------------------------------+----------------------------------------------------------------
main                           |
                          cov5 |    .873843     .03371    -3.50   0.000     .8102086    .9424752
                          cov4 |     .88606   .0638291    -1.68   0.093     .7693875    1.020425
                          cov3 |   .8901085   .0433537    -2.39   0.017     .8090665    .9792682
                          cov2 |   .3641753   .0209824   -17.53   0.000     .3252877     .407712
                          cov1 |          1  (omitted)
                        group2 |   .9030463   .0330724    -2.78   0.005     .8404974    .9702501
                 nicc_mean_inv |   1.099948   .1710667     0.61   0.540     .8109446    1.491946
                  log_SRR_2017 |   1.040903   .2048026     0.20   0.839     .7078353    1.530693
                               |
c.log_SRR_2017#c.nicc_mean_inv |   .9921902   .0344426    -0.23   0.821     .9269292    1.062046
                               |
          avg_sh_cty_broad1or2 |   1.642621   .1439275     5.66   0.000     1.383421    1.950384
                         avgLR |   .7223323   .0475242    -4.94   0.000     .6349423    .8217502
                      abanrate |   .9704623   .0038219    -7.61   0.000     .9630002    .9779821
-------------------------------+----------------------------------------------------------------
tvc                            |
                         avgLR |  -.0013397   .0016238    -0.83   0.409    -.0045223    .0018428
                      abanrate |   .0001455   .0000877     1.66   0.097    -.0000264    .0003173
                  log_SRR_2017 |   .0055575   .0017122     3.25   0.001     .0022016    .0089134
          avg_sh_cty_broad1or2 |  -.0057603   .0022745    -2.53   0.011    -.0102182   -.0013024
-------------------------------+----------------------------------------------------------------
                         theta |   .0791964   .0096021
------------------------------------------------------------------------------------------------
LR test of theta=0: chibar2(01) = 184.95               Prob >= chibar2 = 0.000

Note: Variables in tvc equation interacted with _t.
Note: Standard errors of hazard ratios are conditional on theta.

[Carlo Lazzaro]

Nazli:
I woul stick with step 2 code removing statistically not significant interactions.

[Nazli Uludere Aragon]

Thank you very much for the advice! To confirm: you mean the version with tvc() option but removing the not significant interactions? Nazli

[Carlo Lazzaro]

Nazil:
no, I meant:

 stcox cov5 cov4 cov3 cov2 cov1 group2 nicc_mean_inv log_SRR_2017 c.log_SRR_2017#c.nicc_mean_inv avg_sh_cty_broad1or2 avgLR abanrate , vce(cluster geoid_set) nolog

Sorry for the confusion.

[Nazli Uludere Aragon]

Thanks a lot!