*! version 0.1 22Dec2020
also see: stpm2, strcs

Title

mrsprep -- Prepare data for fitting a marginal relative survival model

Syntax

mrsprep using filename[if] [in], agediag(varname) datediag(varname) breaks(numlist) [options]

Options
agediag(varname) age at diagnosis (in years)

breaks(numlist) break points to calculate time-dependent weights

datediag(varname) date at diagnosis (in years)

by(varlist) calculate mean hazard and weights separately by varlist

indweights(varname) individual level weights

keep(varlist) name of additional variables to keep in dataset

newframe(#) name of new frame to create

pmage(varname) name of age variable in popmort file

pmother(varname) name of other variables in popmort file

pmrate(varname) name of rate variable in popmort file

pmyear(varname) name of calendar year variable in popmort file

pmmaxage(#) maximum age in popmort file

pmmaxyear(#) maximum year in popmort file

verbose more detailed output

Description

mrsprep prepares data for fitting modelling marginal relative survival models using parametric models.
mrsprep obtains the mean expected hazard rate at each event time (weighted by the inverse of expected survival)
and expands the data so that time-dependent weights can be incorporated when subsequently fitting models.

By expanding the data and calculating weighted mean expected mortality rates and time-dependent weights,
it is then possible to use standard relative survival estimation commands (e.g. stpm2 and strcs)
to directly estimate marginal relative survival.
Including individual weights enables external (age) standardization without the need to model the effect of age.

The data must be stset before using mrsprep. As the data is expanded the id() option of stset must be used. There should only be 1 row of data per subject before using mrsprep.

After using using mrsprep a new data frame is created (default name, mrs_data). The original data is still available. It possible to switch between the data sets (e.g. frame change default).

Options

agediag(varname) gives the name of the variable containing age at diagnosis. This should be in years. Note that if possible it is best to
avoid using truncated (integer) age as this assumes that each person was diagnosed on their birthday.

breaks(numlist) break points to split the time-scale for calculation of the time-dependent weights. The narrower the interval, the greater the accuracy.

datediag(varname) gives the name of the variable containing the date at diagnosis.

by(varlist) will calculate mean expected hazard rates separately by variables in varlist.

indweights(varname) incorporates individual level weights to up- or down-weight individuals relative to a reference
population. This is useful for external age standardization.

keep(varlist) keeps the named variables in the expanded dataset.

newframe(framename) gives the name of new frame to create. mrsprep creates a new frame mrs_frame which is the working frame immediately after running the command. Using frame change default will return to the original data frame. If the frame mrs_frame exists and you want to replace it use newframe(,replace).

pmage(varname) gives the name of the age variable in the population mortality file. The default is _age. This variable cannot exist in
the patient data file, but should exist in the population mortality file.

pmother(varname) names additional variables in the population mortality file. Usually this will include sex, but could
additionally be, for example, information on region, deprivation etc. All variables listed should be in both the data and
the population mortality file.

pmrate(varname) gives the name of the rate variable in the population mortality file. The default is rate. The rate should be expressed
per person year. If you only have one year survival probabilities in the population mortality file, then you can obtain
the rate using gen rate = -ln(survprob), where survprob is the one year survival probability.

pmmaxage(#) specifies the maximum age for which general-population mortality rates are provided in the population mortality file.
Rates for individuals older than this value are assumed to be the same as for maximum age #.
The default maximum age is 99.

pmmaxyear(#) specifies the maximum year for which general-population mortality rates are provided in the population mortality file.
Rates for individuals still at risk after this year are assumed to be the same as for maximum year #.

pmyear(varname) name of year variable in the population mortality file. The default is _year. This variable cannot exist in the patient data file, but should exist in the population mortality file.

verbose give some details about the how far the calculations have proceeded.

Examples

All examples use melanoma data available with strs. You can either run the examples and return to your current data or leave the example data in memory.
For the latter, you should not have data in memory in the active frame.

Example 1: Direct modelling of marginal relative survival

Load and stset melanoma data

. use "https://pclambert.net/data/melanoma.dta"

. stset exit, origin(dx) failure(status=1,2) id(id) exit(time dx + 10*365.24) scale(365.24)"{

Use mrsprep to expand data and calculate weights

. mrsprep using "https://pclambert.net/data/popmort.dta", pmother(sex) agediag(age) datediag(dx) ///

verbose breaks(0(0.2)10)

Incorprate weights when using stset

. stset tstop [iw=wt], enter(tstart) failure(event==1)

Fit marginal model using the weighted mean hazard

. stpm2, scale(hazard) df (5) bhazard(meanhazard_wt) vce(cluster id)

Predict marginal relative survival

range tt 0 10 101

predict s_mrs, surv timevar(tt)

Plot results

. twoway (line s_mrs* tt, lcolor(red..) lpattern(solid dash dash)) ///

, legend(off) ///

ylabel(0.6(0.1)1, format(%3.1f)) ///

ytitle("Marginal relative survival") ///

xtitle("Years from diagnosis")

(click to run Example 1)

(click to run Example 1 and leave data in memory)

Example 2: Incorporate individual weights for external standardization

Load and stset melanoma data

. use "https://pclambert.net/data/melanoma.dta"

. stset exit, origin(dx) failure(status=1,2) id(id) exit(time dx + 10*365.24) scale(365.24)"{}

Change age groups to those defined in ICSS

. drop agegrp

. egen agegrp=cut(age), at(0 45 55 65 75 200) icodes

. replace agegrp = agegrp + 1

. label variable agegrp "Age group"

. label define agegrplab 1 "0-44" 2 "45-54" 3 "55-64" 4 "65-74" 5 "75+", replace

. label values agegrp agegrplab

Create weights in reference population (ICSS)

. recode agegrp (1=0.28) (2=0.17) (3=0.21) (4=0.20) (5=0.14), gen(ICSSwt)

Calculate relative weights}

. local total= _N

. bysort agegrp: gen a_age = _N/`total'

. gen double wt_age = ICSSwt/a_age

Prepare data for marginal model

. mrsprep using "https://pclambert.net/data/popmort.dta", pmother(sex) agediag(age) datediag(dx) ///

verbose breaks(0(0.2)10) ///

indweights(wt_age) ///

newframe(mrs_stand, replace)

Incorprate weights when using stset

. stset tstop [iw=wt], enter(tstart) failure(event==1)

Fit marginal model using the weighted mean hazard

. stpm2, scale(hazard) df (5) bhazard(meanhazard_wt) vce(cluster id)

Predict externally age standardized marginal relative survival

. range tt 0 10 101

. predict s_mrs, surv timevar(tt) ci

Plot results

. twoway (line s_mrs* tt, lcolor(red..) lpattern(solid dash dash)) ///

, legend(off) ///

ylabel(0.6(0.1)1, format(%3.1f)) ///

ytitle(Marginal relative survival) ///

xtitle(Years from diagnosis)

(click to run Example 2)

(click to run Example 2 and leave data in memory)

Example 3: Fit a parametric model for the CIF

Load and stset melanoma data

. use "https://pclambert.net/data/melanoma.dta"

. stset exit, origin(dx) failure(status=1,2) id(id) exit(time dx + 10*365.24) scale(365.24)"{}

Change age groups to those defined in ICSS

. drop agegrp

. egen agegrp=cut(age), at(0 45 55 65 75 200) icodes

. replace agegrp = agegrp + 1

. label variable agegrp "Age group"

. label define agegrplab 1 "0-44" 2 "45-54" 3 "55-64" 4 "65-74" 5 "75+", replace

. label values agegrp agegrplab

Create weights in reference population (ICSS)

. recode agegrp (1=0.28) (2=0.17) (3=0.21) (4=0.20) (5=0.14), gen(ICSSwt){p_end

Proportion within each age group by sex to calculate weights

. gen female = sex == 2

. bysort female: egen totalsex = total(sex)

. bysort agegrp female: gen a_age_sex = _N/totalsex

. gen double wt_age_sex = ICSSwt/a_age_sex

Prepare data for marginal model

. mrsprep using popmort.dta, pmother(sex) agediag(age) datediag(dx) ///

pmmaxyear(2000) ///

verbose breaks(0(0.2)10) ///

indweights(wt_age_sex) ///

by(female)

Incorprate weights when using stset

. stset tstop [iw=wt], enter(tstart) failure(event==1)

Fit proportional hazards marginal model

. stpm2 female, scale(hazard) df (5) bhazard(meanhazard_wt) vce(cluster id)

Relax proportional hazards assumption

. stpm2 female, scale(hazard) df (5) bhazard(meanhazard_wt) vce(cluster id) ///

tvc(female) dftvc(3)

predict externally age standardized marginal relative survival by sex

. range tt 0 10 101

. predict s_mrs_male, surv timevar(tt) ci at(female 0)

. predict s_mrs_female, surv timevar(tt) ci at(female 1)

Plot results

. twoway (line s_mrs_male* tt, lcolor(red..) lpattern(solid dash dash)) ///

(line s_mrs_female* tt, lcolor(blue..) lpattern(solid dash dash)) ///

, legend(off) ///

ylabel(0.6(0.1)1, format(%3.1f)) ///

ytitle(Marginal relative survival) ///

xtitle(Years from diagnosis)

(click to run Example 3)

(click to run Example 3 and leave data in memory)

Author

Paul Lambert, University of Leicester, UK.
(paul.lambert@leicester.ac.uk)

References

Lambert PC, Syriopoulou E, Rutherford M.J.

Also see

stpm2 strcs;
[ST] stset