Predicting stcurve survival curve for three-way interaction (with quadratic)

[David Lucas]

Hello! I have a dataset of firms with founding and exit (dissolution) dates as follows:

* Example generated by -dataex-. For more info, type help dataex
clear
input float(founded_date x2_continuous dissolved duration date_exit x_binary companynumber)
21223  35.87123 0 1818 23041 0   1
21267  52.30959 1  948 22215 0   2
19463  44.91233 1 1604 21067 0   3
20587 36.237705 1  886 21473 0   4
20698  46.83601 0 2343 23041 1   5
19709  47.42101 0 3332 23041 0   6
19030  48.43716 1 1904 20934 0   7
19984  24.45114 0 3057 23041 0   8
21283  30.92603 0 1758 23041 0   9
20523  58.78224 0 2518 23041 0  10
19940  39.16256 1 1162 21102 0  11
20368   59.2664 1 1420 21788 0  12
21236  35.60137 1  517 21753 0  13
19681  49.30137 1 1827 21508 0  14
20602  56.44718 1  521 21123 0  15
19966  53.45662 1 1143 21109 0  16
20591  53.66667 1  735 21326 1  17
21271  54.73973 1  944 22215 0  18
21235 38.015068 0 1806 23041 1  19
19109  52.77869 1  691 19800 0  20
19064 37.820583 0 3977 23041 0  21
20075  61.08767 1  747 20822 0  22
21040  42.89589 0 2001 23041 0  23
20364  42.50496 1  479 20843 0  24
19564  42.77206 0 3477 23041 0  25
19117 29.547813 1  802 19919 1  26
20632   73.2795 0 2409 23041 0  27
19520  29.31918 1  238 19758 0  28
20618  50.32513 0 2423 23041 0  29
19746 32.934246 1 2770 22516 0  30
19201  45.32192 1  844 20045 0  31
19715  47.81233 1 3207 22922 0  32
20495  54.61339 1 1699 22194 0  33
20900  53.84932 0 2141 23041 0  34
20207  61.72874 0 2834 23041 0  35
20349   44.5065 0 2692 23041 0  36
20546  43.46175 0 2495 23041 0  37
19710   50.2589 1  888 20598 0  38
20975   52.4137 0 2066 23041 0  39
20039  32.98904 0 3002 23041 1  40
19500  52.62648 1 1783 21283 0  41
21214  34.29041 1  973 22187 0  42
19530  50.17671 0 3511 23041 0  43
21292  55.79452 0 1749 23041 0  44
19985  43.17808 1  732 20717 0  45
20240  44.77468 1 1464 21704 0  46
20376  63.07877 1 1973 22349 0  47
20220  47.23699 0 2821 23041 0  48
19810 36.863014 1  599 20409 0  49
19624   44.7726 1 1415 21039 0  50
20214  38.63562 1 1546 21760 0  51
19884  53.06575 1  525 20409 0  52
19537  43.81004 1  592 20129 0  53
20256  49.62841 1  615 20871 1  54
20708 32.960068 0 2333 23041 0  55
19858  53.46137 1  572 20430 0  56
21077  44.37534 1 1453 22530 0  57
20067  51.17123 1 1924 21991 0  58
19236  47.87432 1 1418 20654 0  59
19841  43.78219 1  589 20430 0  60
19193  48.53279 1 1160 20353 0  61
20257  54.38082 0 2784 23041 1  62
20708 34.682594 1  464 21172 0  63
20173  57.47032 1 1790 21963 0  64
19243  43.14418 0 3798 23041 0  65
21292  47.44749 1 1413 22705 0  66
19396   52.4516 0 3645 23041 0  67
20145  47.86712 0 2896 23041 0  68
19578 36.840183 0 3463 23041 0  69
19584  55.45206 1 1203 20787 0  70
19967  38.46301 1  589 20556 1  71
21256  34.09954 0 1785 23041 0  72
19965  37.37123 1  892 20857 0  73
20877  46.08858 1 1303 22180 1  74
19905  42.20548 0 3136 23041 0  75
20130  38.98959 1 1168 21298 0  76
20880   40.6274 0 2161 23041 0  77
19976  34.69452 0 3065 23041 0  78
19143  58.16301 1  587 19730 0  79
19894  56.00959 1  284 20178 0  80
19871 25.447945 1  727 20598 0  81
19845 28.375343 1  592 20437 0  82
21123  41.50137 1  777 21900 0  83
20944  28.29863 1 1327 22271 0  84
19856 36.197258 1  903 20759 0  85
20472  35.23907 1 2163 22635 0  86
19109   38.1102 1  586 19695 0  87
19410  27.80822 1  845 20255 0  88
19799  66.75069 0 3242 23041 0  89
21227  33.49315 0 1814 23041 0  90
19400  43.76714 1 1198 20598 0  91
19304  49.39294 1 2393 21697 0  92
21203  40.26302 1 1025 22228 0  93
20927 32.946117 1 1834 22761 0  94
20654  47.96586 1 1239 21893 1  95
21267 34.267124 1  983 22250 1  96
20294  56.83473 0 2747 23041 0  97
20818   56.9137 0 2223 23041 1  98
20368  60.10083 1  524 20892 0  99
21167  39.49589 0 1874 23041 0 100
end
format %td founded_date

I am estimating a Cox proportional hazards model, where I want to interact a binary (x2_binary) variable with the quadratic specification of a continuous (x2_continuous) variable. Essentially, the interaction is i.x2_binary##c.x2_continuous##c.x2_continuous.

The stcox estimation seems fine. However, I'm struggling to produce estimates from stcurve at selected levels of this quadratic/three-way interaction. I note that the two-way interaction of x2_binary * x2_continuous runs without issue for both margins (hazard prediction) and stcurve:

//stset the above data
stset date_exit, failure(dissolved==1) id(companynumber) scale(365.25) origin(time founded_date) exit(time mdy(01,31,2023))

//linear interaction - both margins and stcurve produce results
stcox i.x_binary##c.x2_continuous, vce(robust)
margins, at(x2_continuous = (18 64) x_binary = (0 1))
stcurve, survival at(x2_continuous = (18 64) x_binary = (0 1))

So far, so good. However, when I include the quadratic specification (x2_binary * x2_continuous * x2_continuous), the margins calculation continues to work—but the stcurve fails to calculate:

//quadratic (three-way) interaction
stcox i.x_binary##c.x2_continuous##c.x2_continuous, vce(robust)
margins, at(x2_continuous = (18 64) x_binary = (0 1)) //works
stcurve, survival at(x2_continuous = (18 64) x_binary = (0 1)) //fails

As noted, margins will yield a prediction for the hazard. But stcurve does not:

the model contains an interaction involving factor variables; the level for each variable in the interaction must be set with the at() option
r(198);

Is there an at( ) specification that I'm missing? Any insights on this would be appreciated. Thank you!