Survival analysis and proportionality hazard

[Sarina Hamidizade]

I would appreciate if you could let me know whether it is correct to use parametric panel survival regression for my data. If so, could you please let me know how to test proportionality hazard assumption? My data set is as follows:

ID	TIME	EVENT	x1	x2	x3	x4	x5
1	1	0	1.281024	0.022934	0.874538	1.216367	0.06094
1	2	0	1.270012	0.006454	0.820673	1.004861	-0.01408
1	3	0	1.05289	-0.05938	0.922422	0.729264	0.020476
1	4	0	1.113115	-0.01523	0.858878	0.809738	0.076037
1	5	1	1.219644	-0.05868	0.887088	0.484342	0.009778
2	1	0	1.062264	0.107021	0.814602	0.835928	0.19996
2	2	0	1.055786	0.081916	0.879486	0.686728	0.142627
2	3	0	0.970588	0.076064	0.906775	0.809796	0.165915
2	4	0	1.058996	0.130203	0.818112	0.875989	0.234452
2	5	0	1.124525	0.147841	0.75871	1.079925	0.276444
2	6	0	1.599781	0.262461	0.546151	1.312741	0.369479
2	7	0	1.575608	0.262096	0.564481	1.156476	0.348624
2	8	0	1.544273	0.240911	0.590729	1.07697	0.325612
2	9	0	1.721708	0.215246	0.552291	0.841011	0.2935
2	10	0	1.723163	0.20863	0.533981	0.786512	0.293033
2	11	0	1.630677	0.186235	0.547719	0.728193	0.273577
2	12	0	2.172313	0.319455	0.441393	0.946985	0.427395
3	1	0	0.874395	-0.03468	0.793502	0.609515	-0.00263
3	2	1	0.825239	-0.14194	0.952213	0.57288	-0.01915

I used the following commands to do survival analysis:
. xtset id time, yearly
. stset time, failure(event==1)
. xtstreg x1 x2 x3 x4 x5, dist(exponential)
. estat ic

The estimated regression:

Random-effects exponential regression Number of obs = 1,554
Group variable: id Number of groups = 152

Obs per group:
min = 2
avg = 10.2
max = 12

Integration method: mvaghermite Integration pts. = 12

Wald chi2(5) = 84.85
Log likelihood = -209.83151 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
x1 | .7474386 .5523192 -0.39 0.694 .1756228 3.181048
x2 | .0269099 .0935812 -1.04 0.299 .0000295 24.54876
x3 | 814.9069 1289.601 4.24 0.000 36.64971 18119.46
x4 | 1.032383 .4342369 0.08 0.940 .4527012 2.354344
x5 | .0005696 .0019129 -2.22 0.026 7.89e-07 .4112869
_cons | .0000404 .0000722 -5.66 0.000 1.22e-06 .0013406
-------------+----------------------------------------------------------------
/sigma2_u | .4575335 .352393 .1011175 2.070235
------------------------------------------------------------------------------
LR test vs. exponential model: chibar2(01) = 2.69 Prob >= chibar2 = 0.0504


Akaike's information criterion and Bayesian information criterion

-----------------------------------------------------------------------------
Model | Obs ll(null) ll(model) df AIC BIC
-------------+---------------------------------------------------------------
. | 1,554 . -209.8315 7 433.663 471.1031
-----------------------------------------------------------------------------
Note: N=Obs used in calculating BIC; see [R] BIC note.

Best regards,

[Carlo Lazzaro]

Sarina:
the outcome of LR test of random-effects exponential regression vs exponential model favours the latter.
As far as your other (broad) questions are concerned, I woud point you to the following Stata .pdf manual entries: -streg-; -streg postestimation-; -stcox postestimation.

[Sarina Hamidizade]

Dear Mr. Lazzaro
First, I thank you very much for sharing your time and knowledge. Then, as the Stata pdf manual mentions the likelihood-ratio test compares the random-effects
model with a survival model with fixed effects only. In my case, the results support the fixed-effects model. Therefore, I used the following command but I am not sure if it is fixed-effects or not:

. streg x1 x2 x3 x4 x5, dist(exponential)

The output:

Exponential regression -- log relative-hazard form

