stpm2 and mi predictnl

[Laura Myles]

Dear Listers,

I am running survival analysis using flexible parametric models in -stpm2-.

I have missing data in one of my covariates. As part of my sensitivity analyses, I am imputing this missing information using multiple imputation. I would like to produce survival plots based on my imputed results and I have found examples online showing this could be achieved using -mi predictnl- just after sptm2. However, this is not straightforward once I run it in on my data so I thought I tried on the original sample data (melanoma.dta).


I use the following code:

replace stage = . if stage==0

mi set flong
mi register regular subsite agegrp sex
mi register imputed stage
mi impute chained (mlogit) stage = i.subsite sex i.agegrp , add(10)

mi stset surv_mm, fail(status==1 2) id(id) scale(12)
mi stset surv_mm, failure(status=1,2) scale(12) id(id) exit(time 120.5)

mi estimate, dots cmdok sav(mi_stpm2,replace): ///
stpm2 ib1.stage i.agegrp, df(5) bhaz(rate) scale(hazard) nolog eform

mi predictnl survimp2 = predict(survival at(agegrp 2)) using mi_stpm2

On this sample data and my data, I get the following error message:

variable _rcs1 not found
xb(xb) invalid
predict(survival at(agegrp 2)) invalid


_rcs1 is part of the survival analysis output so I am not sure what could be going wrong. I tried updating my stato ado files (update all) but this has not resolved the issue. Any suggestion would be greatly appreciated.

[Miguel Sanchez]

I wrote to Paul Lambert about this error, and he replied:

"You can’t add a ci option here, but should add it as an option of the mi predictnl command, e.g., ci(lci uci)"
I hadn't realized that ci was an option of predictnl. The entire command also requires a force option, so the correct version would be:

mi predictnl survimp2 = predict(survival at(agegrp 2) timevar(_t)) using mi_stpm2, ci(surv_ll surv_ul) force

Hi, following on this one - is it possible here to change the method of the ci calculation? E.g., Wilson or Jeffreys? I have negative values from the CIs in my prediction, and thus I cannot plot them using a log scale. I'm trying alternatives...

Thanks!

[Pavel Roshanov]

just following up on this in case anyone has suggestions

[Pavel Roshanov]

I wonder if anyone has any suggestions for this in a multiple imputation context.

Thanks very much,
Pavel

[Pavel Roshanov]

Absolutely. Here is one modifying your code on the colon cancer dataset:

set more off
use http://pauldickman.com/data/colon.dta, clear

stset surv_mm, failure(status=1 2) scale(12) exit(time 10*12)
gen _age = min(int(age + _t),99)
gen _year = int(yydx + mmdx/12 + _t)
merge m:1 _year sex _age using http://pauldickman.com/data/popmort
keep if _merge==3
replace stage=. if stage==0

sts gen H=na

gen timevar = .
replace timevar = 1 in 1
replace timevar = 2 in 2
replace timevar = 3 in 3
replace timevar = 4 in 4
replace timevar = 5 in 5
replace timevar = 6 in 6
replace timevar = 7 in 7
replace timevar = 8 in 8
replace timevar = 9 in 9
replace timevar = 10 in 10

mi set flong
mi register imputed stage



set seed 29390
mi impute chained (mlogit) stage = i.subsite sex i.agegrp H _d, add(10)

tab(agegrp),gen(agegrp)
tab(stage),gen(stage)

mi estimate, dots cmdok sav(mi_stpm2,replace): ///
    stpm2 stage1 sex, df(5) bhaz(rate) scale(hazard) nolog  ///
    tvc(stage1) dftvc(stage1:2)
    
mi predictnl hrimp = predict(hr hrnum(stage1 1)  hrdenom(stage1 0) timevar(timevar))  using mi_stpm2, ci(mill miul) force
/*output here is:
option hr not allowed
predict(hr hrnum(stage1 1)  hrdenom(stage1 0) timevar(timevar)) invalid
*/

// compare predictions to first imputation
stpm2 stage1 sex if _mi_m==1, df(5) bhaz(rate) scale(hazard) nolog  ///
    tvc(stage1) dftvc(stage1:2)

//using predict
predict hrcc if _mi_m==1, hrnum(stage1 1) hrdenom(stage1 0) timevar(timevar) ci

//using standsurv
 standsurv if _mi_m==1,  at1(stage1 0 ) at2(stage1 1) ci contrast(ratio) hazard  contrastvar(hratio_stage1) timevar(timevar)
 
  list timevar hrcc hrcc_lci hrcc_uci  hratio_stage1 hratio_stage1_lci hratio_stage1_uci  if inlist(timevar,1,2,3,4,5,6,7,8,9,10) & _mi_m==1
  
