Improving fitting and predictions for flexible parametric survival models

Introduction

The first article in The Stata Journal introduced stpm

I wrote stpm2 in 2007.

Here I will introduce stpm3

Flexible Parametric Survival Models

Flexible parametric survival models are used with time-to-event (survival) outcomes
They are “flexible” in that they use spline functions to model the effect of time
Easy to relax proportionality assumptions through interactions between covariates and the effect of time.

Choice of Scale

Log Cumulative Hazard (stpm2)

\[ \ln[H(t|\mathbf{x}_i)] = \ln[-\ln(S(t|\mathbf{x}_i))] = \eta_i(t) = s\left(\ln(t)|\boldsymbol{\gamma}, \mathbf{k}_{0}\right) + \mathbf{x}_i \boldsymbol{\beta} \]

Log odds scale (stpm2)

\[ \ln\left[\frac{1-S(t|\mathbf{x}_i)}{S(t|\mathbf{x}_i)}\right] = \eta_i(t) = s\left(\ln(t)|\boldsymbol{\gamma}, \mathbf{k}_{0}\right) + \mathbf{x}_i \boldsymbol{\beta} \]

Log hazard scale (strcs) ln(time)

\[ \ln\left[h(t|\mathbf{x}_i)\right] = \eta_i(t) = s\left(\ln(t)|\boldsymbol{\gamma}, \mathbf{k}_{0}\right) + \mathbf{x}_i \boldsymbol{\beta} \]

`stcox` vs `stpm2`

stcox and stpm2

. stcox hormon, nolog noshow 
------------------------------------------------------------------------------
          _t | Haz. ratio   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
      hormon |   1.540262    .132659     5.02   0.000     1.301016    1.823503
------------------------------------------------------------------------------

. stpm2 hormon, scale(hazard) df(4) eform nolog
------------------------------------------------------------------------------
             |     exp(b)   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
xb           |
      hormon |   1.540779   .1326967     5.02   0.000     1.301464    1.824099
       _rcs1 |    2.50398    .069006    33.31   0.000     2.372319    2.642949
       _rcs2 |   1.198509   .0330973     6.56   0.000     1.135364    1.265166
       _rcs3 |   1.018274   .0145595     1.27   0.205     .9901346    1.047214
       _rcs4 |   .9961938   .0067963    -0.56   0.576     .9829618    1.009604
       _cons |   .2935573   .0097629   -36.85   0.000     .2750327    .3133296
------------------------------------------------------------------------------. stcox hormon, nolog noshow 
------------------------------------------------------------------------------
          _t | Haz. ratio   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
      hormon |   1.540262    .132659     5.02   0.000     1.301016    1.823503
------------------------------------------------------------------------------

. stpm2 hormon, scale(hazard) df(4) eform nolog
------------------------------------------------------------------------------
             |     exp(b)   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
xb           |
      hormon |   1.540779   .1326967     5.02   0.000     1.301464    1.824099
       _rcs1 |    2.50398    .069006    33.31   0.000     2.372319    2.642949
       _rcs2 |   1.198509   .0330973     6.56   0.000     1.135364    1.265166
       _rcs3 |   1.018274   .0145595     1.27   0.205     .9901346    1.047214
       _rcs4 |   .9961938   .0067963    -0.56   0.576     .9829618    1.009604
       _cons |   .2935573   .0097629   -36.85   0.000     .2750327    .3133296
------------------------------------------------------------------------------. stcox hormon, nolog noshow 
------------------------------------------------------------------------------
          _t | Haz. ratio   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
      hormon |   1.540262    .132659     5.02   0.000     1.301016    1.823503
------------------------------------------------------------------------------

. stpm2 hormon, scale(hazard) df(4) eform nolog
------------------------------------------------------------------------------
             |     exp(b)   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
xb           |
      hormon |   1.540779   .1326967     5.02   0.000     1.301464    1.824099
       _rcs1 |    2.50398    .069006    33.31   0.000     2.372319    2.642949
       _rcs2 |   1.198509   .0330973     6.56   0.000     1.135364    1.265166
       _rcs3 |   1.018274   .0145595     1.27   0.205     .9901346    1.047214
       _rcs4 |   .9961938   .0067963    -0.56   0.576     .9829618    1.009604
       _cons |   .2935573   .0097629   -36.85   0.000     .2750327    .3133296
------------------------------------------------------------------------------

Predictions

Why a new command?

stpm2 started before Stata factor variables
Use better basis functions for splines
Make (conditional) predictions easier
Use frames for predictions
Include splines on log hazard scale
Include functional forms of covariates in linear predictor
Make marginal/standardized predictions much easier.
- This is the main reason

Change of spline basis functions

stpm2 uses restricted cubic splines to model the effect of follow-up time.
- data dependent - problems with out of sample prediction.
New gensplines command used by stpm3
- Natural splines basis functions by default.
- Also possible to use B-splines

