Dear Stata list

I would like to estimate a CI for the smoothed scaled Schoenfeld residuals after fitting a Cox regression model (to estimate the CI for the time varying hazard ratio).

A CI was provided in cox.zph in Splus and is also provided in some R packages (e.g. ggcoxzph) but this doesn't seem to be in the estat phtest options after stcox.

Would anyone be able to advise on how to calculate the CI manually after smoothing the residuals? Can this be done in a straightforward way using predictnl, or similar, as in the example code below (line marked ***)?

As the pattern of the time varying HR in my own dataset is complex, I've been using cubic splines for the smoothing. There are also multiple episodes per participant in my example (but not in the dataset used below).

thanks in advance for any advice

Matt Cairns

Code:
**************************************************************************************   
webuse drugtr, clear        
stcox drug,
estat phtest, plot(drug)
predict scaled*, scaledsch  // get scaled Schoenfeld residuals

// manually generate lowess of residuals to allow joint plotting below:
lowess scaled1 _t, bwidth(0.8) mean noweight generate(lowess)

// generate cubic spline for smoothing
mkspline rcs = _t if _d==1, cubic nknots(3)
        
// Estimate B(t) for drug, from the scaled schoenfeld residual        
regress scaled1  rcs*,    //
predictnl Bt = xb(), ci(Bt_LB Bt_UB)   // ***
                            
twoway scatter    scaled1     _t , mcolor(red) ///
   ||  connected  Bt    _t , sort msymbol(none) lcolor(black) ///            
   ||  connected  Bt_LB _t , sort msymbol(none) lcolor(black) lpattern(dash) ///
   ||  connected  Bt_UB _t , sort msymbol(none) lcolor(black) lpattern(dash) ///
   ||  connected  lowess _t, sort msymbol(none) lcolor(red) ///    
   legend(cols(1))  ///
   name(drug, replace)