/*output here shows that standsurv and predict would give nearly identical results:

        | timevar        hrcc    hrcc_lci    hrcc_uci   hratio_~1   hrati~lci   hrati~uci |
        |---------------------------------------------------------------------------------|
 15565. |       1   .10844999   .09978368   .11786898   .10849877   .09983011   .11792017 |
 15566. |       2   .15890791   .14671461   .17211458   .15901694   .14681975   .17222743 |
 15567. |       3   .21908062   .19810078   .24228233   .21925888   .19826885   .24247105 |
 15568. |       4   .29510389   .25943194   .33568075   .29536435   .25966927   .33596621 |
 15569. |       5    .3795031   .32253015   .44653997   .37985286   .32283667   .44693868 |
        |---------------------------------------------------------------------------------|
 15570. |       6   .46547052   .38154234   .56786045   .46591051    .3819134   .56838174 |
 15571. |       7   .54522721   .43228679   .68767477   .54575086   .43271411   .68831589 |
 15572. |       8   .61209619   .47285049   .79234717   .61269039   .47332405   .79309198 |
 15573. |       9   .66226365   .50314812   .87169787   .66291154   .50365777   .87252045 |
 15574. |      10   .69521221    .5241272   .92214259   .69589645   .52466367   .92301392 |
        +---------------------------------------------------------------------------------+
*/

[Paul Dickman]

Providing a fully-worked example including data would make it easier to answer.

[Pavel Roshanov]

Hello,

I am trying to estimate time-dependent marginal hazard ratios in multiply imputed data.

I have been doing this successfully in complete-case data with standsurv:

standsurv,  at1(var1 0) at2(var1 1) ci contrast(ratio) hazard  contrastvar(hratio_var1) timevar(_t)

I have tried using mi predictnl with no luck:

mi predictnl hrimp = predict(hr hrnum(var1 1)  timevar(_t) )  using mi_stpm2 , ci(lci uci)  force

Can anyone suggest code to do this?

In a related matter, I have also been trying to calculate time-dependent attributable fractions in MI data with no success but can do so in complete-case data.

Surely others would find the answer useful.

Thank you,
Pavel

[Paul Lambert]

Hi,

I agree with Paul Dickman that using the timevar() option is often preferable, but you do need to be careful here. We did not see what you wanted to with the predictions. Your original question was around the following

 
 mi predictnl survimp2 = predict(survival at(agegrp 2) timevar(_t)) using mi_stpm2

and we found a solution on how to replicate this. With this you could plot survimp2 vs _t at each level of stage, e.g.

line survimp2 _t if stage==1

There is an explanation of using the timevar() option here. https://pclambert.net/software/stpm2/stpm2_timevar/

In the recent code you used the meansurv() option. This predicts marginal survival functions rather than conditional survival functions. Essentially it predicts a survival curve for each individual in your study and then takes an average of the N survival curves. Use of the at() option forces the specified covariates to take specific values.

So in most cases when you use the timevar() option, if you want conditional survival curves, (i.e. for a specific covariate pattern), you should specify all covariate values in the at() option or potentially use the zeros option as well to force unspecified covariates to be zero.

When using meansurv you only want to specify covariates which you want to force to take specific values, for example when making contrasts between marginal survival functions.

Paul

[Laura Myles]

Paul Dickman Thanks again for the clarification!

[Paul Dickman]

Paul - I would like to check with you. In a previous thread (https://www.statalist.org/forums/for...r-net-survival), I was advised to create a timevar rather than using _t when using -predict- as below:

range timevar_new 0 10 100
predict s1_new, meansurv (at sex 1) timevar(timevar_new)

Is this still an option when working with imputed data or should I use _t?

Yes, this is the recommended option. By default, -predict- provides predictions for all observations at values of _t. By specifying -timevar(_t)- we are just forcing what should be the default behaviour. I made the predictions at _t to avoid adding extra complexity to the code. The dataset in my example has over 15000 observations and we don't need that many predictions so in practice I would use another temptime variable, such as in your example where we only predict for 100 values of time. One of the reasons it took so long to spot this bug is that most of us routinely use temptime when making predictions.

[Paul Dickman]

Thanks Steve Samuels for the nice words about our course. Thanks Laura Myles for identifying this problem. Here is the link to the exercises for our course. I have updated the exercise (285) on multiple imputation. The do file for that exercise can be found here.

[Laura Myles]

Thank you all again for your help in resolving my query!

Paul - I would like to check with you. In a previous thread (https://www.statalist.org/forums/for...r-net-survival), I was advised to create a timevar rather than using _t when using -predict- as below:
range timevar_new 0 10 100
predict s1_new, meansurv (at sex 1) timevar(timevar_new) Is this still an option when working with imputed data or should I use _t? Thanks again.

[Steve Samuels]

Sorry, try this link.

[Matt Warkentin]

Steve Samuels I have been following this thread closely. The hyperlink for the course notes does not seem to work for me. Could you provide the link again?

[Steve Samuels]

I want to thank Laura, for bringing this problem to the list, and the two Pauls for contributing their great knowledge to this problem. I urge anyone who aspires to use survival techniques to read the notes and exercises for their short course, available here.
Information about the course itself can be found here.