Restricted Cubic Splines in `stpm2`

Natural Splines in `stpm3`

Beyond the last knot

Factor Variables

Main effects

stpm3 i.hormon i.size age,           /// main effects
      scale(lncumhazard) df(3)       //  model scale & dfstpm3 i.hormon i.size age,           /// main effects
      scale(lncumhazard) df(3)       //  model scale & dfstpm3 i.hormon i.size age,           /// main effects
      scale(lncumhazard) df(3)       //  model scale & df

Interactions

stpm3 i.hormon##ib2.size age         /// Interactions
      scale(lncumhazard) df(3)       //  model scale & dfstpm3 i.hormon##ib2.size age         /// Interactions
      scale(lncumhazard) df(3)       //  model scale & dfstpm3 i.hormon##ib2.size age         /// Interactions
      scale(lncumhazard) df(3)       //  model scale & df

...with time-dependent effects

stpm3 i.hormon##i.size age,          /// Interactions
      tvc(i.hormon i.size) dftvc(2)  /// Interactions with time    
      scale(lncumhazard) df(3)       ///  model scale & df
      baselevelsstpm3 i.hormon##i.size age,          /// Interactions
      tvc(i.hormon i.size) dftvc(2)  /// Interactions with time    
      scale(lncumhazard) df(3)       ///  model scale & df
      baselevelsstpm3 i.hormon##i.size age,          /// Interactions
      tvc(i.hormon i.size) dftvc(2)  /// Interactions with time    
      scale(lncumhazard) df(3)       ///  model scale & df
      baselevels

Makes use of predict and standsurv much much easier.

Factor Variables Example


. stpm3 i.hormon##i.size age,          /// Interactions
>       tvc(i.hormon i.size) dftvc(2)  /// Interactions with time    
>       scale(lncumhazard) df(3)       ///  model scale & df
>       baselevels

Iteration 0:   log likelihood = -2909.4852  
Iteration 1:   log likelihood = -2902.0087  
Iteration 2:   log likelihood =  -2901.908  
Iteration 3:   log likelihood =  -2901.908  

                                                        Number of obs =  2,982
                                                        Wald chi2(6)  = 111.41
Log likelihood = -2901.908                              Prob > chi2   = 0.0000

-----------------------------------------------------------------------------------
                  | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
------------------+----------------------------------------------------------------
xb                |
           hormon |
              no  |          0  (base)
             yes  |    .405057    .239317     1.69   0.091    -.0639957    .8741098
                  |
             size |
         <=20 mm  |          0  (base)
       >20-50mmm  |   .4471447   .0961388     4.65   0.000     .2587162    .6355733
          >50 mm  |    1.02456   .1419684     7.22   0.000     .7463075    1.302813
                  |
      hormon#size |
   yes#>20-50mmm  |   -.097008   .2231739    -0.43   0.664    -.5344208    .3404048
      yes#>50 mm  |  -.0525507   .2586502    -0.20   0.839    -.5594958    .4543944
                  |
              age |   .0139055   .0022484     6.18   0.000     .0094988    .0183123
------------------+----------------------------------------------------------------
time              |
             _ns1 |  -28.61436   2.485529   -11.51   0.000    -33.48591   -23.74282
             _ns2 |    9.07063   1.222754     7.42   0.000     6.674077    11.46718
             _ns3 |  -1.746533   .1441763   -12.11   0.000    -2.029113   -1.463953
                  |
hormon#c._ns_tvc1 |
             yes  |   .7800628   1.573629     0.50   0.620    -2.304194     3.86432
                  |
hormon#c._ns_tvc2 |
             yes  |  -.8465712   .6738948    -1.26   0.209    -2.167381    .4742383
                  |
  size#c._ns_tvc1 |
       >20-50mmm  |  -1.370193   1.381062    -0.99   0.321    -4.077025     1.33664
          >50 mm  |   .8596402   1.567863     0.55   0.583    -2.213316    3.932596
                  |
  size#c._ns_tvc2 |
       >20-50mmm  |   1.190279    .420111     2.83   0.005     .3668762    2.013681
          >50 mm  |   1.040614   .5661812     1.84   0.066    -.0690804    2.150309
                  |
            _cons |  -.9955528   .1441542    -6.91   0.000     -1.27809   -.7130158
-----------------------------------------------------------------------------------
Warning: This is a test version of stpm3

Extended Functions

We rarely assume the effect of a covariate is linear.
Use splines, fractional polynomials etc.
May have interactions with these non-linear effects

This makes predictions complicated as you need to use the appropriate values of the spline basis functions.

extended functions - nonlinear function in a varlist.

Makes use of predict and standsurv much much easier.

Extended Functions Example 1

Natural splines - @ns()

