# Multiple imputation of missing covariates when using the Fine–Gray model

Edouard F. Bonneville  Department of Biomedical Data Sciences, Leiden
University Medical Center, the Netherlands Jan Beyersmann  Institute of
Statistics, Ulm University, Germany  Ruth H. Keogh  Department of Medical
Statistics, London School of Hygiene and Tropical Medicine, United Kingdom
Jonathan W. Bartlett  Department of Medical Statistics, London School of
Hygiene and Tropical Medicine, United Kingdom  Tim P. Morris  MRC Clinical
Trials Unit at UCL, United Kingdom  Nicola Polverelli  Unit of Bone Marrow
Transplantation, Division of Hematology, Fondazione IRCCS Policlinico San
Matteo di Pavia, Italy  Liesbeth C. de Wreede  Department of Biomedical Data
Sciences, Leiden University Medical Center, the Netherlands DKMS Clinical
Trials Unit, Germany  Shared senior authorship  Hein Putter  Department of
Biomedical Data Sciences, Leiden University Medical Center, the Netherlands
Mathematical Institute, Leiden University, the Netherlands  Shared senior
authorship

###### Abstract

The Fine–Gray model for the subdistribution hazard is commonly used for
estimating associations between covariates and competing risks outcomes. When
there are missing values in the covariates included in a given model,
researchers may wish to multiply impute them. Assuming interest lies in
estimating the risk of only one of the competing events, this paper develops a
substantive-model-compatible multiple imputation approach that exploits the
parallels between the Fine–Gray model and the standard (single-event) Cox
model. In the presence of right-censoring, this involves first imputing the
potential censoring times for those failing from competing events, and
thereafter imputing the missing covariates by leveraging methodology
previously developed for the Cox model in the setting without competing risks.
In a simulation study, we compared the proposed approach to alternative
methods, such as imputing compatibly with cause-specific Cox models. The
proposed method performed well (in terms of estimation of both subdistribution
log hazard ratios and cumulative incidences) when data were generated assuming
proportional subdistribution hazards, and performed satisfactorily when this
assumption was not satisfied. The gain in efficiency compared to a complete-
case analysis was demonstrated in both the simulation study and in an applied
data example on competing outcomes following an allogeneic stem cell
transplantation. For individual-specific cumulative incidence estimation,
assuming proportionality on the correct scale at the analysis phase appears to
be more important than correctly specifying the imputation procedure used to
impute the missing covariates.

Keywords: Fine–Gray model; subdistribution hazard; competing risks; multiple
imputation; missing covariates; cumulative incidence function

##  1 Introduction

The presence of missing covariate data continues to be pervasive across
biomedical research. Among the many existing approaches for dealing with
missing covariate data, multiple imputation (MI) methods in particular have
become increasingly popular in practice 1. Compared to a complete-case
analysis (CCA), MI can provide inferences that are both less biased and more
efficient, under certain missingness mechanisms 2.

The most common approach to MI is to specify and fit regression models for
partially observed covariates, from which imputations are then generated.
Ideally, each one of these imputation models should be compatible with the
substantive model of interest. That is, the assumptions made by both models
should not conflict with each other, e.g. the imputation model should at least
feature the remaining substantive model covariates, as well as the outcome. We
refer to an imputation model as being “directly specified” when substantive
model covariates and outcome variable(s), or any transformations thereof, are
included explicitly as predictors in the imputation model. The use of these
directly specified imputation models in settings with missingness spanning
multiple covariates is more commonly known as MICE (multivariate imputation by
chained equations 3).

In the context of cause-specific Cox proportional hazards models 4, it has
been shown that the imputation model for a partially observed covariate should
at least include as predictors the other covariates from the substantive
model, together with the cause-specific cumulative hazard and event indicator
for each competing risk 5. Analogously to the standard single-event survival
setting, this directly specified imputation model is generally only
approximately compatible with the proportional hazards substantive model6.
Concretely, when the outcome model assumes proportional hazards, the
conditional distribution of a partially observed covariate modelled using MICE
is only an approximation of the “true” (i.e. implied assuming the substantive
model is correctly specified) conditional distribution of the partially
observed covariate given the outcome and other substantive model covariates.
If imputed values can instead be directly sampled from the latter
distribution, it would ensure compatibility between analysis and imputation
model. This alternative “indirect” way of obtaining imputations is referred to
as the substantive-model-compatible imputation (SMC-FCS7) approach, and it was
extended by Bartlett and Taylor to accommodate cause-specific Cox substantive
models8. In terms of estimating cause-specific hazard ratios, simulation
studies have shown that the SMC-FCS approach tends to outperform MICE in cases
when the substantive model is correctly specified 5, 8.

When a Fine–Gray subdistribution hazard model 9 is the substantive model of
interest, there has to our knowledge been no research on how one should
specify an imputation model for a missing covariate 10. Nevertheless, MICE is
still being used in the presence of missing covariates when the substantive
model is a Fine–Gray model, particularly in the context of prediction models.
While the structure of the imputation model is rarely reported, articles which
do describe their imputation procedure appear to use different approaches. For
example, in the prognostic Fine–Gray model presented by Archer et al. (where
the primary outcome was time to serious fall resulting in hospital admission
or death, with competing death due to other causes), the imputation model for
a missing covariate contained the other substantive model covariates, and the
cause-specific cumulative hazard and event indicator for each competing
risk11. In contrast, the MICE procedure reported as part of the prognostic
models presented by Clift et al. used cumulative subdistribution hazards in
the imputation model12. Heuristically, it would seem the latter approach is
more consistent with the substantive model, as the former imputes
approximately compatibly with a cause-specific Cox model structure rather than
the Fine–Gray model structure.

In this work, we extend the SMC-FCS approach for missing covariates to
accommodate a Fine–Gray substantive model for one of the competing events. In
the presence of random right censoring, the core idea is to multiply impute
the potential censoring times for individuals failing from competing events in
a first step 13, and thereafter use existing SMC-FCS methodology (originally
developed for the standard Cox model7) to impute the missing covariates in a
second step.

The structure of the manuscript is as follows. We introduce competing risks
notation in Section 2. In Section 3 we outline the proposed method, and
thereafter assess its performance in a simulation study in Section 4. We also
provide an illustrative analysis using a dataset from the field of allogeneic
hematopoietic stem cell transplantation (alloHCT) in Section 5. Finally,
findings are discussed in Section 6, together with recommendations on how to
impute covariates in competing risks settings more generally.

##  2 Notation

We consider a setting in which individuals experience only one of $K$
competing events. We denote the event time as $\tilde{T}$, and the competing
event indicator as $\tilde{D}\in\\{1,...,K\\}$. In practice, individuals are
subject to some right-censoring time $C$, meaning we only observe realisations
$(t_{i},d_{i})$ of $T=\min(C,\tilde{T})$ and $D=I(\tilde{T}\leq C)\tilde{D}$,
where $I(\cdot)$ is the indicator function and $D=0$ indicates a right-
censored observation. The cause-specific hazard for the
$k\textsuperscript{th}$ event is defined as

$h_{k}(t)=\lim_{\Delta t\downarrow 0}\frac{P(t\leq\tilde{T}<t+\Delta
t,\tilde{D}=k\,|\,\tilde{T}\geq t)}{\Delta t}.$

These hazards together make up the event-free survival function,

$P(\tilde{T}>t)=\exp\left\\{-\sum_{k=1}^{K}\int_{0}^{t}h_{k}(u)\mathop{}\\!\mathrm{d}u\right\\}=\exp\left\\{-\sum_{k=1}^{K}H_{k}(t)\right\\},$

assuming the distribution of $T$ is continuous, and $H_{k}(t)$ is the cause-
specific cumulative hazard function for the $k\textsuperscript{th}$ event. The
cause-specific cumulative incidence function is then defined as

$F_{k}(t)=P(\tilde{T}\leq
t,\tilde{D}=k)=\int_{0}^{t}h_{k}(u)S(u-)\mathop{}\\!\mathrm{d}u,$

where $S(u-)$ is the event-free survival probability just prior to $u$.

The subdistribution hazard for the $k\textsuperscript{th}$ event is defined as

$\displaystyle\lambda_{k}(t)$
$\displaystyle=\frac{-\mathop{}\\!\mathrm{d}\log\\{1-F_{k}(t)\\}}{\mathop{}\\!\mathrm{d}t},$
$\displaystyle=\frac{\mathop{}\\!\mathrm{d}F_{k}(t)}{\mathop{}\\!\mathrm{d}t}\times\\{1-F_{k}(t)\\}^{-1},$

which can be thought of as the hazard for the improper random variable
$\tilde{V}_{k}=I(\tilde{D}=k)\times\tilde{T}+I(\tilde{D}\neq k)\times\infty$,
for which we can write $F_{k}(t)=P(\tilde{V}_{k}\leq t)$ 14. The probability
mass at infinity makes $\tilde{V}_{k}$ improper, i.e. that its density
function does not integrate to one.

Suppose interest lies in modelling the cumulative incidence of one of the
competing events, say $D=1$, conditional on (time-fixed) covariates
$\mathbf{Z}$. The Fine–Gray model for cause 1 can be written as

$\lambda_{1}(t\,|\,\mathbf{Z})=\lambda_{0}(t)\exp(\boldsymbol{\beta}^{\intercal}\mathbf{Z}),$

with $\lambda_{0}(t)$ being the subdistribution baseline hazard function and
$\boldsymbol{\beta}$ representing the effects of covariates $\mathbf{Z}$ on
the subdistribution hazard. The cumulative incidence function for cause 1 can
then be written as

$F_{1}(t\,|\,\mathbf{Z})=1-\exp\Biggl{\\{}-\exp(\boldsymbol{\beta}^{\intercal}\mathbf{Z})\int_{0}^{t}\lambda_{0}(u)\mathop{}\\!\mathrm{d}u\Biggr{\\}},$

where $\int_{0}^{t}\lambda_{0}(u)\mathop{}\\!\mathrm{d}u=\Lambda_{0}(t)$ is
the cumulative baseline subdistribution hazard. If we define a baseline
cumulative incidence function $F_{0}(t)=1-\exp\\{-\Lambda_{0}(t)\\}$ (i.e. the
cumulative incidence when $\mathbf{Z}=0$), the model can also be written as

$F_{1}(t\,|\,\mathbf{Z})=1-\\{1-F_{0}(t)\\}^{\exp(\boldsymbol{\beta}^{\intercal}\mathbf{Z})}.$
(1)

