I am running a flexible parametric survival model using stpm2.

In my data, I have two variables that will be interacted with each other. Both of these variables have proportional hazards assumption violated. By interacting these two variables with each other, plus time-varying effects for both, I gather this would result in a three-way interaction as mentioned somewhere else on this forum.

How can I get the proper time-varying hazard ratios from stpm2?

I am able to do this with stcox (assuming proportional hazards), but I cannot figure out how to do the analogous version in stpm2, with interaction terms incorporated into the predicted hazard ratios.

Sample code below. Please help!

- Ash

Code:
// load example data

use "http://www.stata-press.com/data/fpsaus/ew_breast_ch7.dta", clear
stset survtime, failure(dead==1) exit(time 5) id(ident)

// binary variables

generate x=0 if dep5==0
replace x=1 if dep5==1

generate z=0 if agediag<40
replace z=1 if agediag>=40

// stratified cox model with different baselines for each level of z

stcox x if z==0
stcox x if z==1

// or, equivalently

stcox x##z, strata(z)
lincom _b[1.x]+_b[1.x#0.z], eform
lincom _b[1.x]+_b[1.x#1.z], eform

// Now, how can I obtain the equivalent stratified hazard ratios in stpm2, allowing for tvc for both the x and z variable?
// I can run separate models...

stpm2 x if z==0, scale(hazard) df(3) tvc(x) dftvc(5) eform
predict hr_at0, hrnumerator(x 1) hrdenominator(x 0)

stpm2 x if z==1, scale(hazard) df(3) tvc(x) dftvc(5) eform
predict hr_at1, hrnumerator(x 1) hrdenominator(x 0)

twoway (line hr_at0 _t, sort) (line hr_at1 _t, sort), yscale(log)

drop hr*

// But, how do I do this using interaction terms in the same model and properly incorporate them into the predict command?
// Is my code below correct? It gives the same graphs as the stratified stpm2 models...?

generate x_by_z=x*z

stpm2 x z x_by_z, scale(hazard) df(3) tvc(x z x_by_z) dftvc(5) eform
predict hr_at0, hrnumerator(x 1 z 0 x_by_z 0) hrdenominator(x 0 z 0 x_by_z 0)
predict hr_at1, hrnumerator(x 1 z 1 x_by_z 1) hrdenominator(x 0 z 1 x_by_z 0)

twoway (line hr_at0 _t, sort) (line hr_at1 _t, sort), yscale(log)

drop hr*