stpm3 i.hormon i.size @ns(age, df(3)), /// natural spline for age
      scale(lncumhazard) df(3)         //  model scale & dfstpm3 i.hormon i.size @ns(age, df(3)), /// natural spline for age
      scale(lncumhazard) df(3)         //  model scale & df


. stpm3 i.hormon i.size @ns(age, df(3)), /// natural spline for age
>       scale(lncumhazard) df(3)         //  model scale & df

Iteration 0:   log likelihood = -2895.6086  
Iteration 1:   log likelihood = -2890.8231  
Iteration 2:   log likelihood = -2890.7951  
Iteration 3:   log likelihood = -2890.7951  

                                                        Number of obs =  2,982
                                                        Wald chi2(6)  = 355.97
Log likelihood = -2890.7951                             Prob > chi2   = 0.0000

------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
xb           |
      hormon |
        yes  |   .1882895    .087913     2.14   0.032     .0159832    .3605957
             |
        size |
  >20-50mmm  |   .6183297   .0634388     9.75   0.000      .493992    .7426674
     >50 mm  |   1.169518   .0871101    13.43   0.000     .9987854    1.340251
             |
 _ns_f1_age1 |  -2.334033   .8595762    -2.72   0.007    -4.018772   -.6492949
 _ns_f1_age2 |   -.926103   .3621613    -2.56   0.011    -1.635926   -.2162799
 _ns_f1_age3 |  -1.546971    .382083    -4.05   0.000     -2.29584   -.7981023
-------------+----------------------------------------------------------------
time         |
        _ns1 |  -29.15713    1.63076   -17.88   0.000    -32.35336    -25.9609
        _ns2 |   9.752766   .8468891    11.52   0.000     8.092894    11.41264
        _ns3 |  -1.494942   .0987252   -15.14   0.000     -1.68844   -1.301444
       _cons |   .7863726   .1988728     3.95   0.000      .396589    1.176156
------------------------------------------------------------------------------
Extended functions
 (1) @ns(age, df(3))
Warning: This is a test version of stpm3

Extended Functions Example 2

Use with interactions

stpm3 i.size                               ///
      i.hormon##@bs(age, df(3) degree(2)), /// interaction with bs function
      scale(lncumhazard) df(3)             //  model scale & df stpm3 i.size                               ///
      i.hormon##@bs(age, df(3) degree(2)), /// interaction with bs function
      scale(lncumhazard) df(3)             //  model scale & df


. stpm3 i.size                               ///
>       i.hormon##@bs(age, df(3) degree(2)), /// interaction with bs function
>       scale(lncumhazard) df(3)             //  model scale & df 

Iteration 0:   log likelihood = -2892.5789  
Iteration 1:   log likelihood = -2887.8413  
Iteration 2:   log likelihood = -2887.8145  
Iteration 3:   log likelihood = -2887.8145  

                                                        Number of obs =  2,982
                                                        Wald chi2(9)  = 362.88
Log likelihood = -2887.8145                             Prob > chi2   = 0.0000

--------------------------------------------------------------------------------------
                     | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
---------------------+----------------------------------------------------------------
xb                   |
                size |
          >20-50mmm  |   .6178289   .0634708     9.73   0.000     .4934283    .7422294
             >50 mm  |   1.166116   .0872509    13.37   0.000     .9951074    1.337125
                     |
              hormon |
                yes  |  -2.043709   1.342542    -1.52   0.128    -4.675042    .5876238
         _bs_f1_age1 |  -.9448818   .2968562    -3.18   0.001    -1.526709   -.3630544
         _bs_f1_age2 |  -.4385819    .198683    -2.21   0.027    -.8279935   -.0491703
         _bs_f1_age3 |   .6777717   .3473155     1.95   0.051    -.0029542    1.358498
                     |
hormon#c._bs_f1_age1 |
                yes  |   3.077341   1.590321     1.94   0.053    -.0396321    6.194313
                     |
hormon#c._bs_f1_age2 |
                yes  |   1.723866    1.26789     1.36   0.174     -.761153    4.208885
                     |
hormon#c._bs_f1_age3 |
                yes  |   2.533026   1.539937     1.64   0.100    -.4851953    5.551248
---------------------+----------------------------------------------------------------
time                 |
                _ns1 |  -29.19083   1.632273   -17.88   0.000    -32.39003   -25.99164
                _ns2 |   9.761409   .8475986    11.52   0.000     8.100146    11.42267
                _ns3 |  -1.502955   .0992601   -15.14   0.000    -1.697501   -1.308409
               _cons |   .2208277   .2201609     1.00   0.316    -.2106796    .6523351
--------------------------------------------------------------------------------------
Extended functions
 (1) @bs(age, df(3) degree(2))
Warning: This is a test version of stpm3

Extended Functions Example 3