In the presence of random right censoring, the Fine–Gray model is usually
fitted by maximising a partial likelihood that uses time-dependent inverse
probability of censoring weights (IPCW) 9.

##  3 MI approaches with a Fine–Gray substantive model

We consider a setting with $p$ partially observed covariates
$X=X_{1},...,X_{p}$, $q$ fully observed covariates $Z=Z_{1},...,Z_{q}$, and
$K=2$ competing events. We assume that (possibly conditional on $Z$) censoring
is independent of both $X$ and the competing risks outcomes
$\tilde{T},\tilde{D}$. We furthermore let $X^{\text{obs}}$ and
$X^{\text{mis}}$ respectively denote the observed and missing components of
$X$ for an individual, and let $R$ be the vector of observation indicators
(equal to 1 if the corresponding element of $X$ is observed, or equal to 0 if
it is missing).

The substantive model of interest is
$\lambda_{1}(t\,|\,X,Z)=\lambda_{0}(t)\exp\\{g(X,Z;\boldsymbol{\beta})\\}$,
which is a Fine–Gray model for cause 1, and where $g(X,Z;\boldsymbol{\beta})$
is a function of $X$ and $Z$, parametrised by $\boldsymbol{\beta}$. In this
section, we provide an overview of possible approaches for imputing each
partially observed $X_{j}$. In addition to an approach which imputes
compatibly with the assumed substantive model, we also consider alternative
methods which are either a) only approximately compatible with the substantive
model; or b) impute assuming a different underlying competing risks structure
(i.e. cause-specific proportional hazards). We require that the proposed
approaches be valid under the missing-at-random (MAR) assumption, that is,
$P(R\,|\,T,D,X,Z)=P(R\,|\,T,D,X^{\text{obs}},Z)$.

###  3.1 MI based on cause-specific hazards models

####  3.1.1 CS-SMC

A first MI approach to consider is to impute compatibly with cause-specific
Cox models, despite the substantive model of interest being a Fine–Gray model
for cause 1. As described by Bartlett and Taylor8, this method relies on the
substantive-model-compatible imputation density for $X_{j}$, given by

$f(X_{j}\,|\,T,D,X_{-j},Z)\propto f(T,D\,|\,X,Z)f(X_{j}\,|\,X_{-j},Z),$ (2)

where $X_{-j}$ refers to the components of $X$ after removing $X_{j}$, and
$f(\cdot)$ is a density function. For example, $f(T,D\,|\,X,Z)$ is used as
shorthand notation for $f_{T,D\,|\,X,Z}(t,d\,|\,x,z)$, that is, the density
function for the conditional distribution $T,D\,|\,X,Z$, evaluated at $(t,d)$
for given values $x$ and $z$.

In practice, the substantive model $f(T,D\,|\,X,Z;\psi)$ assumed for
$f(T,D\,|\,X,Z)$ is a cause-specific Cox model (one for each competing risk).
Therefore, $\psi$ ($\psi\in\Psi$) contains the cumulative baseline hazards and
log hazard ratios for each cause-specific hazard. A model
$f(X_{j}\,|\,X_{-j},Z;\phi)$ indexed by $\phi$ ($\phi\in\Phi$), is also
assumed for $f(X_{j}\,|\,X_{-j},Z)$. The idea is then to sample candidate
imputed values for the missing $X_{j}$ using $f(X_{j}\,|\,X_{-j},Z;\phi)$, and
accept these if they also represent draws from a density proportional to
$f(T,D\,|\,X,Z;\psi)f(X_{j}\,|\,X_{-j},Z;\phi)$. We refer to this method as
the cause-specific SMC-FCS approach (CS-SMC).

####  3.1.2 CS-Approx

The approximately compatible analogue to the cause-specific SMC-FCS approach
is described by Bonneville et al. 5. As briefly described in the introduction,
this approach involves directly specifying an imputation model
$f(X_{j}\,|\,T,D,X_{-j},Z;\alpha)$ for $f(X_{j}\,|\,T,D,X_{-j},Z)$. In order
to ensure approximate compatibility with assumed cause-specific Cox
substantive models, the imputation model should include as predictors
$X_{-j}$, $Z$, $D$ (as a factor variable), and the (marginal, as obtained
using the Nelson–Aalen estimator) cause-specific cumulative hazard for each
cause $\hat{H}_{k}(T)$, evaluated at an individual’s event or censoring time.
We refer to this method as approximately compatible cause-specific MICE (CS-
Approx).

####  3.1.3 MI based on the relation between the cause-specific and
subdistribution hazards

The imputations generated by the CS-SMC and CS-Approx approaches will
typically not be consistent with the assumption of proportional
subdistribution hazards for cause 1 made by the substantive model of interest.
This is because, for cause 1, proportionality on the cause-specific hazard
scale will generally imply non-proportionality on the subdistribution hazard
scale 14. One can derive the functional form of these time-varying covariate
effects on the subdistribution hazard scale by using the relation between the
subdistribution hazard and the cause-specific hazards15. The CS-SMC and CS-
Approx approaches can therefore be thought of as procedures to impute
(approximately) compatibly with a Fine–Gray model with time-varying covariate
effects, the functional form of which are determined by the assumptions made
for the cause-specific Cox models of each competing event.

A relevant question at this point is whether the relation between cause-
specific and subdistribution hazards can instead be used as part of a
procedure to impute compatibly with proportional subdistribution hazards for
cause 1. In order to motivate such a procedure, we first note that the
conditional density of the observed outcome given covariates used in Equation
(2) can be written both in terms of cause-specific hazards, and in terms of
the cumulative incidence functions, as

$\displaystyle f(T,D\mid X,Z)$
$\displaystyle=\\{h_{1}(T\,|\,X,Z)S(T\,|\,X,Z)\\}^{I(D=1)}\\{h_{2}(T\,|\,X,Z)S(T\,|\,X,Z)\\}^{I(D=2)}$
$\displaystyle\qquad\times S(T\,|\,X,Z)^{1-I(D=1)-I(D=2)},$
$\displaystyle=f_{1}(T\,|\,X,Z)^{I(D=1)}f_{2}(T\,|\,X,Z)^{I(D=2)}$
$\displaystyle\qquad\times\\{1-F_{1}(T\,|\,X,Z)-F_{2}(T\,|\,X,Z)\\}^{1-I(D=1)-I(D=2)},$
(3)

with
$f_{k}(t\,|\,X,Z)=\mathop{}\\!\mathrm{d}F_{k}(t\,|\,X,Z)/\mathop{}\\!\mathrm{d}t$
known as the “subdensity” for cause $k$ 16. These subdensities, in turn, can
be expressed in terms of the subdistribution hazard, as

$\displaystyle f_{k}(t\,|\,X,Z)$
$\displaystyle=\lambda_{k}(t\,|\,X,Z)\\{1-F_{k}(t\,|\,X,Z)\\},$
$\displaystyle=\lambda_{k}(t\,|\,X,Z)\exp\\{-\Lambda_{k}(t\,|\,X,Z)\\}.$ (4)

Specifying a Fine–Gray model for cause 1 is an assumption regarding only part
of Equation (3), namely for any terms involving $f_{1}(T\,|\,X,Z)$. The
practical implication of this is that Equation (2) cannot be used to impute
the missing $X_{j}$ without making assumptions about cause 2. One could for
example assume (for imputation purposes) a cause-specific Cox model for cause
2, derive the implied $h_{1}(t\,|\,X,Z)$ using the relation between the
subdistribution hazard and the cause-specific hazards, and then use both
cause-specific hazards to evaluate $f(T,D\,|\,X,Z)$ in (2).

Given that a Fine–Gray model is assumed for cause 1, some computational
difficulties can be encountered while making assumptions for cause 2. For
example, specifying a Fine–Gray model also for cause 2 in the imputation
procedure could result in the total failure probability at an observed event
time $F_{1}(T\,|\,X,Z)+F_{2}(T\,|\,X,Z)$ exceeding 1, meaning we would not be
able to draw imputed values using (2) for high-risk individuals17. An
additional example concerns the approach described in the previous paragraph,
where $h_{1}(t\,|\,X,Z)$ is derived based on $h_{2}(t\,|\,X,Z)$ and
$\lambda_{1}(t\,|\,X,Z)$. The numerical integration step generally needed to
compute $h_{1}(t\,|\,X,Z)$ could make the overall imputation procedure rather
computationally inefficient. More details on potential issues when specifying
a model for cause 2 when a Fine–Gray model is assumed for cause 1 can be found
in Bonneville et al. 18.

The above points mean that it is desirable to use an alternative approach
which avoids having to specify a model for the cause-specific (or
subdistribution) hazard of cause 2. In the next subsection, we propose a SMC-
FCS approach assuming a Fine–Gray substantive model for cause 1, which avoids
making explicit modelling assumptions concerning cause 2.

###  3.2 MI based on the Fine–Gray model

####  3.2.1 FG-SMC

Suppose for now that the potential censoring time $C$ is known for all
individuals. This is for example the case when there is a fixed end of study
date (i.e. “administrative” censoring), and no additional random right-
censoring. Fine and Gray referred to these kind of data as “censoring
complete”, since the subdistribution at-risk process is known9. Equivalently,
the “observed” subdistribution random variable for cause 1 (henceforth
referred to as “subdistribution time”), $V=I(D=1)\times T+I(D\neq 1)\times C$,
is known for all individuals. In turn, this implies that (with complete
covariate data), the Fine–Gray model can be estimated by fitting a standard
Cox model with outcome $V$ and event indicator $I(D=1)$.

Consequently, an intuitive approach to imputing the missing $X_{j}$ in our
setting might therefore be to apply existing SMC-FCS methodology for standard
Cox models (see section 6.3 of Bartlett et al.7), but instead using $V$ and
$I(D=1)$ as our outcome variables. We refer to this method as Fine–Gray SMC-
FCS (FG-SMC). The substantive-model-compatible imputation density is now

$f(X_{j}\,|\,V,D,Z)\propto f(V,D\,|\,X,Z)f(X_{j}\,|\,X_{-j},Z),$ (5)

where the conditional density of the observed outcome given the covariates can
be written as

