/***************************************************************************
This code available at:
http://pauldickman.com/software/stnet/compare_stns_uk_melanoma.do

Comments available at:
http://pauldickman.com/software/stnet/compare_stns/

Compare estimates of net survival (Pohar Perme) between stnet, strs, stpp, and 
stns using the English melanoma data from the LSHTM CSG short course (data are
available to short course participants at the short course home page).

Note: code for stnet not yet available.

1. stnet (https://pauldickman.com/software/stnet/)
2. strs  (https://pauldickman.com/software/strs/)
3. stpp  (https://pclambert.net/software/stpp/)
4. stns  (https://www.stata-journal.com/article.html?article=st0326)

Installation of user-written commands
ssc install stnet
net install https://www.pauldickman.com/strs/strs, all
ssc install stpp
ssc install stns

Paul Dickman, June 2021
***************************************************************************/

// Modify the life table (stns requires days as the time unit)
use http://pauldickman.com/data/Lifetable_2013, clear
gen double rate_day=rate/365.25
label variable rate_day "Expected mortality rate per day"
generate _agedays=_age*365.25
label variable _agedays "Age in days"
generate yearindays=mdy(1,1,_year)
label variable yearindays "Year as a Stata date (days since 1-1-1960)"
sort _year sex _age dep
save Lifetable_2013_stns, replace

// MAIN ANALYSIS
use http://pauldickman.com/data/melanoma_2013, clear
set seed 20210610 // so that the sampling (next line) is reproducible
sample 10 // Work with a 10% sample to save time
gen agediagindays=(diagmdy-birthmdy)

// Using stns
// stset with time in days
stset ftime, failure(dead==1) id(id)
stns list using Lifetable_2013_stns.dta,age(agediagindays=_agedays) ///
period(diagmdy=yearindays) strata(sex dep) rate(rate_day) ///
at(1(1)10, scalefactor(365.25) method(step))

// Using strs
// stset with time in years
stset finmdy, failure(dead==1) origin(diagmdy) scale(365.24) id(id)
strs using Lifetable_2013_stns, br(0(.08333333)10) mergeby(_year sex _age dep) ///
    diagage(agediag) diagyear(ydiag) pohar notables save(replace)

// Using stpp
stpp R_pp using Lifetable_2013_stns, agediag(agediag) datediag(diagmdy) ///
     pmother(sex dep) list(1(1)10)   	

// output summary of results for strs	 
use grouped, clear
list end n d cr_e2 cns_pp lo_cns_pp hi_cns_pp if inlist(end,5,10)