stpm3 i.size                                 /// factor variable for size
      i.hormon##@bs(age, df(3) degree(2))    /// interaction bs function
      @poly(pr,degree(2)),                   /// quadratic for pr
      tvc(i.hormon @ns(age, df(3))) dftvc(2) /// different function of age    
      scale(lncumhazard) df(3)               //  model scale & df stpm3 i.size                                 /// factor variable for size
      i.hormon##@bs(age, df(3) degree(2))    /// interaction bs function
      @poly(pr,degree(2)),                   /// quadratic for pr
      tvc(i.hormon @ns(age, df(3))) dftvc(2) /// different function of age    
      scale(lncumhazard) df(3)               //  model scale & df stpm3 i.size                                 /// factor variable for size
      i.hormon##@bs(age, df(3) degree(2))    /// interaction bs function
      @poly(pr,degree(2)),                   /// quadratic for pr
      tvc(i.hormon @ns(age, df(3))) dftvc(2) /// different function of age    
      scale(lncumhazard) df(3)               //  model scale & df stpm3 i.size                                 /// factor variable for size
      i.hormon##@bs(age, df(3) degree(2))    /// interaction bs function
      @poly(pr,degree(2)),                   /// quadratic for pr
      tvc(i.hormon @ns(age, df(3))) dftvc(2) /// different function of age    
      scale(lncumhazard) df(3)               //  model scale & df stpm3 i.size                                 /// factor variable for size
      i.hormon##@bs(age, df(3) degree(2))    /// interaction bs function
      @poly(pr,degree(2)),                   /// quadratic for pr
      tvc(i.hormon @ns(age, df(3))) dftvc(2) /// different function of age    
      scale(lncumhazard) df(3)               //  model scale & df stpm3 i.size                                 /// factor variable for size
      i.hormon##@bs(age, df(3) degree(2))    /// interaction bs function
      @poly(pr,degree(2)),                   /// quadratic for pr
      tvc(i.hormon @ns(age, df(3))) dftvc(2) /// different function of age    
      scale(lncumhazard) df(3)               //  model scale & df


. stpm3 i.size                                 /// factor variable for size
>       i.hormon##@bs(age, df(3) degree(2))    /// interaction with bs function
>       @poly(pr,degree(2)),                   /// quadratic for pr
>       tvc(i.hormon @ns(age, df(3))) dftvc(2) /// different function of age    
>       scale(lncumhazard) df(3)               //  model scale & df 

Iteration 0:   log likelihood = -2872.9093  
Iteration 1:   log likelihood = -2852.5136  
Iteration 2:   log likelihood = -2851.7415  
Iteration 3:   log likelihood = -2851.7382  
Iteration 4:   log likelihood = -2851.7382  

                                                        Number of obs =  2,982
                                                        Wald chi2(11) = 346.49
Log likelihood = -2851.7382                             Prob > chi2   = 0.0000

------------------------------------------------------------------------------------------
                         | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------------------+----------------------------------------------------------------
xb                       |
                    size |
              >20-50mmm  |   .6002304   .0635988     9.44   0.000      .475579    .7248819
                 >50 mm  |   1.137697   .0872257    13.04   0.000     .9667376    1.308656
                         |
                  hormon |
                    yes  |    -1.7906   1.361312    -1.32   0.188    -4.458724    .8775227
             _bs_f1_age1 |  -.6016323   .4335225    -1.39   0.165    -1.451321    .2480562
             _bs_f1_age2 |   .0686366   .2907389     0.24   0.813    -.5012011    .6384743
             _bs_f1_age3 |   2.223462   .5670446     3.92   0.000     1.112075    3.334849
                         |
    hormon#c._bs_f1_age1 |
                    yes  |   2.803594   1.623107     1.73   0.084    -.3776377    5.984826
                         |
    hormon#c._bs_f1_age2 |
                    yes  |   1.525465   1.291585     1.18   0.238    -1.005995    4.056925
                         |
    hormon#c._bs_f1_age3 |
                    yes  |   2.479399   1.583994     1.57   0.118    -.6251728     5.58397
                         |
            _poly_f1_pr1 |  -.0007745   .0001618    -4.79   0.000    -.0010916   -.0004573
            _poly_f1_pr2 |   1.88e-07   6.10e-08     3.09   0.002     6.87e-08    3.08e-07
-------------------------+----------------------------------------------------------------
time                     |
                    _ns1 |  -40.88186   9.448295    -4.33   0.000    -59.40018   -22.36354
                    _ns2 |    16.3311   4.386648     3.72   0.000     7.733432    24.92878
                    _ns3 |  -.7436135   .4878034    -1.52   0.127    -1.699691    .2124636
                         |
       hormon#c._ns_tvc1 |
                    yes  |  -.9720923   1.567377    -0.62   0.535    -4.044095    2.099911
                         |
       hormon#c._ns_tvc2 |
                    yes  |  -.0983894   .6908374    -0.14   0.887    -1.452406    1.255627
                         |
