A natural history and copula-based joint model for regional and distant breast cancer metastasis

The few existing statistical models of breast cancer recurrence and progression to distant metastasis are predominantly based on multi-state modelling. While useful for summarising the risk of recurrence, these provide limited insight into the underlying biological mechanisms and have limited use for understanding the implications of population-level interventions. We develop an alternative, novel, and parsimonious approach for modelling latent tumour growth and spread to local and distant metastasis, based on a natural history model with biologically inspired components. We include marginal sub-models for local and distant breast cancer metastasis, jointly modelled using a copula function. Different formulations (and correlation shapes) are allowed, thus we can incorporate and directly model the correlation between local and distant metastasis flexibly and efficiently. Submodels for the latent cancer growth, the detection process, and screening sensitivity, together with random effects to account for between-patients heterogeneity, are included. Although relying on several parametric assumptions, the joint copula model can be useful for understanding – potentially latent – disease dynamics, obtaining patient-specific, model-based predictions, and studying interventions at a population level, for example, using microsimulation. We illustrate this approach using data from a Swedish population-based case-control study of postmenopausal breast cancer, including examples of useful model-based predictions.

The tumour volumes at symptomatic detection are simulated using Equation 2 from the manuscript, assuming the parameter η ′ = − log(η) = 9.
The bivariate outcome of affected lymph nodes at diagnosis and time to distant metastasis was simulated by assuming a bivariate Frank copula with θ = −3; intuitively, this corresponds to a negative correlation between the two outcomes. We drew random, correlated bivariate samples from the copula function and consequently applied the inversion method to each marginal model to simulate the observed outcomes; here we assumed the parameter values log(γ 1 ) = −0.5, log(γ 2 ) = 10, log(ω 1 ) = 0, and log(ω 2 ) = 15. Further details on the simulation algorithm for the time to detection of distant metastasis are included elsewhere (Gasparini and Humphreys, 2022).
Finally, we discretized time in years and we applied administrative censoring after 20 years of followup. Discretization of time in years can be seen as a worst-case scenario; in practice, discretizing time in finer intervals (e.g. months) should improve the performance of the algorithm.

A.1.4 Methods
We fitted the true data-generating model for data collected in the absence of screening, as described in the main body of the manuscript, to each independently simulated data set. Standard errors for the model parameters were computed by inverting the Hessian matrix at the optimum.

A.1.5 Performance Measures and Number of Replications
The key performance measure of interest is bias, to quantify whether the estimation procedure can retrieve the true data-generating parameters, on average; in order to quantify the uncertainty in the bias estimation, we report Monte Carlo standard errors as well. We also report empirical and model-based standard errors, including relative percentage errors in standard errors, and coverage probability. Performance measures are described in more detail elsewhere (Morris et al., 2019).
To estimate the key performance measure with a given precision, we first ran 10 replications of the simulation study. Allowing a Monte Carlo standard error for bias of 0.01, we would require n sim = Var/MCSE 2 ≈ 415 repetitions of our full simulation study. This calculation was based on assuming Var = 0.041542, which was taken as the largest value among all empirical and modelbased standard errors for Ψ obtained from the above-mentioned 10 replications. To be conservative we rounded up the number of repetitions to the nearest 100, which yielded a total of n sim = 500. Figure A1: Bias with 95% confidence intervals based on the Monte Carlo standard errors, simulation study for the model in the absence of screening.

A.1.6 Results
Bias with 95% confidence intervals (based on Monte Carlo standard errors) for each model parameter is depicted in Figure A1. All parameters could be estimated with virtually no bias, except for log(ω 1 ) and log(ω 2 ) (parameters of the marginal model for time to detection of distant metastasis) where we could see a small positive bias in absolute terms. Nevertheless, we deemed this bias to be negligible as the relative bias on the original scale of the parameters was close to zero ( Figure A2).
It can also be shown that estimation of the model standard errors using the inverse of the Hessian matrix at the optimum yielded satisfactory results, as highlighted by empirical and model-based standard errors being close to each other ( Figure A3) and by the small relative percentage error in standard errors ( Figure A4).
Finally, we could show that the small, aforementioned bias did not affect coverage probability, which was (overall) close to the optimal value of 95% (Figures A5).
Estimated values for bias, empirical and model-based standard errors, relative percentage errors, and coverage probability, together with their Monte Carlo standard errors, are tabulated for reference in Table A1.
Model-Based S.E. Figure A4: Relative percentage error in standard error with 95% confidence intervals based on the Monte Carlo standard errors, simulation study for the model in the absence of screening.

A.2.1 Aim
This simulation study aims to assess the ability to retrieve the true data-generating parameters of the model for data collected in a screened population.