$\displaystyle f(V,D\mid X,Z)$ $\displaystyle=f_{1}(V\mid
X,Z)^{I(D=1)}\\{1-F_{1}(V\mid X,Z)\\}^{I(D\neq 1)},$
$\displaystyle=[\lambda_{1}(V\,|\,X,Z)\exp\\{-\Lambda_{1}(V\,|\,X,Z)\\}]^{I(D=1)}\exp\\{-\Lambda_{1}(V\,|\,X,Z)\\}^{I(D=0)}$
$\displaystyle\qquad\exp\\{-\Lambda_{1}(V\,|\,X,Z)\\}^{I(D=2)},$
$\displaystyle=\lambda_{1}(V\,|\,X,Z)^{I(D=1)}\exp\\{-\Lambda_{1}(V\,|\,X,Z)\\},$
(6)

using Equation (4) and the fact that
$f_{1}(V\,|\,X,Z)^{I(D=1)}=f_{1}(T\,|\,X,Z)^{I(D=1)}$. Note that while
Equations (5) and (6) depend only on $I(D=1)$, we still use $D$ in the
notation to make the contribution of those failing from cause 2 to the density
explicit, which is relevant for the upcoming sections. Importantly, this
procedure relies on a stronger MAR assumption (compared to the one introduced
at the beginning of Section 3), namely
$P\\{R\,|\,V,I(D=1),X,Z\\}=P\\{R\,|\,V,I(D=1),X^{\text{obs}},Z\\}$. In
essence, we ignore any terms involving $f_{2}(T\,|\,X,Z)$ in Equation (2)
based on the assumption that missingness in $X$ does not depend on either
$I(D=2)$ or the failure time for those failing from cause 2.

####  3.2.2 FG-Approx

The form of Equation (6) mirrors the likelihood in the standard Cox context,
which can be obtained by replacing $\lambda_{1}(V\,|\,X,Z)$ with the hazard of
a single event (in absence of competing risks). The practical implications of
this for our MI context are that the findings of White and Royston6 in the
single-event survival setting should in principle extend to the Fine–Gray
context. Namely, that the (approximately compatible) directly specified
imputation model $f(X_{j}\,|\,V,D,X_{-j},Z;\alpha)$ for a partially observed
$X_{j}$ should contain as predictors at least $X_{-j}$, $Z$, the indicator for
the competing event of interest $I(D=1)$, and the cumulative subdistribution
baseline hazard for the same event $\Lambda_{0}(V)$. Instead of the unknown
true $\Lambda_{0}(V)$, one could use the estimated marginal cumulative
subdistribution hazard $\hat{\Lambda}_{1}(V)$ instead, obtained using the
Nelson–Aalen estimator using $V$ and $I(D=1)$ are outcome variables. We refer
to this approximately compatible MICE approach as FG-Approx.

###  3.3 Accommodating random right-censoring

In addition to (deterministic) administrative censoring, random right-
censoring may occur. In the presence of random right-censoring, the
contribution of those failing from cause 2 to density (6) is no longer
evaluable, since we do not know their potential censoring time. In a sense,
their subdistribution time has been informatively censored by their cause 2
failure.

####  3.3.1 Via imputation of potential censoring times

One approach to estimate the parameters of a Fine–Gray model in the presence
of random right censoring is to consider the potential censoring times for
those failing from cause 2 as missing data, and multiply impute them. To this
end, Ruan and Gray13 suggested the use of Kaplan–Meier (KM) imputation 19.
Specifically, potential censoring times are randomly drawn from the
conditional distribution with distribution function
$1-P(C>t\,|\,C>T)=1-\hat{G}(t-)/\hat{G}(T-)$, where $\hat{G}(t)$ is a KM
estimate of the survival distribution of the censoring times $P(C>t)$. The
imputation of these potential censoring times effectively produces multiple
“censoring complete” datasets, in which a Fine–Gray model can be fit using
standard software. Inference is then based on a pooled model, which combines
the models fitted in each censoring complete dataset using Rubin’s rules 20.

We can make use of the above ideas in order to multiply impute covariates
compatibly with a Fine–Gray model in the presence of random right random
censoring. Specifically, we can apply the FG-SMC (or FG-Approx) method in each
censoring complete dataset obtained after first imputing the potential
censoring times for those failing from cause 2. In order to formalise this
procedure, recall that $\boldsymbol{\beta}$ represents the parameters of the
substantive model, and that $X=\\{X^{\text{obs}},X^{\text{mis}}\\}$. We can
similarly partition $V=\\{V^{\text{obs}},V^{\text{mis}}\\}$, where
$V^{\text{mis}}$ is the vector of missing censoring times for those failing
from cause 2.

From a Bayesian perspective, the goal is to estimate the conditional density
of $\boldsymbol{\beta}$ given the observed data, namely

$\displaystyle f(\boldsymbol{\beta}\,|\,X^{\text{obs}},Z,V^{\text{obs}},D)$
$\displaystyle=\int_{V}\int_{X}f(\boldsymbol{\beta}\,|\,X^{\text{obs}},X^{\text{mis}},Z,V^{\text{obs}},V^{\text{mis}},D)\times$
$\displaystyle\qquad
f(X^{\text{mis}},V^{\text{mis}}\,|\,X^{\text{obs}},Z,V^{\text{obs}},D)\mathop{}\\!\mathrm{d}X^{\text{mis}}\mathop{}\\!\mathrm{d}V^{\text{mis}}.$
(7)

If we can sample imputed values $M$ times from
$f(X^{\text{mis}},V^{\text{mis}}\,|\,X^{\text{obs}},Z,V^{\text{obs}},D)$, the
integral above can be approximated by an average over
$f(\boldsymbol{\beta}\,|\,X^{\text{obs}},X^{\text{mis}},Z,V^{\text{obs}},V^{\text{mis}},D)$
(the “complete data” posterior density) evaluated at those $M$ moments 21.

One option to sample from
$f(X^{\text{mis}},V^{\text{mis}}\,|\,X^{\text{obs}},Z,V^{\text{obs}},D)$, the
joint posterior predictive density, is to iteratively sample from
$f(V^{\text{mis}}\,|\,X,Z,V^{\text{obs}},D)$ and
$f(X^{\text{mis}}\,|\,X^{\text{obs}},Z,V,D)$. This is clearly necessary when
the censoring distribution depends on $X$, since the potential censoring times
will differ depending on the most recently imputed $X$. These imputed
censoring times would have to be based for example on a Cox model for the
censoring hazard instead of a marginal KM estimate.

If however
$C\raisebox{0.40002pt}{\rotatebox[origin={c}]{90.0}{$\models$}}X\,|\,Z$, we
could use a sequential approach, where we factorise

$\displaystyle
f(X^{\text{mis}},V^{\text{mis}}\,|\,X^{\text{obs}},Z,V^{\text{obs}},D)$
$\displaystyle=f(X^{\text{mis}}\,|\,X^{\text{obs}},Z,V,D)f(V^{\text{mis}}\,|\,X^{\text{obs}},Z,V^{\text{obs}},D),$
$\displaystyle=f(X^{\text{mis}}\,|\,X^{\text{obs}},Z,V,D)f(V^{\text{mis}}\,|\,Z,V^{\text{obs}},D).$

Practically speaking, this involves imputing the potential censoring times
(possibly in strata of $Z$) in a first step, and then imputing the missing $X$
in a second step. This can be implemented easily using existing software
packages in R: {kmi} for the imputation of censoring times 22, and {smcfcs}
for the imputation of the missing covariates 23 \- see Figure 1 for an
illustration of the workflow.

Note that the described KM-based procedure for imputing potential censoring
times does not take into account any of the uncertainty in estimating
$P(C>t)$. Ruan and Gray discussed using the non-parametric bootstrap to
account for this uncertainty and improve estimation properties, and found
similar results both with and without a bootstrap step13. In Appendix A, we
visualise and give additional details concerning the imputation of potential
censoring times.

Figure 1: Sequential workflow for (compatible) covariate imputation and
analysis for a Fine–Gray substantive model with two covariates $X$ and $Z$, in
the presence of random right-censoring. In the first step, the potential
censoring times for those failing from cause 2 are multiple imputed using the
{kmi} package. In the second step, the missing covariates are imputed using
the {smcfcs} (or {mice}) package. This workflow is valid when the probability
of being censored is independent of $X$ and any $Z$ related to the censoring
process are modelled in {kmi}.

####  3.3.2 Via censoring weights in the likelihood

Rather than multiply imputing the potential censoring times, an alternative
approach is to incorporate inverse probability of censoring weights directly
in (6). If we define time-dependent weights

$\displaystyle w(t)$ $\displaystyle=1$ $\displaystyle\text{if }t\leq T,$
$\displaystyle w(t)$ $\displaystyle=P(C>t\,|\,C>T)=\frac{G(t-)}{G(T-)}$
$\displaystyle\text{if }t>T,$

then the conditional density of the (subdistribution) outcome given the
covariates can be written as

$\displaystyle f(V,D\mid X,Z)$
$\displaystyle=[\lambda_{1}(V\,|\,X,Z)\exp\\{-\Lambda_{1}(V\,|\,X,Z)\\}]^{I(D=1)}\exp\\{-\Lambda_{1}(V\,|\,X,Z)\\}^{I(D=0)}\times$
$\displaystyle\qquad\exp\Biggl{\\{}-\int_{0}^{\infty}w(u)\lambda_{1}(u\,|\,X,Z)\mathop{}\\!\mathrm{d}u\Biggr{\\}}^{I(D=2)},$
(8)

where the term for those failing from the competing event involves integration
in practice up to a maximum potential follow-up time $t^{*}$. As described by
Lambert et al., this integral can be approximated by splitting time into
intervals, in which the corresponding $w(t)$ is assumed to be constant24.

The integration step needed for those failing from cause 2 in (8) means that
this approach cannot be implemented in a straightforward way with existing
software, unlike the approach described in the previous subsection. The
simulation study in this paper therefore focuses on the approach involving
multiple imputation of potential censoring times.

###  3.4 Implementation of MI approaches

Methods CS-SMC, CS-Approx, FG-SMC, and FG-Approx can all be implemented using
existing software packages in R. In this section, we summarise the steps
needed to apply these methods in a given dataset in the presence of random
right censoring (possibly in combination with administrative censoring). A
minimal R code example can be found in supplementary material S1.

  1. 1.

Add columns $\hat{H}_{1}(T)$ and $\hat{H}_{2}(T)$ to the original data, which
are the marginal cause-specific cumulative hazards for each competing risk
evaluated at an individual’s event or censoring time (obtained using the
Nelson–Aalen estimator).

  2. 2.