c._bs_f1_age1#c._ns_tvc1 |   -.020625   6.570365    -0.00   0.997     -12.8983    12.85705
                         |
c._bs_f1_age1#c._ns_tvc2 |  -.6998829   1.753699    -0.40   0.690    -4.137069    2.737303
                         |
c._bs_f1_age2#c._ns_tvc1 |   9.092159   4.472742     2.03   0.042     .3257446    17.85857
                         |
c._bs_f1_age2#c._ns_tvc2 |  -3.010989   1.192695    -2.52   0.012    -5.348628   -.6733506
                         |
c._bs_f1_age3#c._ns_tvc1 |   10.17703   6.297045     1.62   0.106    -2.164951    22.51901
                         |
c._bs_f1_age3#c._ns_tvc2 |  -6.442536   2.195179    -2.93   0.003    -10.74501   -2.140064
                         |
                   _cons |  -.1301986    .312701    -0.42   0.677    -.7430812     .482684
------------------------------------------------------------------------------------------
Extended functions
 (1) @bs(age, df(3) degree(2))
 (2) @ns(age, df(3))
 (3) @poly(pr, degree(2))
Warning: This is a test version of stpm3

Current Extended Functions

`@bs()`	B-splines
`@fp()`	fractional polynomials
`@ns()`	natural cubic splines
`@poly()`	polynomials
`@rcs()`	restricted cubic splines

Other functions?

Makes use of predict and standsurv much much easier.

Predict to Frames

Save predictions in current frame
Save predictions to a new frame
Merge predictions with existing frame

Type of predictions

Different types of predictions
- Predict at observed values of covariates
- Predict at specific values of covariates
- Take average of predictions (marginal/standardized effects)
Difference choices for time
- Predict at _t
- Predict at single time point for all subjects
- Predict at user specified time values

Standard Predictions

Predict at observed values of covariates (and _t)

predict   xb, xb ci
predict    S, survival ci
predict    h, hazard ci

By default this type of prediction this will be saved in current frame.

Predictions

Model with factor variables and extended function

stpm3 i.hormon i.size @ns(age, df(3)), /// natural spline for age
      scale(lncumhazard) df(3)         //  model scale & dfstpm3 i.hormon i.size @ns(age, df(3)), /// natural spline for age
      scale(lncumhazard) df(3)         //  model scale & df

Predictions

predict S_age50_h0_s2 S_age60_h0_s2 S_age70_h0_s2,  ///
        at1(age 50 hormon 0 size 2)                 ///
        at2(age 60 hormon 0 size 2)                 ///
        at3(age 70 hormon 0 size 2)                 ///
        survival  ci                                ///
        timevar(0 10, step(0.1))                    ///
        frame(f1)                                   
frame f1: list S_age50_h0_s2 S_age60_h0_s2 S_age70_h0_s2 in 1/10predict S_age50_h0_s2 S_age60_h0_s2 S_age70_h0_s2,  ///
        at1(age 50 hormon 0 size 2)                 ///
        at2(age 60 hormon 0 size 2)                 ///
        at3(age 70 hormon 0 size 2)                 ///
        survival  ci                                ///
        timevar(0 10, step(0.1))                    ///
        frame(f1)                                   
frame f1: list S_age50_h0_s2 S_age60_h0_s2 S_age70_h0_s2 in 1/10predict S_age50_h0_s2 S_age60_h0_s2 S_age70_h0_s2,  ///
        at1(age 50 hormon 0 size 2)                 ///
        at2(age 60 hormon 0 size 2)                 ///
        at3(age 70 hormon 0 size 2)                 ///
        survival  ci                                ///
        timevar(0 10, step(0.1))                    ///
        frame(f1)                                   
frame f1: list S_age50_h0_s2 S_age60_h0_s2 S_age70_h0_s2 in 1/10predict S_age50_h0_s2 S_age60_h0_s2 S_age70_h0_s2,  ///
        at1(age 50 hormon 0 size 2)                 ///
        at2(age 60 hormon 0 size 2)                 ///
        at3(age 70 hormon 0 size 2)                 ///
        survival  ci                                ///
        timevar(0 10, step(0.1))                    ///
        frame(f1)                                   
frame f1: list S_age50_h0_s2 S_age60_h0_s2 S_age70_h0_s2 in 1/10predict S_age50_h0_s2 S_age60_h0_s2 S_age70_h0_s2,  ///
        at1(age 50 hormon 0 size 2)                 ///
        at2(age 60 hormon 0 size 2)                 ///
        at3(age 70 hormon 0 size 2)                 ///
        survival  ci                                ///
        timevar(0 10, step(0.1))                    ///
        frame(f1)                                   