A.2.2 Data-Generating Mechanisms
The data-generating mechanism for this simulation study is closely related to that described in Section A.1. We, once again, simulated independent datasets of 1500 individuals each, assuming spherical tumours that originate from a single cell of volume V Cell with d Cell = 0.01 mm and that such tumours are detectable only once they cross the threshold of volume V 0 (corresponding to a diameter d 0 = 0.5 mm). We again assumed the parameters k N = k W = 4 to be fixed.
Compared to the previous simulation study, we do not require the constraint τ 1 = τ 2 (for identifiability purposes) anymore. We, therefore, simulate inverse growth rates from a Gamma distribution with shape parameter τ 1 = 2 and rate τ 2 = 4; this leads to a distribution with a mean of 0.5 and a variance of 0.125.
Age at onset of breast cancer was simulated assuming age-increasing rates based on a procedure described previously (Abrahamsson and Humphreys, 2016;, which use published breast cancer incidence rates (Forastero et al., 2010). We then enforced a screening program starting at the age of 40, with screens every two years. Death was not included in our simulations, analogously as with previous work, see e.g. Isheden et al.
(2019); Gasparini and Humphreys (2022). Tumours were detected according to the logistic screening sensitivity function described in the manuscript, including the diameter of the tumour as a covariate and assuming parameters β 1 = −5 (model intercept) and β 2 = 0.5 (coefficient of size).
The remaining details of the data-generating mechanisms (including assumed parameters) are analogous to what was previously described in Section A.1.

A.2.4 Methods
The method of interest is the true data-generating model for a screened population, which we fit to each independently simulated data set; once again, standard errors for the model parameters were computed by inverting the Hessian matrix at the optimum.

A.2.5 Performance Measures and Number of Replications
The key performance measures of interest are the same ones that were used in the previous simulation study.
To estimate the key performance measure with a given precision, we once again start by running 10 replications of the simulation study. Allowing a Monte Carlo standard error for bias of 0.01, we would require n sim = Var/MCSE 2 ≈ 746 repetitions of our full simulation study. This calculation was based on assuming Var = 0.07464286, which was taken as the largest value among all empirical and model-based standard errors for Ψ obtained from the above-mentioned 10 replications. To be conservative we rounded up the number of repetitions to the nearest 100, which yielded a total of n sim = 800.

A.2.6 Results
Our model implementation converged for 790 repetitions out of 800 (98.75%); performance measures for the iterations that converged are presented below and tabulated in Table A2.

Figures A6 and A7 describe bias and relative bias for each model parameter. Most parameters could
be estimated with no significant absolute bias, or with negligible bias on a relative scale (such as log(τ 2 ), − log(η), and log(ω 2 )). We could however observe significant bias for the parameters of the screening sensitivity model, with overestimation of the intercept β 1 and underestimation of the coefficient for size β 2 : this is in agreement with previous work, see e.g. Gasparini and Humphreys (2022), and not surprising given the small number of tumours of very small size at detection.
Despite this significant bias, the fitted sensitivity lines did not differ too much from the truth, as illustrated in Figure A8.
Our approach for the estimation of model standard errors was deemed adequate, as the modelbased standard errors closely matched empirical standard errors ( Figure A9); consequently, relative Figure A7: Relative bias with inter-quartile intervals, for parameters on the original scale. Simulation study for the model in a screened population.
Finally, coverage probabilities were optimal for all model parameters except for those with the largest absolute and relative biases, namely the coefficients of the screening sensitivity function β 1 and β 2 , for which coverage probability was low at 37.8% and 63.5%, respectively. Figure A10: Relative percentage error in standard error with 95% confidence intervals based on the Monte Carlo standard errors, simulation study for the model in a screened population.

Events Censored
At risk Overall Figure B2 Marginal survival curves (in black) for time to diagnosis of distant metastasis for subjects with zero, one, and two affected lymph nodes at diagnosis. Curves are shown for all models compared in the analysis of the CAHRES data, assuming a Frank, Clayton, AMH, or product copula, and averaged over the observed covariates' distribution in each subgroup. Standard errors are computed using the numerical delta method as implemented in the predictnl function from the rstpm2 package (Liu et al., 2018). Kaplan-Meier curves (in grey) for relevant subgroups of the study population are also included for comparison.

Table B2
Estimated cured proportions under different copula formulations; application to CAHRES data.
Values in each cell are cured proportions with 95% confidence intervals, and are presented overall and by number of affected lymph nodes. Standard errors are computed using the numerical delta method as implemented in the predictnl function from the rstpm2 package (Liu et al., 2018).

Table B3
Estimated survival probabilities for time to distant metastasis according to all models with different copula formulations; application to CAHRES data. Values in each cell are survival probabilities with 95% confidence intervals. Standard errors are computed using the numerical delta method as implemented in the predictnl function from the rstpm2 package (Liu et al., 2018).