Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • Paul Dickman
    replied
    Thanks Steve.

    This code is from an exercise Paul Lambert and I use in short courses. I'm not sure of the problem, but this is the latest version of the code for this exercise and it runs for me. I believe Paul Lambert (author of -stpm2- and co-author of this exercise) is looking into it. For the moment, could you confirm that this code runs for you.

    Neither Paul Lambert or I have great expertise in methods for multiple imputation, so I'm happy to receive any comments on how the code could be improved.

    Code:
    //======================================================================//
    // EXERCISE 285: 
    // From short-course on population-based cancer patient survival
    // http://cansurv.net/
    // October 2015
    // Paul Dickman
    //======================================================================//
    set more off
    version
    which stpm2
    
    use http://pauldickman.com/survival/colon, 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/survival/popmort
    keep if _merge==3
    
    tab stage
    
    /* Check stage distribution over age and gender */
    tab stage agegrp, column
    tab stage sex, column
    
    /* Graphs of survival by age group and stage */
    stpm2 ib1.stage##i.agegrp , df(5) bhaz(rate) scale(hazard) eform nolog
    predict survival, surv
    
    line survival _t if stage==0, lpattern(dash) sort || ///
    line survival _t if stage==1, sort || ///
    line survival _t if stage==2, sort || ///
    line survival _t if stage==3, sort by(agegrp) ///
     legend(order(1 "Unknown" 2 "Localised" 3 "Regional" 4 "Distant"))
    
    /* Fit model using missing indicator approach */
    stpm2 ib1.stage i.agegrp , df(5) bhaz(rate) scale(hazard) eform nolog
    
    /* Refit model using complete records approach */
    replace stage=. if stage==0
    stpm2 ib1.stage i.agegrp , df(5) bhaz(rate) scale(hazard) eform nolog
    
    // The outcome should be included in the imputation model.
    // Falcaro et al suggest the Nelson-Aalen estimate of cum. hazard.
    sts gen H=na
    
    // Declare multiple-imputation data
    mi set wide
    
    //  Declare which variables have missing data
    mi register imputed stage
    
    // Perform imputation. 
    // This creates 10 additional copies of the obs with missing stage.
    // Carpenter and Kenward (2013) suggest 30 imputations. We use 10.
    set seed 29390
    mi impute chained (mlogit) stage = i.subsite sex i.agegrp H _d, add(10)
    
    list id _mi_m agegrp sex stage _t _d if id==2287
    list id _mi_m agegrp sex stage _t _d if id==3362
    list id _mi_m agegrp sex stage _t _d if id==3501
    
    // Note that as stpm2 is not an official Stata command we need to use cmd ok option
    // Save the estimates so that they can be used for making predictions 
    mi estimate, dots cmdok sav(mi_stpm2,replace): ///
        stpm2 ib1.stage i.agegrp, df(5) bhaz(rate) scale(hazard) nolog eform
    
    // predict survival using -mi predictnl-    
    mi predictnl survimp2 = predict(survival at(agegrp 2)) using mi_stpm2
    
    // compare predictions to complete case analysis
    stpm2 ib1.stage i.agegrp if _mi_m==0, df(5) scale(h) bhaz(rate)
    predict surv, survival at(agegrp 2)
    
    line surv survimp2 _t if stage==1 & _mi_m==0, sort || ///
    line surv survimp2 _t if stage==2 & _mi_m==0, sort || ///
    line surv survimp2 _t if stage==3 & _mi_m==0, sort ///
    title("Predicted survival for agegrp==2 (60-74)") ///
    legend(order(1 "Localised (Complete)" 2  "Localised (Imputed)" ///
                  3 "Regional (Complete)" 4 "Regional (Imputed)" ///
                  5 "Distant (Complete)" 6 "Distant (Imputed)")) ///
               name(imputed, replace)
               
    // predict 5-year net survival using meansurv approach
    gen t5=5
    predict ns5, survival timevar(t5)
    mi predictnl ns5imp = predict(survival timevar(t5)) using mi_stpm2
    table stage agegrp, contents(mean ns5 mean ns5imp) format(%5.3f)
    table stage, contents(mean ns5 mean ns5imp) format(%5.3f)
    table agegrp, contents(mean ns5 mean ns5imp) format(%5.3f)

    Leave a comment:


  • Steve Samuels
    replied
    Thank you, Laura. I've written to the authors about this.

    Leave a comment:


  • Laura Myles
    replied
    Hi Steve,

    Thanks for your suggestion. I will have a go!

    The example I found is at the end of this document (page 250 onwards).

    http://www.pauldickman.com/survival/.../solutions.pdf

    Laura



    Leave a comment:


  • Steve Samuels
    replied
    Sorry, mi estimate isn't necessary and will trigger an error, because the extracted data set is not mi set. This should be okay.
    Code:
    tempfile t1
    save `t1', replace
    save results, emptyok replace
    forvalues i = 1/10 {
        use `t1', clear
        mi extract `i'
        stpm2 ib1.stage i.agegrp, df(5) bhaz(rate) scale(hazard) nolog eform
        predict surv, survival at(agegrp  2)
        gen dset = `i'
        keep  dset surv valve _t agegrp
        append using results
        save results , replace
    }
    use results, clear
    collapse (mean) surv,  by (agegrp _t)
    sort_t
    twoway connect surv _t

    Leave a comment:


  • Steve Samuels
    replied
    Here's a round-about approach that gets the prediction for every completed data imputation, then takes the mean:
    Code:
    tempfile t1
    save `t1', replace
    save results, emptyok replace
    forvalues i = 1/10 {
                use `t1', clear
                mi extract `i'
                mi estimate, dots cmdok sav(mi_stpm2,replace): ///
                     stpm2 ib1.stage i.agegrp, df(5) bhaz(rate) scale(hazard) nolog eform
                predict surv, survival at(agegrp  2)
                gen dset = `i'
                keep  dset surv valve _t agegrp
                append using results
                save results , replace
    }
    use results, clear
    collapse (mean) surv,  by (agegrp _t)
    sort_t
    twoway connect surv _t

    Leave a comment:


  • Steve Samuels
    replied
    Perhaps predict will work (not tried):
    Code:
    mi predict survimpt using mi_stpm2, survival  at(agegrp 2)
    ADDED: predict doesn't work either. The error message is "option survival" not found.
    I'm curious. Can you provide a link to one of the online examples that use predictnl? Thanks.
    Last edited by Steve Samuels; 01 Nov 2018, 13:27.

    Leave a comment:

Working...
X