frame f1: list S_age50_h0_s2 S_age60_h0_s2 S_age70_h0_s2 in 1/10predict S_age50_h0_s2 S_age60_h0_s2 S_age70_h0_s2,  ///
        at1(age 50 hormon 0 size 2)                 ///
        at2(age 60 hormon 0 size 2)                 ///
        at3(age 70 hormon 0 size 2)                 ///
        survival  ci                                ///
        timevar(0 10, step(0.1))                    ///
        frame(f1)                                   
frame f1: list S_age50_h0_s2 S_age60_h0_s2 S_age70_h0_s2 in 1/10predict S_age50_h0_s2 S_age60_h0_s2 S_age70_h0_s2,  ///
        at1(age 50 hormon 0 size 2)                 ///
        at2(age 60 hormon 0 size 2)                 ///
        at3(age 70 hormon 0 size 2)                 ///
        survival  ci                                ///
        timevar(0 10, step(0.1))                    ///
        frame(f1)                                   
frame f1: list S_age50_h0_s2 S_age60_h0_s2 S_age70_h0_s2 in 1/10

List Predictions


. frame f1: list tt S_age50_h0_s2 S_age60_h0_s2 S_age70_h0_s2 in 1/10,  ab(13)

     +----------------------------------------------------+
     | tt   S_age50_h0_s2   S_age60_h0_s2   S_age70_h0_s2 |
     |----------------------------------------------------|
  1. |  0               1               1               1 |
  2. | .1       .99996383       .99995863       .99994405 |
  3. | .2       .99973446       .99969626       .99958923 |
  4. | .3       .99916874        .9990492       .99871432 |
  5. | .4       .99817701       .99791498        .9971812 |
     |----------------------------------------------------|
  6. | .5        .9967111       .99623878       .99491659 |
  7. | .6       .99475256       .99399981       .99189377 |
  8. | .7       .99230356       .99120106       .98811855 |
  9. | .8        .9893804       .98786173       .98361906 |
 10. | .9       .98600871       .98401177        .9784382 |
     +----------------------------------------------------+

Plot Predictions

Predictions

frame f1: {
  twoway (line S_age50_h0_s2 S_age60_h0_s2 S_age70_h0_s2 tt), ///
         ytitle("Survival function")                          ///
         xtitle("Years since surgery")                        ///
         legend(order(1 "Age 50" 2 "Age 60" 3 "age 70")       ///
                ring(0) pos(1) cols(1))
}

Merge to Existing Frame

Predictions

predict S_age50_h1_s2 S_age60_h1_s2 S_age70_h1_s2,  ///
        at1(age 50 hormon 1 size 2)                 ///
        at2(age 60 hormon 1 size 2)                 ///
        at3(age 70 hormon 1 size 2)                 ///
        survival ci                                 ///
        frame(f1, merge)   predict S_age50_h1_s2 S_age60_h1_s2 S_age70_h1_s2,  ///
        at1(age 50 hormon 1 size 2)                 ///
        at2(age 60 hormon 1 size 2)                 ///
        at3(age 70 hormon 1 size 2)                 ///
        survival ci                                 ///
        frame(f1, merge)   predict S_age50_h1_s2 S_age60_h1_s2 S_age70_h1_s2,  ///
        at1(age 50 hormon 1 size 2)                 ///
        at2(age 60 hormon 1 size 2)                 ///
        at3(age 70 hormon 1 size 2)                 ///
        survival ci                                 ///
        frame(f1, merge)


. frame f1 {
.   list tt S_age50_h0_s2 S_age50_h1_s2 in 1/5, abbrev(13)

     +------------------------------------+
     | tt   S_age50_h0_s2   S_age50_h1_s2 |
     |------------------------------------|
  1. |  0               1               1 |
  2. | .1       .99996383       .99995634 |
  3. | .2       .99973446       .99967945 |
  4. | .3       .99916874       .99899661 |
  5. | .4       .99817701       .99779973 |
     +------------------------------------+
. }

Predict at single timepoint

Predictions

predict S10,           ///
        at1(.)         /// observed values
        survival       ///
        timevar(10)    /// single time value
        merge          //  merge in current frame predict S10,           ///
        at1(.)         /// observed values
        survival       ///
        timevar(10)    /// single time value
        merge          //  merge in current frame predict S10,           ///
        at1(.)         /// observed values
        survival       ///
        timevar(10)    /// single time value
        merge          //  merge in current frame

Type of Predictions

`cumhazard`	-	cumulative hazard function
`failure`	-	failure function
`hazard`	-	hazard function
`lnhazard`	-	ln(hazard) function
`rmft`	-	restricted mean failure time
`rmst`	-	restricted mean survival time
`survival`	-	survival function
`xb`	-	the linear predictor

Use ci or se for confidence intervals/standard errors.
use contrast(difference) or contrast(ratio) for contrasts.

Contrasts