Multiply impute the potential censoring for those failing from cause 2 using
{kmi}, yielding $m$ censoring complete datasets (i.e. with “complete” $V$).
The censoring distribution has support at both random and administrative
censoring times. Any completely observed covariates that are known to affect
the probability of being censored should be included as predictors in the
model for the censoring process. {kmi} imputes based on stratified KM when $Z$
are categorical, and based on a Cox model at least one of $Z$ is continuous.
If for example an individual’s time of entry into a study determines their
maximum follow-up duration, this should be accounted for in the imputation
procedure (e.g. by stratifying by year of entry).

  3. 3.

In each censoring complete dataset, add an additional column
$\hat{\Lambda}_{1}(V)$. This takes the value of the marginal cumulative
subdistribution hazard for cause 1 at an individual’s observed or imputed
subdistribution time, obtained with the Nelson–Aalen estimator based on
$I(D=1)$ and imputed $V$.

  4. 4.

In each censoring complete dataset (each with different $V$ and
$\hat{\Lambda}_{1}(V)$, but same $\hat{H}_{1}(T)$ and $\hat{H}_{2}(T)$),
create a single imputed dataset using the desired covariate imputation
method(s):

     * •

CS-SMC: use {smcfcs} to impute the missing covariate(s) compatibly with cause-
specific Cox models. All covariates used in the Fine–Gray substantive model
should feature in at least one of the specified cause-specific models.

     * •

CS-Approx: use {mice} to impute the missing covariate(s), where the imputation
model contains as predictors the remaining substantive model covariates, $D$
(as a factor variable), and both $\hat{H}_{1}(T)$ and $\hat{H}_{2}(T)$.

     * •

FG-SMC: use {smcfcs} to impute the missing covariate(s) compatibly with the
Fine–Gray substantive model. This is done by using the imputation methods
developed for the standard Cox model, but with as outcome variables $I(D=1)$
and imputed $V$.

     * •

FG-Approx: use {mice} to impute the missing covariate(s), where the imputation
model contains as predictors the remaining substantive model covariates,
$I(D=1)$, and $\hat{\Lambda}_{1}(V)$.

  5. 5.

Fit the Fine–Gray substantive model in each imputed dataset (using standard
Cox software with $I(D=1)$ and imputed $V$ as outcome variables), and pool the
estimates using Rubin’s rules.

##  4 Simulation study

We aim to evaluate the performance of different MI methods in the presence of
missing covariate data when specifying a Fine–Gray model for the
subdistribution hazard for one event of interest in the presence of one
competing event. Specifically, we assume interest lies in the estimation (for
cause 1) of both subdistribution hazard ratios, and the cumulative incidence
for a particular individual at some future time horizon. We follow the ADEMP
structure for the reporting of the simulation study25.

###  4.1 Data-generating mechanisms

We generate datasets of $n=2000$ individuals, with two covariates $X$ and $Z$.
We assume $Z\sim\mathcal{N}(0,1)$ and
$X\,|\,Z\sim\text{Bernoulli}\\{(1+e^{-Z})^{-1}\\}$.

We let $h_{k}(t\,|\,X,Z)$, $\lambda_{k}(t\,|\,X,Z)$ and
$F_{k}(t\,|\,X,Z)=P(\tilde{T}\leq t,\tilde{D}=k\,|\,X,Z)$ respectively denote
the cause-specific hazards, subdistribution hazards and cumulative incidence
functions for cause $k$, conditional on $X$ and $Z$. The competing event times
will be generated following two mechanisms: one where the Fine–Gray model for
cause 1 is correctly specified, and another where it is misspecified. These
are detailed below, together with assumptions concerning both censoring and
the missing data mechanisms.

####  4.1.1 Correctly specified Fine–Gray

For this mechanism, we simulate data using the ‘indirect’ method described in
Beyersmann et al.14, and originally used in the simulations by Fine and Gray
9. This approach involves first drawing the competing event indicator
$\tilde{D}$, and then generating an event time for those with $\tilde{D}=1$.
The final step is to generate times of the competing event for the remaining
individuals, who were assigned $\tilde{D}=2$.

Here, we directly specify the cumulative incidence of cause 1 as

$F_{1}(t\,|\,X,Z)=1-\big{[}1-p\\{1-\exp(-b_{1}t^{a_{1}})\\}\big{]}^{\exp(\beta_{1}X+\beta_{2}Z)}.$

The above expression corresponds to a Fine–Gray model, with as baseline
cumulative incidence function a Weibull cumulative distribution function with
shape $a_{1}$ and rate $b_{1}$ (parametrisation used in Klein and
Moeschberger26) multiplied by a probability $p$. Explicitly,

$F_{0}(t)=p\\{1-\exp(-b_{1}t^{a_{1}})\\}.$

With $\lim_{t\to\infty}F_{0}(t)=p$, we have that
$P(\tilde{D}=1\,|\,X,Z)=1-(1-p)^{\exp(\beta_{1}X+\beta_{2}Z)}$, and
$P(\tilde{D}=2\,|\,X,Z)=1-P(\tilde{D}=1\,|\,X,Z)=(1-p)^{\exp(\beta_{1}X+\beta_{2}Z)}$.
These are the individual-specific cumulative incidences for each event at time
infinity. Also note that the baseline subdistribution hazard for this
mechanism can be obtained by
$\\{\mathop{}\\!\mathrm{d}F_{0}(t)/\mathop{}\\!\mathrm{d}t\\}\times\\{1-F_{0}(t)\\}^{-1}$.

The idea then is to generate the event times for cause 1 conditionally on the
event indicator and covariates, using

$\displaystyle P(\tilde{T}\leq t\,|\,\tilde{D}=1,X,Z)$
$\displaystyle=\frac{P(\tilde{T}\leq
t,\tilde{D}=1\,|\,X,Z)}{P(\tilde{D}=1\,|\,X,Z)}$
$\displaystyle=\frac{1-\big{[}1-p\\{1-\exp(-b_{1}t^{a_{1}})\\}\big{]}^{\exp(\beta_{1}X+\beta_{2}Z)}}{1-(1-p)^{\exp(\beta_{1}X+\beta_{2}Z)}}.$
(9)

To sample from the above, we first need to draw
$\tilde{D}\sim\text{Bernoulli}\\{(1-p)^{\exp(\beta_{1}X+\beta_{2}Z)}\\}+1$. We
can then use inverse transform sampling to draw failure times within the
subset of individuals with $\tilde{D}=1$. Shortening
$\exp(\beta_{1}X+\beta_{2}Z)=\exp(\eta)$, and with $u\sim\mathcal{U}(0,1)$, we
can invert (9) as

$t=\Biggl{[}-\frac{1}{b_{1}}\log\Bigg{[}1-\frac{1-\big{[}1-u\\{1-(1-p)^{\exp(\eta)}\\}\big{]}^{1/\exp(\eta)}}{p}\Bigg{]}\Biggr{]}^{1/{a_{1}}}.$

For the competing event, we can factorise the cumulative incidence function as

$P(\tilde{T}\leq t,D=2\,|\,X,Z)=P(\tilde{T}\leq
t\,|\,\tilde{D}=2,X,Z)P(\tilde{D}=2\,|\,X,Z).$

A proportional hazards model can then be specified (for convenience) for

$P(\tilde{T}\leq
t\,|\,\tilde{D}=2,X,Z)=1-\exp\Bigl{\\{}-H_{20}^{*}(t)\exp(\beta_{1}^{*}X+\beta_{2}^{*}Z)\Bigr{\\}},$

where $H_{20}^{*}(t)$ is the cumulative baseline hazard associated to the
cumulative incidence function conditional on $\tilde{D}=2$. Since the event
indicator is already drawn, the failure times can be drawn using standard
methods within the subset with $\tilde{D}=2$. Here, we specify a Weibull
baseline hazard as $h_{20}^{*}(t)=a_{2}b_{2}t^{a_{2}-1}$.

We fix
$\\{\beta_{1},\beta_{2},\beta_{1}^{*},\beta_{2}^{*}\\}=\\{0.75,0.5,0.75,0.5\\}$,
and the Weibull parameters used for both events as shape
$\\{a_{1},a_{2}\\}=0.75$ and rate $\\{b_{1},b_{2}\\}=1$. We vary
$p=\\{0.15,0.65\\}$, which is the expected proportion of event 1 failures for
individuals with $X=0$ and $Z=0$.

####  4.1.2 Simulation based on cause-specific hazards (misspecified
Fine–Gray)

In this data-generating mechanism (DGM), we assume proportionality on the
cause-specific hazard scale, and simulate using latent failure times 27. We
specify baseline Weibull hazards for both cause-specific hazards as

$\displaystyle h_{1}(t\,|\,X,Z)$
$\displaystyle=a_{1}b_{1}t^{a_{1}-1}\exp(\gamma_{11}X+\gamma_{12}Z),$
$\displaystyle h_{2}(t\,|\,X,Z)$
$\displaystyle=a_{2}b_{2}t^{a_{2}-1}\exp(\gamma_{21}X+\gamma_{22}Z),$

where $\\{a_{1},a_{2}\\}$ and $\\{b_{1},b_{2}\\}$ are respectively the shape
and rate parameters. Under this DGM, a Fine–Gray model for cause 1 will be
misspecified. Nevertheless, the coefficients resulting from the misspecified
Fine–Gray model could still be interpreted as time-averaged effects on the
(complementary log-log transformed) cumulative incidence function28.

We aim to have a scenario close to the one described in 4.1.1 (in terms of
event proportions), where the main difference is that proportionality now
holds on the cause-specific hazard scale. To fix the parameters in this DGM,
we first simulate a large dataset of one million individuals following the
mechanism described in the previous subsection, where proportional
subdistribution hazards hold. Parametric cause-specific proportional hazards
models assuming baseline Weibull hazards are then fitted for each failure
cause. The point estimates obtained from these models are used as the cause-
specific data-generating parameters
$\\{a_{1},b_{1},\gamma_{11},\gamma_{12}\\}$ and
$\\{a_{2},b_{2},\gamma_{21},\gamma_{22}\\}$. These parameters will of course
differ depending on $p=\\{0.15,0.65\\}$, and also depending on the censoring
distribution. While the cause-specific models fitted on this large dataset
will be misspecified (cause-specific baseline hazards are not of Weibull
shape, and covariates effects on the cause-specific hazards are non-
proportional), the resulting “least false” parameters are still useful.

Figure 2 summarises the DGMs, prior to the addition of any censoring. In the
correctly specified Fine–Gray scenarios, the subdistribution log hazard ratio
$\lambda_{1}(t\,|\,X=1,Z)/\lambda_{1}(t\,|\,X=0,Z)$ is time constant, while
the cause-specific log hazard ratios are time-dependent. The reverse is true
for the misspecified Fine–Gray scenarios. Overall, the correctly specified and
misspecified Fine–Gray scenarios are very comparable in terms of (true)
baseline hazards and cumulative incidences, for both values of $p$.