No. of subjects = 1,554 Number of obs = 1,554
No. of failures = 50
Time at risk = 9490
LR chi2(5) = 202.24
Log likelihood = -124.91894 Prob > chi2 = 0.0000

------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
x1 | .568155 .3714036 -0.86 0.387 .157774 2.045965
x2 | .355751 .9217623 -0.40 0.690 .0022164 57.1001
x3 | 456.0766 608.5318 4.59 0.000 33.36572 6234.118
x4 | 1.029038 .3758595 0.08 0.938 .5029531 2.105405
x5 | .0002391 .0006645 -3.00 0.003 1.03e-06 .0554472
_cons | .0001065 .0001614 -6.04 0.000 5.46e-06 .0020764
------------------------------------------------------------------------------

Unfortunately, I couldn't find how to test PH even though I look at it's post estimation manual.

Besides, I tried cox regression as follows:
. stset time, failure(event==1)
. stcox x1 x2 x3 x4 x5

Cox regression -- Breslow method for ties

No. of subjects = 1,554 Number of obs = 1,554
No. of failures = 50
Time at risk = 9490
LR chi2(5) = 203.00
Log likelihood = -219.66919 Prob > chi2 = 0.0000

------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
x1 | .4790793 .3186303 -1.11 0.269 .1301017 1.764136
x2 | 1.263872 3.29019 0.09 0.928 .0076878 207.7801
x3 | 601.7359 842.9844 4.57 0.000 38.63135 9372.857
x4 | 1.025448 .3465663 0.07 0.941 .5287334 1.988798
x5 | .0000954 .0002607 -3.39 0.001 4.48e-07 .0202743
------------------------------------------------------------------------------
. estat phtest, detail

Test of proportional-hazards assumption

Time: Time
----------------------------------------------------------------
| rho chi2 df Prob>chi2
------------+---------------------------------------------------
x1 | -0.27573 4.04 1 0.0445
x2 | -0.17849 2.07 1 0.1506
x3 | -0.18528 1.89 1 0.1696
x4 | -0.08159 0.19 1 0.6588
x5 | 0.17416 2.00 1 0.1577
------------+---------------------------------------------------
global test | 7.23 5 0.2044
----------------------------------------------------------------
Best regards,

[Carlo Lazzaro]

Sarina:
- fixed effect usually entails a factor variable among predictors, like in the following toy-example (factor variable: -wbc3-):

. use http://www.stata-press.com/data/r14/leukemia
(Leukemia Remission Study)


. stset weeks, failure(relapse==1)

     failure event:  relapse == 1
obs. time interval:  (0, weeks]
 exit on or before:  failure

------------------------------------------------------------------------------
         42  total observations
          0  exclusions
------------------------------------------------------------------------------
         42  observations remaining, representing
         30  failures in single-record/single-failure data
        541  total analysis time at risk and under observation
                                                at risk from t =         0
                                     earliest observed entry t =         0
                                          last observed exit t =        35

. stcox treatment1 i.wbc3, basel

         failure _d:  relapse == 1
   analysis time _t:  weeks

Iteration 0:   log likelihood =  -93.98505
Iteration 1:   log likelihood = -80.258598
Iteration 2:   log likelihood = -80.054222
Iteration 3:   log likelihood =  -80.05387
Refining estimates:
Iteration 0:   log likelihood =  -80.05387

Cox regression -- Breslow method for ties

No. of subjects =           42                  Number of obs    =          42
No. of failures =           30
Time at risk    =          541
                                                LR chi2(2)       =       27.86
Log likelihood  =    -80.05387                  Prob > chi2      =      0.0000

------------------------------------------------------------------------------
          _t | Haz. Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
  treatment1 |   .2883715   .1241628    -2.89   0.004     .1240094    .6705794
             |
        wbc3 |
          0  |          1  (base)
          1  |   4.956401   2.311923     3.43   0.001     1.986652    12.36548
------------------------------------------------------------------------------

- estat-phtest- tells you that the proportional hazards assumption is not violated in your -stcox-.

[Sarina Hamidizade]