You can have as many at options as you want.
Obtain differences or ratios between covariate patterns.

stpm3 i.hormon @ns(age, df(3)), scale(lncumhazard) df(5) 
predict S70h0 S70h1, survival ci                 ///
                     at1(age 70 hormon 0)        ///
                     at2(age 70 hormon 1)        ///
                     timevar(0 5, step(0.1))     ///
                     contrast(difference)        ///
                     contrastvar(Sdiff)          ///
                     frame(f1)stpm3 i.hormon @ns(age, df(3)), scale(lncumhazard) df(5) 
predict S70h0 S70h1, survival ci                 ///
                     at1(age 70 hormon 0)        ///
                     at2(age 70 hormon 1)        ///
                     timevar(0 5, step(0.1))     ///
                     contrast(difference)        ///
                     contrastvar(Sdiff)          ///
                     frame(f1)stpm3 i.hormon @ns(age, df(3)), scale(lncumhazard) df(5) 
predict S70h0 S70h1, survival ci                 ///
                     at1(age 70 hormon 0)        ///
                     at2(age 70 hormon 1)        ///
                     timevar(0 5, step(0.1))     ///
                     contrast(difference)        ///
                     contrastvar(Sdiff)          ///
                     frame(f1)stpm3 i.hormon @ns(age, df(3)), scale(lncumhazard) df(5) 
predict S70h0 S70h1, survival ci                 ///
                     at1(age 70 hormon 0)        ///
                     at2(age 70 hormon 1)        ///
                     timevar(0 5, step(0.1))     ///
                     contrast(difference)        ///
                     contrastvar(Sdiff)          ///
                     frame(f1)

Set which at option is the reference with atreference()

Plot Contrasts

frame f1 {
  twoway (rarea Sdiff_lci Sdiff_uci tt, color(red%30)) ///
         (line Sdiff tt, color(red))                   ///
         , xtitle("Years since surgery")               ///
         ytitle("Difference in S(t)")                  ///
         ylabel(,format(%3.2f))                        ///
         legend(off)
}

Standardization

This was the main motivation for introducing factor variables and extended functions.
standsurv is margins for survival data.
Interactions with exposures made code complicated and prone to error when using standsurv

Regression Standardization

Fit regression model with exposure, \(X\), and confounders, \(Z\)
Predict outcome if all individuals were exposed (\(X=0\)) and unexposed (\(X=1\))
Perform contrast \[ \frac{1}{N} \sum_{i=1}^N {\widehat{S}(t|X=1,Z=z_i)} - \frac{1}{N} \sum_{i=1}^N {\widehat{S}(t|X=0,Z=z_i)} \]
Similar ideas in competing risks, relative survival, multistate models

Models with interactions (`stpm2`)

use https://www.pclambert.net/data/rott2b, clear
stset os, f(cause==1,2) scale(12) exit(time 365.24*10)

// generate splines
rcsgen age,  df(3) gen(agercs)
rcsgen pr,   df(3) gen(prrcs)
// interactions
forvalues i = 1/3 {
  gen hagercs`i'  = hormon*agercs`i'
  gen hprrcs`i'   = hormon*prrcs`i'
}

// dummy variables for size
tab size, gen(size)
// hormon - size interactions
gen hsize2 = hormon*size2
gen hsize3 = hormon*size3

Fit Model in `stpm2`

stpm2 agercs1 agercs2 agercs3      /// age splines
      prrcs1 prrcs2 prrcs3         /// nodes splines
      size2 size3                  /// size dummy variables
      hormon                       /// exposure
      hagercs1 hagercs2 hagercs3   /// exposure age interactions
      hprrcs1 hprrcs2 hprrcs3      /// exposure nodes interactions
      hsize2 hsize3                /// exposure size interactions
      , scale(hazard)              /// 
        df(4)                      ///
        tvc(hormon                 /// time interactions
        agercs1 agercs2 agercs3    ///
        prrcs1 prrcs2 prrcs3       ///
        size2 size3)               /// 
        dftvc(3)

Using `standsurv` after `stpm2`

standsurv is a postestimation command that estimates marginal effects and contrasts.

standsurv, at1(hormon 1                                            ///
               hagercs1=agercs1 hagercs2=agercs2 hagercs3=agercs3  ///
               hprrcs1=prrcs1 hprrcs2=prrcs2 hprrcs3=prrcs3        ///
               hsize2=size2 hsize3=size3)                          ///
           at2(hormon 0                                            ///
               hagercs1 0 hagercs2 0 hagercs3 0                    ///
               hprrcs1 0 hprrcs2 0 hprrcs3 0                       ///
               hsize2 0 hsize3 0)                                  ///
           atvar(Sh1 Sh0)                                          ///    
           contrast(difference)                                    ///
           contrastvar(Sdiff)                                      ///
           timevar(tt) ci

