Cox-regression:proportional hazard test by Schoenfeld residuals

[Anete Kaldal]

Hi!
I have performed cox-regression on a data set with both continuous and categorical variables. In the cox-regr. command I specified categorical values by adding i. before variable. Shoenfeld residuals were saved by following command stcox i.spiders age i.sex i.ascites albumin bilirubi i.edema1 choleste i.stage,schoenfeld(sc*) scaledsch(ssc*). When assessing proportionality I used estat phtest, detail. All these went well. When I attempted to make plots with Schoenfeld residuals, I used estat phtest, plot(var) command which went well for continuous variables such as age, albumin etc., but when assessing categorical variables such as sex estat phtest, plot(sex), I got error response that this variable is not found in model. Could you give an advice what is wrong?

[Carlo Lazzaro]

Anete:
welcome to this forum.
What if you type:

stphplot, by(sex)

[Anete Kaldal]

Thanks for quick response!
Well, the command you adviced worked well, just this was not the same graph I was supposed to make...
Here I have copied the commands from .do file:

*Cox-regression*
stcox i.spiders age i.sex i.ascites albumin bilirubi i.edema1 choleste i.stage

*proportional hazard test*
estat phtest, detail

*saving Schoenfeld residuals*
stcox i.spiders age i.sex i.ascites albumin bilirubi i.edema1 choleste i.stage,schoenfeld(sc*) scaledsch(ssc*)

estat phtest, plot(age)
estat phtest, plot(albumin)
estat phtest, plot(bilirubi)
estat phtest, plot(choleste)
estat phtest, plot(edema1)
estat phtest, plot(sex)

For all the continuous variables there was no problem but for all the categorical variables the error message was "sex not found in model" , r(198).

[Anete Kaldal]

I mean, the error message contained the variable that I tried to analyze, like "sex not found in model" or "edema1 not found in model".

[Carlo Lazzaro]

Anete:
yes, I was able to replicate your problem:

. use "http://www.stata-press.com/data/r16/kva.dta"
(Generator experiment)

. stcox load i.bearings,schoenfeld(sc*) scaledsch(ssc*)

         failure _d:  1 (meaning all fail)
   analysis time _t:  failtime

Iteration 0:   log likelihood = -20.274897
Iteration 1:   log likelihood = -10.515114
Iteration 2:   log likelihood = -8.8700259
Iteration 3:   log likelihood = -8.5915211
Iteration 4:   log likelihood = -8.5778991
Iteration 5:   log likelihood =  -8.577853
Refining estimates:
Iteration 0:   log likelihood =  -8.577853

Cox regression -- Breslow method for ties

No. of subjects =           12                  Number of obs    =          12
No. of failures =           12
Time at risk    =          896
                                                LR chi2(2)       =       23.39
Log likelihood  =    -8.577853                  Prob > chi2      =      0.0000

------------------------------------------------------------------------------
          _t | Haz. Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        load |    1.52647   .2188172     2.95   0.003     1.152576    2.021653
  1.bearings |   .0636433   .0746609    -2.35   0.019     .0063855    .6343223
------------------------------------------------------------------------------

. estat phtest, plot(load)

. estat phtest, plot(bearings)
bearings not found in model
r(198);

. estat phtest, plot(i.bearings)
0.bearings 1.bearings not found in model
r(198);


. estat phtest, detail

      Test of proportional-hazards assumption

      Time:  Time
      ----------------------------------------------------------------
                  |       rho            chi2       df       Prob>chi2
      ------------+---------------------------------------------------
      load        |     -0.01171         0.00        1         0.9707
      0b.bearings |            .            .        1             .
      1.bearings  |     -0.25239         0.87        1         0.3500
      ------------+---------------------------------------------------
      global test |                      1.92        2         0.3835
      ----------------------------------------------------------------

.

Actually, the issue creeps up with categorical variables only.
I cannot tell whether it depends on the -plot- option that calls -lowess- in turn or else.
However, what sounds weird is the error message stating that the categorical predictor was not found in the model.

[Anete Kaldal]

Yes, exactly. The continuous variables didn't cause any problem. I tried to define dichotomous variables as continuous, and it went well then, but in my data set there are some categorical variables with more than 2 values and then it causes wrong results.
Thank you indeed for your response, and please let me know if you come across solution for this.

[Paul Dickman]

Maybe I'm stating the obvious, but the following is probably what you want

estat phtest, plot(1.bearings)

[Carlo Lazzaro]

Anete:
Paul's code works for me (Thanks Paul).
However, the reason of the error message is still obscure to me.

[Paul Dickman]

However, the reason of the error message is still obscure to me.

My understanding is that the plot() option specifies that a plot be produced for a specific covariate in the model.

. stcox load i.bearings, schoenfeld(sc*) scaledsch(ssc*)
. estat phtest, plot(bearings)
bearings not found in model
r(198);

. estat phtest, plot(i.bearings)
0.bearings 1.bearings not found in model
r(198);

bearings is not a covariate in the model, so an error is returned.

i.bearings resolves to "0.bearings 1.bearings" which doesn't specify a single covariate in the model, so an error is returned.

Anete, if stage has 4 categories (coded 1, 2, 3, 4) then you could/should use:

estat phtest, plot(2.stage)
estat phtest, plot(3.stage)
estat phtest, plot(4.stage)

[Carlo Lazzaro]

Thanks, Paul.

[Anete Kaldal]

Thank you very much, Paul!
Now I understood the problem, and your advice helped to solve the problem.

Kind regards,
Anete