Figure 2: Summary of data-generating mechanisms prior to the addition of any
censoring. For each value of $p$, both the correctly specified and
misspecified Fine–Gray (FG) scenarios are very comparable in terms of (true)
baseline hazards and cumulative incidences.

####  4.1.3 Censoring

The DGMs outlined above assume no loss to follow-up. As additional scenarios,
we consider independent (i.e. not conditional on any covariates) right
censoring where the censoring times are simulated from an exponential
distribution with rate $\lambda_{C}=0.49$, resulting in approximately 30% of
censored observations. These censoring times will be considered as either: a)
known (administrative censoring); or b) unknown (random censoring).

####  4.1.4 Covariate missingness

Missingness is induced in $X$, while $Z$ remains fully observed. Let $R_{X}$
be a binary variable indicating whether $X$ is missing ($R_{X}=0$) or observed
($R_{X}=1$). We use a missing at random (MAR) mechanism conditional on $Z$,
which was defined as
$\operatorname{logit}P(R_{X}=0\,|\,Z)=\eta_{0}+\eta_{1}Z$. We take
$\eta_{1}=1.5$, a rather strong mechanism where higher values of $Z$ are
associated with more missingness in $X$. The value of $\eta_{0}$ is found via
standard root solving, such that the average probability
$P(R_{X}=0)=\mathbb{E}\\{P(R_{X}=0\,|\,Z)\\}$ of being missing in a given
dataset equals 0.4.

####  4.1.5 Summary

In summary, the simulation study varied

  * •

Censoring type: no censoring, administrative and random censoring

  * •

Relative occurrence of event 1, as low or high. This is done by varying the
baseline cumulative incidence of event 1 (as $t\rightarrow\infty$) as
$p=\\{0.15,0.65\\}$.

  * •

Failure time simulation methods, with a) directly specified cumulative
incidence cause 1 (correctly specified Fine–Gray); b) cause-specific
proportional hazards for both causes (misspecified Fine–Gray).

This adds up to 3 (censoring types) $\times$ 2 (relative occurence event 1)
$\times$ 2 (failure time simulation methods) = 12 scenarios.

###  4.2 Estimands

The first estimands of interest are the subdistribution log hazard ratios
$\beta_{1}$ and $\beta_{2}$ for $X$ and $Z$, respectively. In the correctly
specified Fine–Gray scenarios, these simply correspond to the data-generating
parameters $\\{\beta_{1},\beta_{2}\\}=\\{0.75,0.5\\}$. In the misspecified
Fine–Gray scenarios however, the target values (the “least-false parameters”;
time averaged subdistribution log hazard ratios
$\\{\tilde{\beta}_{1},\tilde{\beta}_{2}\\}$) are obtained by fitting a
Fine–Gray model on a large simulated dataset of one million individuals,
simulated as under the second data-generating mechanism, after applying any
censoring. For computational efficiency, the censoring times are assumed to be
known when fitting the Fine–Gray model on this large dataset.

The second estimands of interest are the conditional cumulative incidence of
event 1 at a grid of timepoints (between timepoints 0 and 5) for reference
individuals $\\{X,Z\\}=\\{0,0\\}$ (baseline) and $\\{X,Z\\}=\\{1,1\\}$. In the
correctly specified Fine–Gray scenarios, this corresponds to

$F_{1}(t\,|\,X,Z)=1-\big{[}1-p\\{1-\exp(-b_{1}t^{a_{1}})\\}\big{]}^{\exp(\beta_{1}X+\beta_{2}Z)},$

while for the misspecified Fine–Gray scenarios, this corresponds to

$\displaystyle F_{1}(t\,|\,X,Z)$
$\displaystyle=\int_{0}^{t}h_{1}(u\,|\,X,Z)\exp\Bigl{\\{}{-H_{1}(u\,|\,X,Z)-H_{2}(u\,|\,X,Z)\Bigr{\\}}}\mathop{}\\!\mathrm{d}u,$
$\displaystyle=\int_{0}^{t}a_{1}b_{1}u^{a_{1}-1}\exp(\gamma_{11}X+\gamma_{12}Z)$
$\displaystyle\qquad\times\exp\Bigl{\\{}{-b_{1}u^{a_{1}}\exp(\gamma_{11}X+\gamma_{12}Z)-b_{2}u^{a_{2}}\exp(\gamma_{21}X+\gamma_{22}Z)\Bigr{\\}}}\mathop{}\\!\mathrm{d}u,$

which is obtained via numerical integration.

###  4.3 Methods

The assessed methods are

  * •

Full: analysis run on full data prior to missing values, as a benchmark for
the best possible performance.

  * •

CCA: complete-case analysis, as a “lower” benchmark that the imputation
methods need to outperform in order to be worthwhile.

  * •

CS-SMC: MI, imputing compatibly with cause-specific Cox proportional hazards
models. This method is described in Section 3.1.1. Both $X$ and $Z$ are used
as predictors in each cause-specific model assumed by this procedure.

  * •

CS-Approx: MI with both marginal cumulative cause-specific hazards (evaluated
at the individual observed event or censoring time) and competing event
indicator included as predictors in the imputation model, in addition to $Z$.
This method is described in Section 3.1.2.

  * •

FG-SMC: MI, imputing compatibly with a Fine–Gray model for cause 1 that has as
covariates $X$ and $Z$. This is the method described in Section 3.2.1.

  * •

FG-Approx: MI with marginal cumulative subdistribution hazard (evaluated at
the individual observed or imputed subdistribution time $V$) and indicator for
event 1 included as predictors in the imputation model, in addition to $Z$.
This method is described in Section 3.2.2.

The imputation methods are run with 30 imputed datasets. This was fixed
following a pilot set of simulations with 50 imputed datasets, which showed
that there was little reduction in empirical standard errors for the
subdistribution log hazard ratios (and their Monte Carlo standard errors)
beyond 30 imputed datasets. Approximately compatible MI methods CS-Approx and
FG-Approx only require a single iteration because there is just one variable
with missing values, while substantive-model-compatible (SMC) MI methods CS-
SMC and FG-SMC are run with 20 iterations. The method used to model
$f(X\,|\,V,D,Z;\alpha)$ for approximately compatible methods is logistic
regression, while for SMC methods $f(X\,|\,Z;\psi)$ is specified as a logistic
regression. We note that $X$ was chosen to be binary as SMC methods do not
require rejection sampling for variables with discrete sample space, thereby
reducing simulation time.

For the scenarios with no or administrative censoring, the subdistribution
time $V$ is fully observed. While $V=T$ for those failing from cause 1, for
those failing from cause 2, $V$ is first set to either a) a large value
greater than the largest observed event 1 time (in absence of censoring); or
b) the known potential censoring time $C$ (administrative censoring). The
marginal cumulative subdistribution hazard used for the approximate
subdistribution MI method is obtained using a marginal model with $I(D=1)$ and
the resulting $V$ as outcome variables. The covariate MI methods are run once
these $V$ and $I(D=1)$ variables have been created. In scenarios with random
censoring, the potential censoring times for those failing from cause 2 are
multiply imputed using the {kmi} R package with default settings: marginal
non-parametric model for the censoring distribution, and no additional
bootstrap layer. This yields 30 imputed datasets, each with a different $V$.
In each of these datasets, the marginal cumulative subdistribution hazard is
estimated in the same way as described above. Thereafter, the covariate MI
methods are run in each of these datasets, yielding one imputed dataset for
each imputed $V$ (total of 30 imputed datasets), corresponding to the workflow
in Figure 1.

For all methods, the Fine–Gray model for cause 1 is estimated using a Cox
model with (known or imputed) $V$ and $I(D=1)$ as outcome variables. When the
imputation methods are used (and for all methods when there is random right
censoring), the estimated $\hat{\beta}_{1}$ and $\hat{\beta}_{2}$ are the
results of coefficients pooled using Rubin’s rules. Confidence intervals
around these estimates are built as described in Section 2.4.2 in the text by
van Buuren 29. For the cumulative incidences, the estimates for the two sets
of reference values of $X$ and $Z$ are first made in each imputed dataset
using Equation (1), and thereafter pooled using Rubin’s rules after
complementary log-log transformation - as described in Morisot et al.30 and
recommended by Marshall et al.31. This predict-then-pool approach (rather than
predicting using a pooled model) has been recommended by multiple authors 32,
33.

####  4.3.1 Performance measures

The primary measure of interest was bias in the estimated subdistribution log
hazard ratios. In order to keep the Monte Carlo standard error (MCSE) of bias
under a desired threshold of 0.01, we require
$n_{\text{sim}}=0.2^{2}/0.01^{2}=400$ replications per scenario, as we expect
empirical standard errors to be under 0.2 for all scenarios (based on a pilot
run). This number was rounded up to $n_{\text{sim}}=500$. In addition to bias,
we recorded empirical and estimated standard errors, and coverage
probabilities. For the cumulative incidence estimates, we focused on both bias
and root mean square error (RMSE).

####  4.3.2 Software

Analyses were performed using R version 4.3.1 34. Core packages used were:
{survival} version 3.5.7 35, {mice} version 3.16.0 3, {smcfcs} version 1.7.1
23, {kmi} version 0.5.5 22, and {rsimsum} version 0.11.3 36.

###  4.4 Results

We summarise the main findings in this section, with full results available in
a markdown file on the Github repository linked at the end of the present
manuscript.

####  4.4.1 Subdistribution log hazard ratios

We focus on the results for $\beta_{1}$, together with its time-averaged
analogue $\tilde{\beta}_{1}$ in the scenarios with time-dependent
subdistribution hazard ratios. Results concerning bias are summarised in
Figure 3 for all 12 scenarios, and presented on the relative scale (Monte
Carlo standard errors were below the desired 0.01 for both bias and relative
bias, across all methods and scenarios).