Fit Model in `stpm3`

stpm3 i.hormon##(@ns(age, df(3))      ///
                 @ns(pr, df(3))       ///
                 i.size),             ///
                 scale(lncumhazard)   ///
                 df(4)                ///
                 tvc(i.hormon         ///
                     @ns(age, df(3))  ///
                     @ns(pr, df(3))   ///
                     i.size)          ///
                 dftvc(3)

Using `standsurv` after `stpm3`

standsurv, at1(hormon 1)            ///
           at2(hormon 0)            ///
           atvar(Sh1 Sh0)           ///    
           contrast(difference)     ///
           contrastvar(Sdiff)       ///
           timevar(tt) ci

Competing Risks Example

// Load Data
use https://www.pclambert.net/data/rott2b

// Cancer Model
stset os, f(cause==1) scale(12) exit(time 365.24*10)
stpm3 i.hormon i.size enodes @ns(age,df(3))  ///
      i.hormon#i.size i.hormon#c.enodes      ///
      i.hormon#@ns(age,df(3)),               ///
      scale(logcumhazard) df(4) 
estimates store cancer

Competing Risks Example

// Load Data
use https://www.pclambert.net/data/rott2b

// Cancer Model
stset os, f(cause==1) scale(12) exit(time 12*10)
stpm3 i.hormon i.size enodes @ns(age,df(3))  ///
      i.hormon#i.size i.hormon#c.enodes      ///
      i.hormon#@ns(age,df(3)),               ///
      scale(logcumhazard) df(4) 
estimates store cancer           
      
// Other cause Model
stset os, f(cause==2) scale(12) exit(time 12*10)
stpm3 i.hormon i.size enodes @ns(age,df(3))  ///
      i.hormon#i.size i.hormon#c.enodes      ///
      i.hormon#@ns(age,df(3)),               ///
      scale(loghazard) df(4) 
estimates store other

Using `standsurv`

range tt 0 10 101
standsurv, cif crmodels(cancer other) ci    ///
           timevar(tt)                      ///
           at1(hormon 0)                    ///
           at2(hormon 1)                    ///
           atvar(CIF0 CIF1)                 ///
           contrast(difference)             ///
           contrastvar(CIFdiff)  range tt 0 10 101
standsurv, cif crmodels(cancer other) ci    ///
           timevar(tt)                      ///
           at1(hormon 0)                    ///
           at2(hormon 1)                    ///
           atvar(CIF0 CIF1)                 ///
           contrast(difference)             ///
           contrastvar(CIFdiff)  range tt 0 10 101
standsurv, cif crmodels(cancer other) ci    ///
           timevar(tt)                      ///
           at1(hormon 0)                    ///
           at2(hormon 1)                    ///
           atvar(CIF0 CIF1)                 ///
           contrast(difference)             ///
           contrastvar(CIFdiff)  range tt 0 10 101
standsurv, cif crmodels(cancer other) ci    ///
           timevar(tt)                      ///
           at1(hormon 0)                    ///
           at2(hormon 1)                    ///
           atvar(CIF0 CIF1)                 ///
           contrast(difference)             ///
           contrastvar(CIFdiff)  range tt 0 10 101
standsurv, cif crmodels(cancer other) ci    ///
           timevar(tt)                      ///
           at1(hormon 0)                    ///
           at2(hormon 1)                    ///
           atvar(CIF0 CIF1)                 ///
           contrast(difference)             ///
           contrastvar(CIFdiff)

Competing Risks

Installation

Still running various checks etc.
A test version is available

Installation

net install stpm3_test, from (https://www.pclambert.net/downloads/stpm3_test)

net install standsurv, from (https://www.pclambert.net/downloads/standsurv_test)

net install gensplines, from (https://www.pclambert.net/downloads/gensplines)

See pclambert.net/software/stpm3

References

Lambert, P. C., and P. Royston. 2009. “Further Development of Flexible Parametric Models for Survival Analysis.” The Stata Journal 9: 265–90.

Royston, P. 2001. “Flexible Parametric Alternatives to the Cox Model, and More.” The Stata Journal 1: 1–28.

Royston, P., and P. C. Lambert. 2011. Flexible Parametric Survival Analysis in Stata: Beyond the Cox Model. Stata Press.

Syriopoulou, Elisavet, Sarwar I. Mozumder, Mark J. Rutherford, and Paul C. Lambert. 2022. “Estimating Causal Effects in the Presence of Competing Events Using Regression Standardisation with the Stata Command standsurv.” BMC Medical Research Methodology 22 (August): 226. https://doi.org/10.1186/s12874-022-01666-x.

Wang, Wenjie, and Jun Yan. 2021. “Shape-Restricted Regression Splines with R Package splines2.” Journal of Data Science 19 (3): 498–517. https://doi.org/10.6339/21-JDS1020.