Improving fitting and predictions for flexible parametric survival models

Northern European Stata Conference 2022, Oslo

Paul Lambert

University of Leicester / Karolinska Institutet

October 12, 2022

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
------------------------------------------------------------------------------

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 & df


Interactions
stpm3 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
      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 & 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 

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 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 & 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/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)   

. 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

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)
  • 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)

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)

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.