Dear Mr. Lazzaro
Greetings
I am so grateful for your detailed example. However, in this example wbc3cat is represented as normal (wb1), moderate (wb2) and high (wb3) but I don't have such a variable so how should I deal with this situation? Really, I am investigating the impact of some financial ratios on the bankruptcy of companies.
Besides, could you please let me know how to test PH if I want to use for example . xtstreg x1 x2 x3 x4 x5, dist(loglogistics) while Prob >= chibar2 = 0.0001?

Best regards,

[Carlo Lazzaro]

Sarina:
I've previously posted a toy-example to get you familiar with the code: obviously, your categorical variable may well differ (eg.: country);
-dist(loglogistic) is available in accelerated failure time (AFT) metric only; hence, no PH and, consequently, no PH test.

[Sarina Hamidizade]

Dear Mr. Lazzaro
I am so sorry for asking you a lot. Unfortunately, when I use

. stcox x1 x2 x3 x4 x5 i.id, basel

this result is reported:

Iteration 22:  log likelihood = -65.042584  (backed up)
flat region resulting in a missing likelihood

In fact, id represents 152 companies. Therefore, I found this link but I am not sure if it helps.
Furthermore, I can't test PH after

. streg x1 x2 x3 x4 x5, dist(exponential).

Thanks in advance.
Best regards,

[Carlo Lazzaro]

Sarina:
-your first code makes the model unable to converge; probably the sort of misspecification is -i.id- (the informativeness of whic I fail to get);
-you can compare the HRs of -stcox- (for which postestimation commands aimed at checking the PH feature are available, such as -estat phtest-) with those obtained from the parametric regression that you ran. If they're not that different, you're on the right track.

[Sarina Hamidizade]

Dear Mr. Lazzaro
Thanks a lot for your time and consideration. Then, since there is a problem with the first link which is refered to by the following link I decided to use mixed-effects survival regression but I didn't fill out the options related to the random-effects.
http://www.statalist.org/forums/foru...fects-in-stcox

. stset time, failure(event==1) scale(1)

     failure event:  event == 1
obs. time interval:  (0, time]
 exit on or before:  failure

------------------------------------------------------------------------------
       1554  total observations
          0  exclusions
------------------------------------------------------------------------------
       1554  observations remaining, representing
         50  failures in single-record/single-failure data
       9490  total analysis time at risk and under observation
                                                at risk from t =         0
                                     earliest observed entry t =         0
                                          last observed exit t =        12

. mestreg x1 x2 x3 x4 x5, dist(exponential)

         failure _d:  event == 1
   analysis time _t:  time

Iteration 0:   log likelihood = -70315.658  
Iteration 1:   log likelihood = -352.53502  
Iteration 2:   log likelihood = -224.77077  
Iteration 3:   log likelihood = -214.16464  
Iteration 4:   log likelihood = -211.19934  
Iteration 5:   log likelihood = -211.17739  
Iteration 6:   log likelihood = -211.17739  

Exponential regression                          Number of obs     =      1,554

                                                Wald chi2(5)      =     196.32
Log likelihood = -211.17739                     Prob > chi2       =     0.0000
------------------------------------------------------------------------------
          _t | Haz. Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
          x1 |   .5681548   .3714034    -0.86   0.387      .157774    2.045964
          x2 |    .355751   .9217622    -0.40   0.690     .0022164    57.10008
          x3 |   456.0763   608.5313     4.59   0.000     33.36571    6234.112
          x4 |   1.029038   .3758594     0.08   0.938      .502953    2.105404
          x5 |   .0002391   .0006645    -3.00   0.003     1.03e-06    .0554472
       _cons |   .0001065   .0001614    -6.04   0.000     5.46e-06    .0020765
------------------------------------------------------------------------------
. estat phtest
subcommand estat phtest is unrecognized

However, I am not sure if it is right to estimate fixed effects or not?
Best regards,

[Carlo Lazzaro]

Sarina:
admittedly, I do not follow you.
You've proposed different models in different posts. Those models are not, in general, exchangeable.
I would recommend you to skim through the literature of your research field and see what others did in the past when presented with the same research topic.

[Sarina Hamidizade]

Dear Mr. Lazzaro
Thanks a lot for your help.
Best regards,