// Load the Melanoma data, keep those with localized stage
frames reset
use melanoma, clear
gen female = sex == 2
stset surv_mm, failure(status==1,2) exit(time 120.5) scale(12)

// (a) Fit a model with 4 df
stpm3 @ns(age,df(3)) i.stage i.female, scale(lncumhazard) df(5)


// There will be a different survival function for each combination of
// age, stage and sex.
// We will look at the variation in predicted survival at 1, 5 and 10 years.
// Predict survival for all subjects given their covariate pattern
// at 1, 5 and 10 years
gen t1  = 1
gen t5  = 5
gen t10 = 10


// We use the a "." in teh at options which means to predict at the 
// observed values of the covariates.
// Note we can pass a dofferent timevar to each at() option.
predict S1 S5 S10, surv frame(S_compare,replace) ///
                   at1(., attimevar(t1))         ///
                   at2(., attimevar(t5))         ///
                   at3(., attimevar(t10))            

// We can then plot a histogram to show the varition in predicted probabilities                   
frame S_compare {
  hist S1, name(S1, replace)
  hist S5, name(S5, replace)
  hist S10, name(S10, replace)
}
graph combine S1 S5 S10, xcommon nocopies


// 
frame S_compare: summ S*
sts list, at(1 5 10)

range tt 0 10 101
standsurv, at1(.) timevar(tt) atvar(S_stpm3)

sts graph, addplot(line S_stpm3 tt)
list tt S_stpm3 if inlist(tt,1,5,10)
frame S_compare: summ S*


frame put age stage female, into(meansurv)

frame meansurv {
  gen female_copy = female
  gen tt10 = 10
}


frame meansurv {
  replace female = 0
  predict S10M, surv timevar(tt10) merge
}


frame meansurv {
  replace female = 1
  predict S10F, surv timevar(tt10) merge
}


standsurv , surv ci timevar(tt)    ///
            atvar(S10M_ss S10F_ss) ///
            at1(female 0)          ///
            at2(female 1) 

list tt S10F_ss  S10M_ss if  tt==10
frame meansurv: summ S10F S10M        
            
            
// stpm3km is useful to assess how well the model fit to data.
// This could be the data used to develop the model or it 
// could be data external to the study.
// The basic idea is to obtain predicted marginal survival in subgroups
// and compare these to non-parametric Kaplan-Meier estimates.

// This can be done using stpm3km.
// For example we can form 6 equally sized groups based on the 
// prognostic index.

stpm3km, groups(5)

// We can also look in other groups
stpm3km age, groups(5)

     // We can also look in other groups
stpm3km stage, factor