//==================// // EXERCISE 288 // New June 2019 //==================// // Using standsurv for life expectancy // (a) Load melanoma data, restrict to the 1985-1994 calendar period, // those aged over 18 and merge in the expected mortality rates. use melanoma, clear keep if year8594 == 1 & age>=18 stset surv_mm, failure(status=1 2) scale(12) exit(time 120.5) id(id) gen _age = min(int(age + _t),99) /*merge on expected rates at exittime*/ gen _year = int(yydx + _t) sort _year sex _age merge m:1 _year sex _age using popmort, keep(match master) keepusing(rate) drop _age _year _merge // (b) Fit a flexible parametric model including continuous age using // restricted cubic splines with 3 knots and and sex. Allow the effect // of age and sex to have time-dependent effects. rcsgen age, df(3) gen(agercs) global ageknots `r(knots)' gen female = sex==2 // create dummy variable for females stpm2 agercs* female, scale(hazard) df(5) bhazard(rate) /// tvc(agercs* female) dftvc(3) // (c) We will do some predictions for a 70 year old male // Store the values of the age spline variables for a 70 year old rcsgen, scalar(70) gen(a) knots(${ageknots}) // (i) Predict relative survival for a 70 year old male up to 10 years and plot range tt 0 10 101 predict rs70, at(agercs1 `=a1' agercs2 `=a2' agercs3 `=a3' female 0) /// surv timevar(tt) ci twoway (line rs70* tt, lcolor(red..) lpattern(solid dash dash)), /// ylabel(0(0.1)1, format(%3.1f) angle(h)) /// legend(off) /// ytitle("Relative Survival") /// xtitle("Years from diagnosis") /// ylabel(0(0.1)1) // (ii) The area under the survival curve gives the restricted mean relative survival time // We can approximate this using the integ command integ rs70 tt // (iii) Alternatively we can use the \stcmd{rmst} option of \stcmd{stpm2}'s \stcmd{predict} command. gen t10 = 10 in 1 predict rs70b in 1, at(agercs1 `=a1' agercs2 `=a2' agercs3 `=a3' female 0) zeros rmst tmax(10) list rs70b* in 1, noobs // Interpret the restricted mean relative survival. // This is rather an awkward interpretation. "If it was not possible to die // from other causes, average life expectency over the next 10 years is // 7.62 years. // (d) Note that we have not extrapolated beyond the range of follow-up, but as this // is a parametric model we can extrapolate if we want to. // We will extrapolate to 40 years. This is when someone 70 at diagnosis // would be 110, so a close to zero probability of being alive. range tt_long 0 40 101 predict rs70_long, at(agercs1 `=a1' agercs2 `=a2' agercs3 `=a3' female 0) surv /// timevar(tt_long) ci twoway (line rs70_long* tt_long, lcolor(red..) lpattern(solid dash dash)), /// ylabel(0(0.1)1, format(%3.1f) angle(h)) /// legend(off) /// ytitle("Relative Survival") /// xtitle("Years from diagnosis") /// ylabel(0(0.1)1) // Explain why this is really not a very interesting extrapolation. // This is not a useful extrapolation we are in a hypothetical world. It is // not really of interest to look at the 40 year survival probability for a // 70 year old in a world where it is not possible to die from other causes. // (e) We will now combine relative survival with expected survival so we // can estimate the all cause probability of death for a 70 year old man. // We will switch to standsurv for these predictions. We are interested // in the prediction of a single individual age 70, so we restict the // prediction to the first row and feed in the relevant covariate values // using the at1() option. gen age70 = 70 in 1 gen dx70 = mdy(1,1,1990) standsurv if _n==1, /// at1(agercs1 `=a1' agercs2 `=a2' agercs3 `=a3' female 0) /// timevar(tt) ci /// atvar(s70_all) /// expsurv(using(popmort.dta) /// popmort file. agediag(age70) /// age at diagnosis datediag(dx70) /// date at diagnosis pmage(_age) /// age variable in popmort file pmyear(_year) /// year variable in popmort file pmother(sex) /// other variables in popmort file pmrate(rate) /// rate varible in popmort file pmmaxyear(2000) /// maximum year in popmort file at1(sex 1) /// variables to match with main at options expsurvvar(expsurv70)) // output expected survival twoway (line s70_all tt) /// (line rs70 tt) /// (line expsurv70 tt) /// , ylabel(0(0.1)1, format(%3.1f) angle(h)) /// legend(order(1 "All cause" 2 "Relative" 3 "Expected") cols(1) ring(0) pos(7)) /// ytitle("Survival") /// xtitle("Years from diagnosis") // Why is the all-cause survival decreasing when the relative survival is flat? // This is because the all-cause survival depends on on expected survival, which // is clearly (and obviously) decreasing over time. // Generate the resticted mean all-cause and expected survival time using // integ and interpret integ s70_all tt integ expsurv70 tt // "In the 10 years after diagnosis the expected year alive for a 70 year old man // diagnosed with Melanoma is 6.03 years. This can be contrasted with the // expected survival of a man in the general population of 7.67 years." // (f) Now repeat, but use standsurv. This will perform a more accuarte numerical // integration and also give 95% confidence intervals. standsurv if _n==1, /// at1(agercs1 `=a1' agercs2 `=a2' agercs3 `=a3' female 0) /// timevar(t10) rmst ci /// atvar(rmst10) /// expsurv(using(popmort.dta) /// popmort file. agediag(age70) /// age at diagnosis datediag(dx70) /// date at diagnosis pmage(_age) /// age variable in popmort file pmyear(_year) /// year variable in popmort file pmother(sex) /// other variables in popmort file pmrate(rate) /// rate varible in popmort file pmmaxyear(2000) /// maximum year in popmort file at1(sex 1) /// expsurvvar(exp_rmst)) list rmst10* in 1, noobs list exp_rmst* in 1, noobs // (g) If we use tt_long for our time variable then we will extrapolate both // relative and expected survival as a way to extrapolate all-cause survival. // In most situations this will work better than direct extrapolation of // all cause survival. standsurv if _n==1, /// at1(agercs1 `=a1' agercs2 `=a2' agercs3 `=a3' female 0) /// timevar(tt_long) ci /// atvar(s70_long) /// expsurv(using(popmort.dta) /// popmort file. agediag(age70) /// age at diagnosis datediag(dx70) /// date at diagnosis pmage(_age) /// age variable in popmort file pmyear(_year) /// year variable in popmort file pmother(sex) /// other variables in popmort file pmrate(rate) /// rate varible in popmort file pmmaxyear(2000) /// maximum year in popmort file at1(sex 1) /// expsurvvar(expsurv70_long)) twoway (line s70_long* tt_long, lcolor(blue..) lpattern(solid dash dash)) /// (line rs70_long tt_long) /// (line expsurv70_long tt_long) /// ,ylabel(0(0.1)1, format(%3.1f) angle(h)) /// legend(order(1 "All cause" 4 "Relative" 5 "Expected") cols(1) ring(0) pos(1)) /// ytitle("Survival") /// xtitle("Years from diagnosis") /// xline(10, lpatter(dash)) // Why can we now integrate the all cause survival curve to obtain an estimate // of life expectancy? Perform the integration using -integ- and interpret. integ s70_long tt_long integ expsurv70_long tt_long // The life expectency of a 70 year old man in the general population is 11.32 // years. For a 70 year old mean diagnosed with Melanoma, the life expectency is // 8.48 years. The loss in expectation of life is thus 2.85 years. // (h) Rather than use -integ- we can do the estimation using standsurv with // the rmst option. Note this is strictly a resticted mean, but we are // estimating up to a time point when we know the all-cause survival will // be effectively zero. gen t40 = 40 in 1 standsurv if _n==1, /// at1(agercs1 `=a1' agercs2 `=a2' agercs3 `=a3' female 0) /// timevar(t40) rmst ci /// atvar(lifeexp70) /// expsurv(using(popmort.dta) /// popmort file. agediag(age70) /// age at diagnosis datediag(dx70) /// date at diagnosis pmage(_age) /// age variable in popmort file pmyear(_year) /// year variable in popmort file pmother(sex) /// other variables in popmort file pmrate(rate) /// rate varible in popmort file pmmaxyear(2000) /// maximum year in popmort file at1(sex 1) /// expsurvvar(exp_lifeexp70)) list lifeexp70* exp_lifeexp70 in 1, noobs // (i) So far we have just estimated for a 70 year old. We can use standsurv // to obtain standardized survival curves over all (or a subset) of // individuals. We need to use the age at diagnosis and date of diagnosis // of each individual in the data set. // The youngest person in our analysis is 18 and so extrapolation for just 40 // years would only make them 58, so we create a new variable t100 with 100 years // of follow-up. This is longer than necessary for a 80 year old, but is needed // as we are now including much younger individuals. gen t100 = 100 in 1 standsurv, /// at1(agercs1 `=a1' agercs2 `=a2' agercs3 `=a3' female 0, atif(female==0)) /// timevar(t100) rmst ci /// atvar(lifeexp_stand) /// expsurv(using(popmort.dta) /// popmort file. agediag(age) /// age at diagnosis datediag(dx) /// date at diagnosis pmage(_age) /// age variable in popmort file pmyear(_year) /// year variable in popmort file pmother(sex) /// other variables in popmort file pmrate(rate) /// rate varible in popmort file pmmaxyear(2000) /// maximum year in popmort file at1(sex 1) /// expsurvvar(exp_lifeexp_stand)) list lifeexp_stand* exp_lifeexp_stand in 1, noobs // The average life expectency for men diagnosed with Melanoma is 15.82 years. // This compares with 21.96 years for a similar group in the general population. // The average loss in expectation of life is 21.96-15.82 = 6.14 years. // (j) We need to remember the above is an average over all individuals. // This is fine as a summary measure, but hides the fact that life expectency // is highly dependent on age. Thus it can be useful to look at life expetency // over a range of ages. // We will predict life expectency and the loss in expectation of life for // men aged 40, 50, 60, 70 and 80. We just need to extract the values of the // age spline variables at these ages. foreach age in 40 50 60 70 80 { qui rcsgen, scalar(`age') gen(a`age'_) knots(${ageknots}) local atopt at1(agercs1 `=a`age'_1' agercs2 `=a`age'_2' agercs3 `=a`age'_3' female 0) gen tmpage`age' = `age' in 1 gen tmpdx`age' = mdy(1,1,1990) in 1 standsurv if _n==1, /// `atopt' /// list of at options timevar(t100) rmst ci /// atvar(le`age') /// new variable names expsurv(using(popmort.dta) /// popmort file. agediag(tmpage`age') /// age at diagnosis datediag(tmpdx`age') /// date at diagnosis pmage(_age) /// age variable in popmort file pmyear(_year) /// year variable in popmort file pmother(sex) /// other variables in popmort file pmrate(rate) /// rate varible in popmort file pmmaxyear(2000) /// maximum year in popmort file at1(sex 1) /// expsurvvar(exp`age')) drop tmpage`age' tmpdx`age' } // Display estimates for selected ages foreach age in 40 50 60 70 80 { di "Age: `age', Cancer = " %5.3f le`age'[1] ", Expected = " %5.2f exp`age' } // (k) In question 284 we calculated the average loss in expectation of life // if males had the relative survival of females. This can also been done // using -standsurv-. We need two at options, one where we predict life // expectancy for males and one where we predict life expectancy for // males, if they had the excess mortality rates of females. // Note we restrict the standardisation to males. As we have no interaction // between age and sex, the at option is simple. See question 287 on how // to code when there are interactions between covariates. // Also note how we need to specify the at suboptions within the expsurv option. standsurv, /// at1(female 0, atif(female==0)) /// prediction for males at2(female 1, atif(female==0)) /// prediction for males (female rs) timevar(t100) rmst ci /// atvar(lifeexp_m lifeexp_m_withf) /// expsurv(using(popmort.dta) /// popmort file. agediag(age) /// age at diagnosis datediag(dx) /// date at diagnosis pmage(_age) /// age variable in popmort file pmyear(_year) /// year variable in popmort file pmother(sex) /// other variables in popmort file pmrate(rate) /// rate varible in popmort file pmmaxyear(2000) /// maximum year in popmort file at1(sex 1) /// expected rates for males at2(sex 1) /// expected rates for males expsurvvar(exp_m1 exp_m2)) /// contrast(difference) /// contrastvar(lel_stand) list lifeexp_m lifeexp_m_lci lifeexp_m_uci in 1 list lifeexp_m_withf* in 1 list lel_stand* in 1 // The average gain in life expectency for males if they had female excess // mortality rates is 2.11 years.