When the Fine–Gray model for cause 1 was correctly specified, the proposed FG-
SMC approach was unbiased regardless of censoring type or (baseline)
proportion of cause 1 failures. In contrast, imputing compatibly with the
(incorrect) assumption of proportional cause-specific hazards showed strong
biases, particularly when $p=0.15$ in the absence of censoring (25% biased).
In the presence of censoring however, this bias dropped to approximately 10%.
The CS-Approx method showed consistent downward biases regardless of $p$ and
censoring type, while the FG-Approx method was only biased when $p=0.65$. The
latter finding is consistent with previous research in the simple survival
setting; namely that the approximately compatible MI approach is expected to
work well when cumulative incidence is low6. When the DGM generated event
times under proportional cause-specific hazards, the magnitude of any biases
present were in general smaller (e.g. closer to the 5% mark for approximate MI
approaches when $p=0.65$). For the FG-SMC approach, bias was most noticeable
when $p=0.65$, and in the absence of censoring. CS-SMC was unbiased throughout
these misspecified Fine–Gray scenarios.

Figure 4 summarises empirical and model-based standard errors, together with
coverage probabilities for $\beta_{1}$ and $\tilde{\beta}_{1}$. The model-
based standard errors were on average close to their empirical counterparts.
CS-SMC appears to have a slight variance advantage over competing approaches,
mainly when $p=0.15$. Interestingly, there was no gain in efficiency when the
censoring times were known compared to when they needed to be imputed. This is
in line with simulation results in both Fine and Gray9 and Ruan and Gray13,
that compared the censoring complete variance estimator (of subdistribution
log hazard ratios) to estimators based on the weighted score function and KM
imputation method, respectively. The FG-SMC approach showed good coverage
(near the nominal 95% mark) when the Fine–Gray model was correctly specified,
although there was slight over-coverage when imputation of censoring times was
required. Using the non-parametric bootstrap when estimating $P(C>t)$, which
was not investigated in the simulation study, is unlikely to correct for this
over-coverage. Under-coverage showed by competing approaches was primarily due
to biased estimates.

####  4.4.2 Individual-specific cumulative incidences

Figure 5 shows the true and average estimated baseline cumulative incidence
function $F_{0}(t)$, the average difference between true and estimated
$F_{0}(t)$, and the RMSE of the estimates. Figure 6 presents the same
information instead for a patient with $\\{X,Z\\}=\\{1,1\\}$. Scenarios where
the censoring times are known are omitted from the Figure, as results were
indistinguishable from scenarios where the censoring times needed to be
imputed.

The cost of imputing compatibly with the wrong model (using CS-SMC when the
Fine–Gray model was correctly specified, or FG-SMC when the DGM was based on
cause-specific proportional hazards) when estimating $F_{0}(t)$ was only
noticeable for the CS-SMC approach in the absence of censoring when $p=0.15$,
in terms of both absolute bias and RMSE. On the whole, the approximately
compatible MI approaches performed comparably in terms of RMSE to the SMC
approaches. In scenarios where the Fine–Gray model was misspecified, the
effect of substantive model misspecification (post-imputation) was clear to
see in terms of estimating $F_{0}(t)$ (over- and underestimation at different
points in time). When $\\{X,Z\\}=\\{1,1\\}$, all imputation approaches
outperformed CCA in terms of RMSE when estimating $F_{1}(t\,|\,X=1,Z=1)$,
though to a lesser extent when $p=0.65$. This can presumably be attributed to
the efficiency gain in estimating $\beta_{2}$.

Figure 3: Relative bias (%) in estimating $\beta_{1}$. For the correctly
specified Fine–Gray (FG) scenarios, $\beta_{1}=0.75$. In the misspecified FG
scenarios, the value of the “least-false” $\tilde{\beta}_{1}$ (time-averaged
log subdistribution hazard ratio) depended on both $p$ and the
presence/absence of censoring. For $p=0.15$, $\tilde{\beta}_{1}\approx 0.76$
without censoring, and $\tilde{\beta}_{1}\approx 0.93$ with censoring. For
$p=0.65$, $\tilde{\beta}_{1}\approx 0.75$ both with and without censoring. The
rectangles drawn for each method and scenario represent the estimated relative
bias (vertical line in the middle of a rectangle), and its 95% Monte Carlo
confidence interval (lower and upper limits are represented by the left and
right rectangle borders, respectively). The confidence interval is constructed
using the standard normal approximation. Figure 4: Summary of empirical and
model-based standard errrors, together with coverage probabilities for
$\beta_{1}$ (or $\tilde{\beta}_{1}$ in scenarios with non-proportional
subdistribution hazards). Monte Carlo standard errors (MCSEs) are numerically
given in brackets for the coverage probabilities, while the MCSEs for model-
based/empirical standard errors are shown by 95% confidence intervals (using
standard normal approximation) around a given point. The MCSEs for model-based
standard errors are smaller than the graphical width of the point itself, and
thus cannot be seen. Figure 5: Per scenario (row) for a baseline individual
$\\{X,Z\\}=\\{0,0\\}$: true (black line) versus estimated cumulative incidence
over time, averaged across the 500 replications per scenario (left column);
difference between estimated and true (middle column); root mean square error
(RMSE) of these estimates (right column). Results for scenarios with
administrative censoring are omitted since they were indistinguishable from
those with random censoring. Figure 6: Per scenario (row) for individual
$\\{X,Z\\}=\\{1,1\\}$: estimated versus true (black line) cumulative incidence
over time, averaged across the 500 replications per scenario (left column);
difference between estimated and true (middle column); root mean square error
(RMSE) of these estimates (right column). Results for scenarios with
administrative censoring are omitted since they were indistinguishable from
those with random censoring.

##  5 Applied data example

We illustrate the methods assessed in the simulations study on a dataset of
3982 adult patients with primary and secondary myelofibrosis undergoing
alloHCT between 2009 and 2019, and registered with the EBMT 37. Myelofibrosis
is a rare and chronic myeloproliferative neoplasm characterised by bone marrow
fibrosis and extramedullary hematopoiesis, for which an alloHCT is the only
curative treatment 38. In the original study, the primary objective was to
evaluate the association between comorbidities at time of alloHCT and (cause-
specific) death without prior relapse of the underlying disease, the so-called
non-relapse mortality. In the present illustration, we instead assume that
interest lies in developing a prognostic model for time to disease relapse in
the first 60 months following an alloHCT. To this end, we developed a
Fine–Gray model for relapse, with death prior to relapse as sole competing
risk.

A set of 18 baseline predictors were chosen on the basis of substantive
clinical knowledge, many of which had a considerable proportion of missing
data (see supplementary material S2.1). These predictors included the 13
variables used in the multivariable models from the original study, and 5
additional variables that were either known to be predictive of disease
relapse (use of T-cell depletion; presence of cytogenetic abnormalities), or
provided relevant auxiliary information regarding the missing values (year of
transplantation; time between diagnosis and transplantation; and whether
diagnosis was primary or secondary myelofibrosis). Note that since this is a
model for (complementary log-log transformed) cumulative incidence of relapse,
we want to make sure to include predictors known to be associated with the
cause-specific hazards of both relapse and non-relapse mortality.

Since around 45% of patients were either event-free or censored within the
first 60 months (see supplementary material S2.2, non-parametric curves),
potential censoring times for those experiencing non-relapse mortality were
first multiply imputed using the {kmi} package in strata defined by
(completely observed) year of transplantation, yielding 100 datasets with
“complete” subdistribution time $V$ but with partially observed covariate
information. In each of these datasets, covariates were imputed once using
each of the four imputation methods used in the simulation study, after 20
cycles across the covariates. The choice of 100 imputed datasets was motivated
using von Hippel’s quadratic rule (i.e. number of imputed datasets needed
should increase approximately quadratically with increasing fraction of
missing information), based on an initial set of 30 imputed datasets 39.
Essentially, we sought to control the MCSEs of the standard errors of the
estimated subdistribution log hazard ratios. Default imputation methods were
used depending on the type of covariate: binary covariates using logistic
regression, ordered categorical using proportional odds regression and nominal
categorical using multinomial logistic regression. For continuous covariates,
the default in {mice} is predictive mean matching, while linear regression is
used for $f(X_{j}\,|\,X_{-j},Z;\psi)$ in {smcfcs}. The imputation model for a
given partially observed variable therefore contained as predictors all
remaining fully and partially observed variables from the substantive model,
together with the outcome. Each imputation approach differs mainly in how they
incorporate the outcome in the imputation model: either by sampling directly
from an assumed substantive model compatible distribution (FG-SMC and CS-SMC),
or by including event indicator(s) and marginal cumulative hazard(s)
explicitly as additional predictors (FG-Approx and CS-Approx).

Figure 7 shows for all methods the estimated baseline cumulative incidence
function, and the width of the corresponding confidence interval at each
timepoint. As was the case in the simulation study, cumulative incidences are
estimated in each imputed dataset, and pooled after complementary log-log
transformation. The estimation procedure used for the standard errors of the
cumulative incidences is described by Ozenne et al.40. The estimates using
both FG-SMC and FG-Approx are virtually overlapping, which is consistent with
the simulation study results when $p=0.15$. Both CS-SMC and CS-Approx also
yielded cumulative incidences that were close to those obtained by the
subdistribution hazard based imputation approaches, which is in line with the
results of the simulation study under random right censoring. The most stark
differences were between CCA (which only uses 20% of patients) and the
imputation approaches: the cumulative incidence of relapse at 60 months was
almost 5% lower than the nearest MI-based curve, with confidence intervals
that were over twice as wide. For completeness, in supplementary material S2.3
we report the pooled subdistribution log hazard ratios, in addition to the
pooled coefficients of cause-specific Cox models for relapse and non-relapse
mortality (each containing the same predictors as the Fine–Gray model for
relapse). The pooled coefficients of the Fine–Gray models were extremely
similar between imputation approaches, and all differed considerably from the
(much more variable) CCA. There were some noticeable differences between
subdistribution hazard based and cause-specific hazard based imputation
approaches when estimating the cause-specific Cox model for non-relapse
mortality (see e.g. pooled coefficients for weight loss prior to
transplantation, hemoglobin or high risk comorbidity score). Furthermore, the
pooled subdistribution log hazard ratios were generally small in magnitude
(none exceeding 0.5), a setting in which both SMC and approximately compatible
approaches are expected to perform similarly.

The differences observed between point estimates obtained using the imputation
based approaches and CCA are in large part explainable by the gulf in
efficiency between the two approaches. Nevertheless, there are indications
that the estimates obtained using imputation methods would be less biased than
their CCA counterparts in this example. An exploratory logistic model showed
that the observed time to competing event and competing event indicator were
both predictive of the probability of being a complete-case, after adjusting
for other known important predictors of missingness such as year of
transplantation (many variables recorded more often later on in time as their
clinical relevance became clearer). Upon closer inspection, it appears that
the probability of being a complete-case is significantly lower only for those
censored earlier on in time. This seemingly unlikely association between
future outcome and baseline complete-case indicator (outcome-dependent MAR,
under which CCA is biased) is likely confounded by transplant centre. That is,
shorter follow-up times and missing values in covariates may both be
symptomatic of a given centre’s overall quality of data collection. Although
ignored in the present analysis for simplicity, there is indeed heterogeneity
in data completeness between EBMT affiliated transplant centres across and
within different countries. The MI of potential censoring times would allow to
model centre effects using standard software, for example by means of
stratification or use of a frailty term.

