\subsection{Stata code for obtaining predictions}\label{code}

This appendix outlines Stata code used to obtain predictions presented in the paper. Some user-defined Stata commands are required which can be installed from the Boston College Statistical Software Components (SSC) archive by calling,

\begin{stsyntax}
{ssc install} \LB\texttt{\it command}\RB
\end{stsyntax}

The following must be installed before running the code:
\begin{itemize}
\item \texttt{stpm2}: To fit the flexible parametric models described in Section \ref{fpm}.
\item \texttt{rcsgen}: To generate the restricted cubic spline functions.
\item \texttt{stcrprep}: To restructure data and calculate time-dependent censoring weights in order to fit models on the subdistribution hazards scale using standard Stata commands.
\end{itemize} 

To obtain marginal (and non-marginal) estimates using standardisation, the \texttt{standsurv} command must be installed. This will be released on SSC soon, however, in the meantime, it can be installed by running,

 \begin{stsyntax}
{net from https://www.pclambert.net/downloads/standsurv}
\end{stsyntax}

\subsubsection{Preparing the data for analysis}
To prepare the data for a survival analysis in Stata, we must first run the \texttt{stset} command. We identify the variable that records survival time (in days), \texttt{exit2}, the indicator variable for cause of death, \texttt{cod}, where death from \texttt{cancer = 1}, \texttt{CVD = 2} and \texttt{other causes = 3} and finally the variable for date of diagnosis, \texttt{dx}. The \texttt{scale} option is used to transform the survival time into years from days and we use the \texttt{exit} option to restrict follow-up time to 10 years from diagnosis and censor those still alive at 2014. In order to ensure that the death indicator, \texttt{\_d}, generated after \texttt{stset} matches the death indicator for cause of death, we create a new cause of death indicator, \texttt{cod2}, so that those who die either after 10 years from diagnosis or 2014 are administratively censored. Finally, to generate restricted cubic spline variables for the non-linear effect of age centred at 45 years old at diagnosis, we use \texttt{rcsgen}. For 3 degrees of freedom, 3 new age spline variables are created, \texttt{rcsage1 \(-\) rcsage3}, and we store knot positions and matrix for orthogonalization which are required for post-estimation predictions at specific ages.  

{\singlespacing
\begin{stlog}
\normalsize{\input{stset.log.tex}\nullskip}
\end{stlog}}
\doublespacing

To restructure the data and calculate the time-dependent censoring weights so that we may fit a model on the subdistribution hazards scale, we use \texttt{stcrprep}\cite{Lambert2016a}. Here, we specify \texttt{wtstpm2} to estimate the censoring distribution using a Royston-Parmar flexible parametric model with covariates included in the \texttt{censcov} option. The data is restructured based on the variable \texttt{failcode}, which splits the data according to the cause of interest. This is used to fit identify for which cause the model is to be fitted for. For clarity, we create dummy variables for each of the causes of death from \texttt{failcode} and generate \texttt{\_cancer}, \texttt{\_cvd} and \texttt{\_other}. Another indicator variable, \texttt{event}, is also created to identify at which split time interval, or row, death (from any cause) is observed for that patient. To incorporate the calculated weights from \texttt{stcrprep}, we must \texttt{stset} the data again with \texttt{tstart} and \texttt{tstop}. These are also provided by \texttt{stcrprep} and give the times at which an individual starts and stops being at risk.   

{\singlespacing
\begin{stlog}
\normalsize{\input{stcrprep.log.tex}\nullskip}
\end{stlog}}
\doublespacing


\subsubsection{Model}
The model described in Section \ref{fpm} can be fitted in two ways after preparing the data. We can either fit separate models for each of the causes of death, or fit a single model to cancer, CVD and other causes simultaneously. Here, we demonstrate for the latter to make illustration of the code for obtaining predictions post-estimation easier. However, in order to fit the equivalent single model with coefficients comparable to the models fitted individually to each of the causes of death, the knot locations on the cause-specific survival time distributions must be stored. These are stored in global macros for each of the causes of death.

{\singlespacing
\begin{stlog}
\normalsize{\input{model.log.tex}\nullskip}
\end{stlog}}
\doublespacing

Here we define a global macro of the list of covariates to be included in the single model. As the data is stacked, interactions need to be created between the covariates and the indicator variable for each cause of death. See Lunn and McNeil\cite{Lunn1995} for further details. The baseline coefficient, i.e. the constant in the cause-specific model, is calculated in \texttt{\_cancer}, \texttt{\_cvd} and \texttt{\_other}. We therefore fit a model for each of the causes of death simultaneously without a constant using \texttt{nocons} and the baseline splines using \texttt{rcsbaseoff}. Instead, the baseline splines are specified as time-dependent splines for the coefficient that corresponds to the constant in its respective model for that particular cause of death. These were stored in the global macro \texttt{bknotstvc\_opt}. Since knots are specified according to the time scale, rather than the log-time scale, the \texttt{knscale(time)} option is used. 

{\singlespacing
\begin{stlog}
\normalsize{\input{model2.log.tex}\nullskip}
\end{stlog}}
\doublespacing

\subsubsection{Predictions}
Although \texttt{standsurv} was written for obtaining marginalised predictions, it can also be used to obtain non-marginalised estimates. This is done by simply specifying the entire covariate pattern so that the predictions are not averaged over any covariate distribution. To obtain predictions at a specific age, we need to calculate the spline variables at that particular age centred at 45 years old with the same knot locations and projection matrix as before. The spline variables are stored in the local macros \texttt{c1}, \texttt{c2} and \texttt{c3}. An example is given below when the cause of interest is cancer and we want to make comparisons between the most and least deprived male patients aged either 50, 65, or 80 years old at diagnosis. 

{\singlespacing
\begin{stlog}
\normalsize{\input{cov_age.log.tex}\nullskip}
\end{stlog}}
\doublespacing

As we do not average over each observation, we must tell \texttt{standsurv} to only take the first observation in the stacked data to calculate non-marginalised predictions. This is done using \texttt{if \_n == 1}. The \texttt{failure} option is used to obtain the cumulative incidence functions that is specified in each \texttt{at} option. To calculate the difference between \texttt{at1} and \texttt{at2}, we use \texttt{contrast(difference)}. 

{\singlespacing
\begin{stlog}
\normalsize{\input{cif_age.log.tex}\nullskip}
\end{stlog}}
\doublespacing

Since we are making predictions at particular covariate patterns for each of the causes separately, specifying \texttt{rmft} gives us estimates of the expected life-years lost due to a particular cause of death. To calculate RMLT, we need to take the sum of all of the \texttt{at} options, where the expected life-years lost due to cancer, CVD and other causes is specified in each. We do this by creating our own contrast in a user-defined mata function which can be called in the option \texttt{userfunction}. An example of this is also given below.    

{\singlespacing
\begin{stlog}
\normalsize{\input{rmft_age.log.tex}\nullskip}
\end{stlog}}
\doublespacing

In order to obtain marginalised estimates, in each \texttt{at} option, only the covariate pattern for the group of interest need to be given. For the covariate distribution that we want to average over, as we have created interactions between the covariates and the causes of death, these must be mapped to each covariate e.g. \texttt{sex\_cancer = sex}. The others are excluded from the \texttt{at} option for the other causes of death. In this case, because we want to average over covariates that we wish to standardise by, we need to identify the row for each patient in the stacked data that corresponds to the failure time of that individual. This is done by creating the indicator variable \texttt{first} and using it as an \texttt{if} condition in \texttt{standsurv}. As before, we give an example for specifying macros for use in the \texttt{at} options for deaths due to cancer. 

{\singlespacing
\begin{stlog}
\normalsize{\input{cov_stand2.log.tex}\nullskip}
\end{stlog}}
\doublespacing

The cause-specific CIF differences are thus calculated as follows,

{\singlespacing
\begin{stlog}
\normalsize{\input{cif_stand2.log.tex}\nullskip}
\end{stlog}}
\doublespacing

As highlighted above, we can write user-functions to define our own contrasts. Below is an example for when interest is in calculating the difference in RMLT between the most and least deprived patients. 

{\singlespacing
\begin{stlog}
\normalsize{\input{rmft_stand2.log.tex}\nullskip}
\end{stlog}}
\doublespacing