Generating multiple estat phtest plots for categorical predictors

[Michael Middletown]

Good morning all,

I am running estat phtest following stcox (Stata 15.1, updated)

I want to generate the plot of scaled Schoenfeld residuals for each categorical predictor as part of assumption testing (I have 20 categorical predictors against 500,000 observations). Categorical predictors have differing number of levels (from 2 to 10).

My initial hope had been to run a for loop for with estat phtest, plot(var) for each of the predictors contained in a var local . This is not possible however, because the plots only succeed when running estat phtest manually for each explicity declared level of the predictors (as per code example for deprivation).

estat phtest, plot(1.imd5)
graph export Schoenfeld_1imd5.svg, as(svg) replace
estat phtest, plot(2.imd5)
graph export Schoenfeld_2imd5.svg, as(svg) replace
estat phtest, plot(3.imd5)
graph export Schoenfeld_3imd5.svg, as(svg) replace
estat phtest, plot(4.imd5)
graph export Schoenfeld_3imd5, as(svg) replace

Keeping IMD5 as an example, none of the options below work (all generate error of 'not found in model'). I can understand why these options do not work (plot requires varname, rather than a list etc), but tried them anyway!

estat phtest, plot(imd5)
estat phtest, plot(i.imd5)
estat phtest, plot(1.imd5 2.imd5 3.imd5 4.imd5)

Is there are a way to request phtest, plot() to iterate through each level of predictor, for each variable, and generate a plot for each?

Or does it boil down to manually listing each level for each variable as the initial code chunk above?

Many thanks,

Michael

[Clyde Schechter]

levelsof imd5 if e(sample), local(levels)
foreach x of local levels {
    estat phtest, plot(`x'.imd5)
    graph export Schoenfeld_`x'imd5.svg, as(svg) replace
}

As an aside, I don't think it is necessary to both use the .svg filename extension and -as(svg)- in -graph export-. Either one will suffice.

[Michael Middletown]

Thank you very much for such a rapid solution Clyde, that works for a single predictor.

However, I am trying to apply your solution such that it can run through a number of predictors within a single step (rather than taking your block of code, and changing the predictor for 20 variables).

I have tried nesting your loop within another loop, as per code chunk below. However, this generates an error (r(101) factor-variable and time-series operators not allowed). I haven't been about to figure out why it is failing - but am assuming I've constructed my loop incorrectly.

Any further advice would be appreciated.

local predictors i.sex i.imd5 i.ethnic
foreach v of local predictors {
    levelsof `v' if e(sample), local(levels)
    foreach x of local levels {
        estat phtest, plot(`x'.`v')
        graph export graph_`x'`v'.svg, replace
        }
        }

[Clyde Schechter]

The problem is in the line -local predictors i.sex i.imd5 i.ethnic-. Get rid of the i.'s there. (Of course, don't get rid of i. prefixes in your -stcox- command.)

To see why, put yourself in Stata's shoes. With your definition of local predictors, when we reach the -levelsof- command on the first iteration of the outer loop, that command reads as -levelsof i.sex if e(sample), local(levels)-. That's what's throwing your error: -levelsof- requires a single variable name after it. i.sex is not a variable name. It's an expression that potentially represents many variables.

If you make that change, all goes well.

[Michael Middletown]

Hi Clyde, thank you!

Was focused on the loop and never clocked the error just above!