Figure 7: Pooled baseline cumulative incidence functions for relapse in the
applied data example (left panel), and width of corresponding confidence
intervals (right panel). These are the estimates for a patient aged 60,
transplanted in 2019 immediately after diagnosis, with 10g/dL hemoglobin, 15
x109/L white blood cells, no peripheral blood blasts, and with reference
levels for all categorical predictors (see supplementary Table 1).

##  6 Discussion

In this paper, we extended the SMC-FCS approach in order to impute missing
covariates compatibly with a Fine–Gray substantive model. For a given
competing event, the theory relies on using the subdistribution time $V$ and
the corresponding event-specific indicator as outcome variables. In the
presence of random right-censoring, $V$ is only partially observed, as the
potential censoring times for those failing from competing events are unknown.
These can be multiply imputed in a first step, after which covariates can be
imputed by conditioning on the “complete” outcome variables. The approach is
straightforward to implement in R by making use of existing software packages
{kmi} and {smcfcs}. While the imputation of potential censoring times appears
underused in the subdistribution hazard modelling literature (relative to
weighted approaches), it has inspired other methodological extensions e.g.
enabling the use of deep learning in discrete time after single imputation of
potential censoring times41.

The simulation study compared the performance of the proposed method to
competing MI approaches, including imputing compatibly with cause-specific
proportional hazards models. The FG-SMC approach performed optimally (in terms
of estimating both subdistribution log hazard ratios, and cumulative
incidences) when the assumption of proportional subdistribution hazards held,
and performed satisfactorily when this assumption did not hold. For cumulative
incidence estimation, the choice of substantive model (i.e. cause-specific Cox
vs. Fine–Gray) at the analysis phase appears to be more important than the
procedure used to impute the missing covariates. In terms of RMSE of these
predictions, most imputation approaches outperform CCA. The applied data
example also demonstrated the possible gain in efficiency when using MI
instead of CCA.

One counterintuitive finding was that the presence of censoring seems to
improve the performance of the misspecified SMC-FCS procedure (e.g. use of CS-
SMC when underlying DGM assumes proportional subdistribution hazards). An
explanation for this phenomenon is that the time-dependent factor relating the
cause-specific and subdistribution hazards for cause 1 (the “reduction
factor”15) is closer to 1 earlier in time. Therefore (in the example with DGM
assuming proportional subdistribution hazards), the violation of
proportionality on the cause-specific hazard scale will appear to be less
severe in earlier time-periods, thereby improving the performance of the
misspecified SMC-FCS approach. This is also in line with earlier findings
showing how similar the results of subdistribution and cause-specific hazards
models can be in presence of heavy censoring28, 42.

An advantage to the proposed approach is that it can be extended in various
ways. First, while not yet possible using existing {smcfcs} software, the
censoring times can depend on partially observed $X$ (rather than just on
complete $Z$), meaning $X$ and $V$ will need to be imputed iteratively. The
FG-Approx method with censoring depending on $X$ can be easily implemented in
{mice} using custom imputation methods. Second, the approach can account for
time-dependent effects, by making direct use of existing approaches developed
in the context of standard Cox models43. Third, the proposed approach can be
extended to accommodate interval censored outcomes, using the methodology
described by Delord and Génin, which relies on analogous principles: multiply
impute interval censored $V$ in order to work with simpler censoring complete
data44.

There are multiple limitations to the present work. The first is that the
proposed SMC-FCS approach does not accommodate delayed entry (left
truncation). Our current recommendation to impute approximately compatibly
with a Fine–Gray model subject to delayed entry and right-censoring is to
include $I(D=1)$ and $\hat{\Lambda}_{1}(T)$ as predictors in the imputation
model, in addition to other substantive model covariates. Here,
$\hat{\Lambda}_{1}(t)$ is the estimated cumulative subdistribution hazard
based on a marginal model that uses time-dependent weights in order to
accommodate both left-truncation and right-censoring 45. Note the proposed
imputation model uses $\hat{\Lambda}_{1}(T)$ and not $\hat{\Lambda}_{1}(V)$,
and therefore some downward bias is to be expected, as explained in Appendix
A. Second, while FG-SMC does not require an explicit model for the competing
risks, it does require the censoring distribution to be specified explicitly
(e.g. non-parametrically using KM, or using a Cox model). Third, the proposed
approach is geared towards imputing missing covariates when only one competing
event is of interest. More generally, the strategy of estimating a Fine–Gray
for each cause in turn is not an approach the current authors endorse, based
on both theoretical 14, 17 and simulation-based arguments 18. When multiple
competing events are of interest, we would instead recommend modelling the
cause-specific hazards, or using the semiparametric approach suggested by Mao
and Lin for joint inference on the cumulative incidence functions46.

In conclusion, the proposed approach is most appropriate for imputing missing
covariates in the context of prognostic modelling of only one event of
interest. Based on the simulation study, imputing compatibly with cause-
specific proportional hazards seems to be a good all-round strategy for a
“complete” competing risks analysis (investigating both the cause-specific
hazards and cumulative incidence functions47), and can at the same time be
used for prognostic modelling based on the cause-specific Cox models.

## Software

All R code (needed to reproduce simulation study, applied data example, and
manuscript figures) is available at https://github.com/survival-
lumc/FineGrayCovarMI. In addition to the minimal R code provided in the
supplementary materials, a wrapper function for the proposed SMC-FCS Fine–Gray
method is available inside the {smcfcs} R package.

## Acknowledgements

The authors would like to thank the Chronic Malignancies Working Party (CMWP)
of the EBMT for allowing the use of the myelofibrosis dataset, as well as the
patients and centres involved in the original study. We also thank Linda
Koster (EBMT Leiden Study Unit) for continued support with this dataset. This
work was supported by Grant BE 4500/4-1 from the German Research Foundation
(JB), Grant W222158-2-50 from the Leiden University Fund (EFB), and UK MRC
grant MR/T023953/2 (JWB).

## Conflicts of interest

The authors declare no conflicts of interest.

## References

  * 1 Carpenter JR, Smuk M. Missing Data: A Statistical Framework for Practice. Biom J. 2021;63(5):915–947.  doi: 10.1002/bimj.202000196
  * 2 Sterne JAC, White IR, Carlin JB, et al. Multiple Imputation for Missing Data in Epidemiological and Clinical Research: Potential and Pitfalls. BMJ. 2009;338:b2393.  doi: 10.1136/bmj.b2393
  * 3 van Buuren S, Groothuis-Oudshoorn K. Mice: Multivariate Imputation by Chained Equations in R. J Stat Softw. 2011;45(3):1–67.  doi: 10.18637/jss.v045.i03
  * 4 Prentice RL, Kalbfleisch JD, Peterson AV, Flournoy N, Farewell VT, Breslow NE. The Analysis of Failure Times in the Presence of Competing Risks. Biometrics. 1978;34(4):541–554.  doi: 10.2307/2530374
  * 5 Bonneville EF, Resche-Rigon M, Schetelig J, Putter H, de Wreede LC. Multiple Imputation for Cause-Specific Cox Models: Assessing Methods for Estimation and Prediction. Stat Methods Med Res. 2022;31(10):1860–1880.  doi: 10.1177/09622802221102623
  * 6 White IR, Royston P. Imputing Missing Covariate Values for the Cox Model. Stat Med. 2009;28(15):1982–1998.  doi: 10.1002/sim.3618
  * 7 Bartlett JW, Seaman SR, White IR, Carpenter JR. Multiple Imputation of Covariates by Fully Conditional Specification: Accommodating the Substantive Model. Stat Methods Med Res. 2015;24(4):462–487.  doi: 10.1177/0962280214521348
  * 8 Bartlett JW, Taylor JMG. Missing Covariates in Competing Risks Analysis. Biostatistics. 2016;17(4):751–763.  doi: 10.1093/biostatistics/kxw019
  * 9 Fine JP, Gray RJ. A Proportional Hazards Model for the Subdistribution of a Competing Risk. J Am Stat Assoc. 1999;94(446):496–509.  doi: 10.1080/01621459.1999.10474144
  * 10 Lau B, Lesko C. Missingness in the Setting of Competing Risks: From Missing Values to Missing Potential Outcomes. Curr Epidemiol Rep. 2018;5(2):153–159.  doi: 10.1007/s40471-018-0142-3
  * 11 Archer L, Koshiaris C, Lay-Flurrie S, et al. Development and External Validation of a Risk Prediction Model for Falls in Patients with an Indication for Antihypertensive Treatment: Retrospective Cohort Study. BMJ. 2022;379:e070918.  doi: 10.1136/bmj-2022-070918
  * 12 Clift AK, Coupland CAC, Keogh RH, et al. Living Risk Prediction Algorithm (QCOVID) for Risk of Hospital Admission and Mortality from Coronavirus 19 in Adults: National Derivation and Validation Cohort Study. BMJ. 2020;371:m3731.  doi: 10.1136/bmj.m3731
  * 13 Ruan PK, Gray RJ. Analyses of Cumulative Incidence Functions via Non-Parametric Multiple Imputation. Stat Med. 2008;27(27):5709–5724.  doi: 10.1002/sim.3402
  * 14 Beyersmann J, Allignol A, Schumacher M. Competing Risks and Multistate Models with R.  Use R!New York: Springer-Verlag, 2012
  * 15 Putter H, Schumacher M, van Houwelingen HC. On the Relation between the Cause-Specific Hazard and the Subdistribution Rate for Competing Risks Data: The Fine–Gray Model Revisited. Biom J. 2020;62(3):790–807.  doi: 10.1002/bimj.201800274
  * 16 Gray RJ. A Class of K-Sample Tests for Comparing the Cumulative Incidence of a Competing Risk. Ann Stat. 1988;16(3):1141–1154.
  * 17 Austin PC, Steyerberg EW, Putter H. Fine-Gray Subdistribution Hazard Models to Simultaneously Estimate the Absolute Risk of Different Event Types: Cumulative Total Failure Probability May Exceed 1. Stat Med. 2021;40(19):4200–4212.  doi: 10.1002/sim.9023
  * 18 Bonneville EF, de Wreede LC, Putter H. Why You Should Avoid Using Multiple Fine–Gray Models: Insights from (Attempts at) Simulating Proportional Subdistribution Hazards Data. J R Stat Soc Series A Stat Soc. 2024 (in press).
  * 19 Taylor JMG, Murray S, Hsu CH. Survival Estimation and Testing via Multiple Imputation. Statistics & probability letters. 2002;58(3):221–232.
  * 20 Rubin D. Multiple Imputation for Nonresponse in Surveys.  New York: John Wiley & Sons, 1987
  * 21 Molenberghs G, Fitzmaurice G, Kenward MG, Tsiatis A, Verbeke G. Handbook of Missing Data Methodology.  CRC Press, 2014.
  * 22 Allignol A, Beyersmann J. Software for Fitting Nonstandard Proportional Subdistribution Hazards Models. Biostatistics. 2010;11(4):674–675.  doi: 10.1093/biostatistics/kxq018
  * 23 Bartlett JW, Keogh RH, Bonneville EF. Smcfcs: Multiple Imputation of Covariates by Substantive Model Compatible Fully Conditional Specification. 2022\.
  * 24 Lambert PC, Wilkes SR, Crowther MJ. Flexible Parametric Modelling of the Cause-Specific Cumulative Incidence Function. Stat Med. 2017;36(9):1429–1446.  doi: 10.1002/sim.7208
  * 25 Morris TP, White IR, Crowther MJ. Using Simulation Studies to Evaluate Statistical Methods. Stat Med. 2019;38(11):2074–2102.  doi: 10.1002/sim.8086
  * 26 Klein J, Moeschberger M. Survival Analysis: Techniques for Censored and Truncated Data (2nd Edition).  New York: Springer-Verlag, 2003.
  * 27 Beyersmann J, Latouche A, Buchholz A, Schumacher M. Simulating Competing Risks Data in Survival Analysis. Stat Med. 2009;28(6):956–971.  doi: 10.1002/sim.3516
  * 28 Grambauer N, Schumacher M, Beyersmann J. Proportional Subdistribution Hazards Modeling Offers a Summary Analysis, Even If Misspecified. Stat Med. 2010;29(7-8):875–884.  doi: 10.1002/sim.3786
  * 29 van Buuren S. Flexible Imputation of Missing Data. Second Edition.  Boca Raton, FL: Chapman & Hall/CRC Press, 2018.
  * 30 Morisot A, Bessaoud F, Landais P, Rébillard X, Trétarre B, Daurès JP. Prostate Cancer: Net Survival and Cause-Specific Survival Rates after Multiple Imputation. BMC Med Res Methodol. 2015;15(1):54.  doi: 10.1186/s12874-015-0048-4
  * 31 Marshall A, Altman DG, Holder RL, Royston P. Combining Estimates of Interest in Prognostic Modelling Studies after Multiple Imputation: Current Practice and Guidelines. BMC Med Res Methodol. 2009;9(1):57.  doi: 10.1186/1471-2288-9-57
  * 32 Wood AM, Royston P, White IR. The Estimation and Use of Predictions for the Assessment of Model Performance Using Large Samples with Multiply Imputed Data. Biom J. 2015;57(4):614–632.  doi: 10.1002/bimj.201400004
  * 33 Mertens BJA, Banzato E, de Wreede LC. Construction and Assessment of Prediction Rules for Binary Outcome in the Presence of Missing Predictor Data Using Multiple Imputation and Cross-Validation: Methodological Approach and Data-Based Evaluation. Biom J. 2020;62(3):724–741.  doi: 10.1002/bimj.201800289
  * 34 R Core Team . R: A Language and Environment for Statistical Computing. Vienna, Austria: 2023.
  * 35 Therneau TM. A Package for Survival Analysis in R. 2023.
  * 36 Gasparini A. Rsimsum: Summarise Results from Monte Carlo Simulation Studies. J Open Source Softw. 2018;3(26):739.  doi: 10.21105/joss.00739
  * 37 Polverelli N, Bonneville EF, de Wreede LC, et al. Impact of Comorbidities and Body Mass Index on the Outcomes of Allogeneic Hematopoietic Cell Transplantation in Myelofibrosis: A Study on Behalf of the Chronic Malignancies Working Party of EBMT. Am J Hematol. 2024;99(5):993—996.  doi: 10.1002/ajh.27262
  * 38 Kröger N, Bacigalupo A, Barbui T, et al. Indication and Management of Allogeneic Haematopoietic Stem-Cell Transplantation in Myelofibrosis: Updated Recommendations by the EBMT/ELN International Working Group. Lancet Haematol. 2024;11(1):e62-e74.  doi: 10.1016/S2352-3026(23)00305-8
  * 39 von Hippel PT. How Many Imputations Do You Need? A Two-stage Calculation Using a Quadratic Rule. Sociol Methods Res. 2020;49(3):699–718.  doi: 10.1177/0049124117747303
  * 40 Ozenne B, Sørensen AL, Scheike T, Torp-Pedersen C, Gerds TA. riskRegression: Predicting the Risk of an Event Using Cox Regression Models. The R Journal. 2017;9(2):440–460.  doi: 10.32614/RJ-2017-062
  * 41 Gorgi Zadeh S, Behning C, Schmid M. An Imputation Approach Using Subdistribution Weights for Deep Survival Analysis with Competing Events. Sci Rep. 2022;12(1):3815.  doi: 10.1038/s41598-022-07828-7
  * 42 van der Pas S, Nelissen R, Fiocco M. Different Competing Risks Models for Different Questions May Give Similar Results in Arthroplasty Registers in the Presence of Few Events. Acta Orthop. 2018;89(2):145–151.  doi: 10.1080/17453674.2018.1427314
  * 43 Keogh RH, Morris TP. Multiple Imputation in Cox Regression When There Are Time-Varying Effects of Covariates. Stat Med. 2018;37(25):3661–3678.  doi: 10.1002/sim.7842
  * 44 Delord M, Génin E. Multiple Imputation for Competing Risks Regression with Interval-Censored Data. J Stat Comput Simul. 2016;86(11):2217–2228.  doi: 10.1080/00949655.2015.1106543
  * 45 Geskus RB. Cause-Specific Cumulative Incidence Estimation and the Fine and Gray Model Under Both Left Truncation and Right Censoring. Biometrics. 2011;67(1):39–49.  doi: 10.1111/j.1541-0420.2010.01420.x
  * 46 Mao L, Lin DY. Efficient Estimation of Semiparametric Transformation Models for the Cumulative Incidence of Competing Risks. J R Stat Soc Series B Stat Methodol. 2017;79(2):573–587.  doi: 10.1111/rssb.12177
  * 47 Latouche A, Allignol A, Beyersmann J, Labopin M, Fine JP. A Competing Risks Analysis Should Report Results on All Cause-Specific Hazards and Cumulative Incidence Functions. J Clin Epidemiol. 2013;66(6):648–653.  doi: 10.1016/j.jclinepi.2012.09.017

##  Appendix A Imputed censoring times, and resulting cumulative
subdistribution hazards

As described in Section 3.3.1, the subdistribution time $V$ is only partially
observed in the presence of random right-censoring. Thus, the potential
censoring times for those failing from cause 2 should first be multiply
imputed, before imputing any missing covariates. This imputation of partially
observed $V$ is visualised more closely in Figure 8, using a simulated dataset
of 2000 individuals following the parametrisation used in the simulation study
scenario with correctly specified Fine–Gray, $p=0.65$, and random exponential
censoring. In this example, the potential censoring times for those failing
from cause 2 were imputed $m=10$ times.

The upper panel shows the imputed potential censoring times for a random
selection of 20 individuals failing from cause 2, in addition to their cause 2
failure time and their true eventual censoring time. The lower panel shows the
estimated marginal cumulative subdistribution hazard function for
$\hat{\Lambda}_{1}(t)$ resulting from using $I(D=1)$ together with either the
imputed or true $V$ as outcomes in a marginal model. We used
$\hat{\Lambda}_{1}(t)$ estimated using the true $V$ to create the secondary
x-axis in the upper panel, which shows the value of this function at a given
timepoint. For example, the marginal cumulative subdistribution hazard was
1.146 at timepoint 2.5, and stayed constant at 1.322 after the last cause 1
event in this sample.

The upper panel in particular gives additional insights regarding the FG-
Approx method, where $I(D=1)$ and $\hat{\Lambda}_{1}(V)$ are included as
predictors in the imputation model. Namely, the secondary x-axis shows the
value of $\hat{\Lambda}_{1}(V)$ used in the imputation model for a missing
$X_{j}$, for given imputed $V$. A first key point is that one should always
use $\hat{\Lambda}_{1}(V)$ in the imputation model, and not
$\hat{\Lambda}_{1}(T)$. Since the observed cause 2 failure time occurs before
the eventual censoring time, $\hat{\Lambda}_{1}(T)$ will always be smaller
than the marginal cumulative subdistribution hazard at the eventual censoring
time, and so the imputation model will incur some downward bias. A second
point is that in settings with fewer event 1 failures (e.g. $p=0.15$ scenario
in the simulation study), the corresponding secondary x-axis will have a
smaller range, since the subdistribution hazard will be lower overall. The
imputed potential censoring times will therefore have a relatively smaller
influence on the imputed $X_{j}$.

The lower panel shows that the estimated $\hat{\Lambda}_{1}(t)$ varies very
little between imputed datasets, with differences only being noticeable later
on in follow-up as risk sets become smaller and associated cumulative hazard
jumps more pronounced. Note also that while $\hat{\Lambda}_{1}(t)$ based on
the true $V$ appears in this dataset to be a kind of “average” of the
functions based on imputed $V$, this will not be the case in general,
especially with smaller sample sizes. The $\hat{\Lambda}_{1}(t)$ based on the
weighted estimator45 will however coincide with the “average” of the functions
based on imputed $V$, as will using the negative log of one minus the Aalen-
Johansen estimate of the marginal cumulative incidence function.

Figure 8: Based on a simulated dataset of $n=2000$ (correctly specified
Fine–Gray, $p=0.65$, random censoring), we show the imputed ($m=10$
imputations) potential censoring times for a random selection of 20
individuals failing from cause 2 (upper panel); and the estimated marginal
cumulative subdistribution hazard function for cause 1 based on true $V$, and
based on imputed $V$ (lower panel).