A review of the use of time-varying covariates in the Fine-Gray subdistribution hazard competing risk regression model.

In survival analysis, time-varying covariates are covariates whose value can change during follow-up. Outcomes in medical research are frequently subject to competing risks (events precluding the occurrence of the primary outcome). We review the types of time-varying covariates and highlight the effect of their inclusion in the subdistribution hazard model. External time-dependent covariates are external to the subject, can effect the failure process, but are not otherwise involved in the failure mechanism. Internal time-varying covariates are measured on the subject, can effect the failure process directly, and may also be impacted by the failure mechanism. In the absence of competing risks, a consequence of including internal time-dependent covariates in the Cox model is that one cannot estimate the survival function or the effect of covariates on the survival function. In the presence of competing risks, the inclusion of internal time-varying covariates in a subdistribution hazard model results in the loss of the ability to estimate the cumulative incidence function (CIF) or the effect of covariates on the CIF. Furthermore, the definition of the risk set for the subdistribution hazard function can make defining internal time-varying covariates difficult or impossible. We conducted a review of the use of time-varying covariates in subdistribution hazard models in articles published in the medical literature in 2015 and in the first 5 months of 2019. Seven percent of articles published included a time-varying covariate. Several inappropriately described a time-varying covariate as having an association with the risk of the outcome.

INTRODUCTION

Survival analysis is concerned with outcomes that occur over time. Examples include time to all‐cause death and time to cause‐specific death (eg, time to death due to cardiovascular causes). A key feature of survival analysis is that of censoring: the study may be terminated or subjects may be lost to follow‐up before the event of interest has been observed to occur for all subjects. Another important concept is that of competing risks. A competing risk is an outcome whose occurrence precludes the occurrence of the event of interest. For instance, subjects who die of noncardiovascular causes are no longer at risk of death due to cardiovascular causes.

Time‐varying (or time‐dependent) covariates occur frequently in biomedical and epidemiological research. These are covariates whose values change over the duration of follow‐up. Examples include biomarkers that change over follow‐up (eg, prostate specific antigen [PSA]) or cumulative exposure to medications.1 A time‐invariant (or time‐fixed) covariate is a covariate whose value remains unchanged or fixed throughout the entire duration of follow‐up.

The use of hazard‐based regression models is ubiquitous in medical and epidemiological research. The most common of these is the Cox proportional hazards regression model that allows one to estimate the effect of covariates on the hazard of the occurrence of the outcome.2 In the absence of competing risks and with only time‐invariant covariates, one is able to estimate the probability of the occurrence of the event within any duration of time from the fitted Cox proportional hazards model. However, this may no longer be possible when the model incorporates time‐varying covariates.

In the presence of competing risks with either time‐varying or time‐invariant covariates, the Cox proportional hazards regression model for the cause‐specific hazard function does not permit estimation of the effect of covariates on the cumulative incidence of the outcome or estimation of the probability of the event occurring within a given duration of follow‐up time. The Fine‐Gray subdistribution hazard regression model permits estimation of the effect of time‐invariant covariates on the cumulative incidence of the event in the presence of competing risks.3 Similarly to the Cox proportional hazards model without competing risks, the Fine‐Gray subdistribution hazard model must be interpreted carefully when time‐varying covariates are included in the model specification. In many scenarios, one loses the ability to estimate the effect of the covariate on the cumulative incidence function (CIF) and to estimate the probability of the event occurring in a given time period.

The objectives of the current paper are twofold. First, to describe the different types of time‐dependent covariates and clarify when hazard‐based regression models can be interpreted in terms of the cumulative incidence, both with and without competing risks. Special emphasis is given to the Fine‐Gray model. Second, to review the medical literature to determine the frequency with which time‐varying covariates are incorporated into Fine‐Gray subdistribution hazard regression models and how authors interpret the resultant models.

TIME‐VARYING COVARIATES AND CUMULATIVE INCIDENCE

In this section, we introduce notation and definitions. We then summarize existing knowledge on the effect of including time‐varying covariates in the conventional Cox model and the Fine‐Gray subdistribution hazard model. We discuss the effect of their inclusion on the interpretation of the resultant regression coefficients.

Definitions and notation

The survival function, S(t), is the probability that an individual survives to time t: S(, where T denotes the time to the occurrence of the event of interest. The complement of the survival function is the cumulative distribution function, which describes the cumulative incidence of the event of interest up to time t: F(.

The hazard function, λ(t), is the instantaneous rate of the occurrence of the event of interest in subjects who are currently at risk of the event: . The survival function and the hazard function are related through the following relationship: , where denotes the cumulative hazard function (alternatively, Λ().4

The Cox proportional hazards regression model allows one to estimate the relative effect of covariates on the hazard function.2 The Cox model can be formulated in multiplicative form as , where denotes a vector of covariates, denotes a vector of regression coefficients, and denotes the baseline hazard function for a subject whose covariates are all equal to zero.

A time‐varying covariate, X(t), is a covariate whose value can change over the duration of follow‐up (eg, time‐varying biomarkers [eg, PSA], current use of medication, and cumulative dose of medication). Kalbfleisch and Prentice distinguish between two different types of time‐varying covariates: external and internal covariates.5 They define these as “an external covariate is one that is not directly involved with the failure mechanism, and an internal covariate is a time measurement that is taken on the individual” (page 123). External covariates can effect the failure process directly, but are not otherwise involved in the failure mechanism, while internal covariates can effect the failure process, but can also be impacted by the failure mechanism (eg, there may be feedback between the time‐varying covariate and failure process). There are two subfamilies of external covariates: external defined covariates and external ancillary covariates. External defined covariates are those covariates whose values are defined or determined for all subjects prior to the commencement of the study. An example would be a controlled experiment studying the effect of cumulative radiation dose, where the (cumulative) radiation dose, as a function of time, is known and defined prior to the beginning of the study. Thus, for external defined covariates, the function X(t) would be defined for all subjects and for all times t, regardless of whether the subject was currently under observation. An external ancillary covariate is a covariate that is the result of a stochastic or random process that is external to the subjects in the study. Kalbfleisch and Prentice suggest that airborne pollution would be an external ancillary covariate in a study on asthma attacks, as the probability distribution of airborne pollution is likely independent of subject characteristics that would be included in the hazard regression model. In contrast to an external covariate, Kalbfleisch and Prentice define an internal covariate to be “the output of a stochastic process that is generated by the individual under study and so is observed only so long as the individual survives and is uncensored. In consequence, its observed value carries information about the survival time of the corresponding individual” (page 124). It can be argued that the large majority of time‐varying covariates encountered in biomedical research are internal time‐varying covariates. Examples include organ transplantation, time‐varying biomarkers (eg, PSA), presence of infection, current use of medication in observational studies, and cumulative dose of medication in observational studies where the dose is not actively manipulated by the investigator according to a previously defined mechanism.

It is well known that the Cox proportional hazards model can incorporate time‐varying covariates, with the usual formulation being: . As with the conventional Cox model with time‐invariant covariates, the inclusion of time‐varying covariates allows one to estimate the effect of covariates on the instantaneous hazard of the event. The model can also incorporate time‐varying covariate effects (ie, the effect of a given covariate on the instantaneous hazard of the event changes over time). However, we do not consider time‐varying covariate effects in this paper.

Relationship of hazard and survival functions in the absence of competing risks and with time‐invariant covariates

Assume that we are in a setting without competing risks and such that all covariates are fixed at baseline. From the conventional Cox model for the hazard function, we can derive the following relationship to the survival function: (see Appendix). Applying a log‐log transformation, we have that log(log(. Thus, the regression coefficients from the Cox regression model can be interpreted as having an effect on the survival function. In particular, if a covariate increases the hazard of the occurrence of the event (β > 0), it will also increase the incidence of the event.6

Survival and hazards in the absence of competing risks and with time‐varying covariates

In the previous section, we demonstrated that the survival function could be estimated as a function of the baseline survival function and the linear predictor from the Cox proportional hazards model in the absence of competing risk and with time‐invariant covariates. Thus, the regression coefficients can be interpreted as having an effect on the survival function. However, with time‐varying covariates, this interpretation only holds when the hazard regression model includes external time‐varying covariates and not with internal time‐varying covariates. To see how this arises more clearly, as above, without covariates, we have that the survival function is related to the cumulative hazard function: . In the absence of competing risks, the CIF is equal to 1 − exp(−Λ(. Given the Cox regression model with time‐varying covariates, we have that Note that, in the above integral, the term exp( cannot be brought outside of the integral because it is a function of the variable of integration. This complicates the interpretation of the regression parameters in the proportional hazards model.

There are two different scenarios that we must consider. The first scenario is when the time‐varying covariate is an external defined covariate. In this setting, the usual relationship between the hazard and the CIF holds, such that (note that, in the presence of time‐varying covariates, the survival function of interest is defined formally as ; however, we use the simpler notation when it is clear from the context what is intended). For external time‐varying covariates, X(t) is known even when a subject is not under observation and the above integral can be evaluated and has a clear probabilistic interpretation. However, one cannot bring the term exp( outside of the integral and have it multiply the resultant cumulative baseline hazard function. Consequently, in general, we can no longer make simple claims that a covariate that has an effect on the hazard of the outcome has an effect of the same direction on the cumulative incidence of the outcome, as it depends on the entire history of the time‐varying covariate. Greater care is needed in such statements. In comparing two individuals with covariates and , with for at all t, then implies that the cumulative incidence for individual 2 with covariate is greater than for individual 1 with covariate at all time points, and vice versa if . When and are not strictly ordered, it is challenging to compare the CIFs for the two subjects, which depends not only on and the history of the time‐varying covariates but also the baseline hazard function.

The second scenario is when the time‐varying covariate is an internal covariate. As such there is almost never a defined stochastic model or known functional form for the covariate. Furthermore, even if this were known, it is questionable whether it would be valid when the individual is no longer under observation. In particular, once the individual has experienced the event, the value of the covariate is no longer known. Thus, in this setting, the required integral cannot be evaluated and the usual relationship between the hazard and survival function no longer exists. Thus, one cannot estimate the effects of internal time‐varying covariates on the survival function via a proportional hazards model. Accordingly, one cannot obtain an estimate of the survival function or of the CIF.

Cumulative incidence and hazards in the presence of competing risks and with time‐invariant covariates

As shown in Section 2.2, in the absence of competing risks, there is a one‐to‐one correspondence between the effect of a covariate on the hazard function and its effect on the survival function. In particular, a covariate that increases the hazard of an event will also increase the incidence of the event. In the presence of competing risks, this is no longer necessarily true.3 Specifically, when one naively treats competing events as censoring events and implicitly fits a model to the cause‐specific hazard function for the event of interest. To address this issue, Fine and Gray introduced a proportional hazards model for the subdistribution hazard function for the kth type of event where D is a variable denoting the type of event that occurred. There is a different subdistribution hazard function for each of the K different types of events. The subdistribution hazard function for a given type of event is defined as the instantaneous rate of occurrence of the given type of event in subjects who have not yet experienced an event of that type.3 Note that the risk set consists of those subjects who are either currently event‐free or who have previously experienced a competing event. In subsequent discussions, it will be important to recall that subjects who experience a competing event remain in the risk set.

The Fine‐Gray subdistribution hazard model estimates the effect of covariates on the subdistribution hazard function via the model , where denotes the baseline subdistribution hazard function for the kth event type. In the absence of time‐varying covariates, there is a one‐to‐one relationship between the subdistribution hazard function and the CIF, which describes the incidence of the occurrence of an event while taking competing risks into account. Thus, the subdistribution hazard model allows one to estimate the effect of covariates on the CIF for the event of interest. In particular, one may recover a relationship similar in form to that described above where CIF, and where CIF denotes the baseline CIF: CIF. By using the log‐log transformation, this can be written as a linear model: log (log(1 − CIF. Thus, if a covariate is associated with an increase in the subdistribution hazard function, it will also be associated with an increase in the cumulative incidence of the event. Thus, in the absence of time‐varying covariates, the use of the Fine‐Gray subdistribution hazard model allows for inferences about the effect of a covariate on the incidence of the outcome.

Incidence and hazards in the presence of competing risks with time‐varying covariates

In the original paper developing the subdistribution hazard regression model, the model was extended to include time‐varying covariates. However, inclusion of time‐varying covariates requires caution in all instances and can be problematic in some. Furthermore, the resultant regression coefficients must be interpreted carefully.

We define the CIF in the presence of time‐varying covariates as CIF. Then, the subdistribution hazard regression model with time‐varying covariates can be written as . Consequently, we have that CIF .3

Note that in the above integral, the term exp( cannot be brought outside of the integral because it is a function of the variable of integration. This complicates the interpretation of the regression parameters in the proportional subdistribution hazard model.

There are two different scenarios that we must consider. The first scenario is when the time‐varying covariate is an external defined covariate. In this setting, the usual relationship between the subdistribution hazard and the CIF holds, such that . For external time‐varying covariates, X(t) is known even when a subject is not under observation and the above integral can be evaluated and has a clear probabilistic meaning. However, one cannot bring the term exp( outside of the integral. Consequently, in general, we can no longer make simple claims that a covariate that has an effect on the subdistribution hazard of the outcome has an effect of the same direction on the cumulative incidence of the outcome, as it depends on the entire history of the time‐varying covariate. We emphasize that the evaluation of the integral is always possible for external time‐varying covariate. In practice, an internal time‐varying covariate may not be fully observed, complicating the evaluation of the integral. This is similar to the conventional Cox model without competing risks, where the evaluation of the integral requires that the time‐dependent covariate is completely observed over time. It follows heuristically that the connection between the subdistribution hazard and the CIF is only valid with external time‐dependent covariates and the interpretation of the Fine‐Gray model with internal time‐dependent covariates is only valid for the subdistribution hazard and not the CIF. Such interpretation is addressed more carefully below.

A different problem presents itself in the setting of internal time‐varying covariates. Several sets of authors have noted that an unconventionality in the definition of the risk set for the subdistribution hazard function makes the inclusion of internal time‐varying covariates problematic.7, 8, 9 As noted above, the risk set for the subdistribution hazard function contains those subjects who have failed due to competing events. Incorporating time‐varying covariates into the Fine‐Gray subdistribution hazard model requires that we know the value of these covariates for the entire time that a subject remains in the risk set. When a subject has failed from a competing event, the model requires that the value of the time‐varying covariate can be fully specified over time. For competing events that include death (or cause‐specific death), one will, in general, not know the value of the time‐varying covariate once the subject has died. For instance, consider the setting in which the event of interest is death due to cardiovascular causes and death due to noncardiovascular causes is the competing event. Elevated blood pressure (hypertension) is a known risk factor for cardiovascular disease.10 To include blood pressure as a time‐varying covariate would require knowing a subject's blood pressure after he or she had experienced a competing event (death due to noncardiovascular causes), as he or she remains in the risk set. However, it is not clear how blood pressure should be defined after noncardiovascular death. For this reason, defining internal time‐dependent covariates after the occurrence of the competing event (noncardiovascular death) may be difficult or impossible. This issue has been considered in the literature, for example, in the work of Beyersmann and Schumacher,8 where ad hoc approaches to extrapolating the internal time‐dependent covariate, such as last value carried forward, have been proposed. The use of such techniques leads to an implicit definition of the subdistribution hazard that may be difficult to interpret and should be used with great caution.

LITERATURE REVIEW OF THE USE OF TIME‐VARYING COVARIATES WITH THE FINE‐GRAY SUBDISTRIBUTION HAZARD MODEL

In the previous section, we discussed the inclusion of time‐varying covariates in the Fine‐Gray competing risk regression model. In this section, we report on a literature review that examined the frequency with which time‐varying covariates were included in the subdistribution hazard model in two different periods in the medical literature. We also report on how the resultant model was interpreted.

We used a search strategy based on that used in a recently‐published review examining how authors in the biomedical literature interpreted the regression coefficients from a Fine‐Gray regression model.11 The search was conducted in two distinct time periods. We searched the PubMed database (https://www.ncbi.nlm.nih.gov/pubmed) using two search strategies to identify papers that used the Fine‐Gray subdistribution hazard model: (i) (“subdistribution hazard”[All Fields] OR “Fine‐Gray”[All Fields]) AND (“2015/01/01”[PDAT]: “2015/12/31”[PDAT]) and (ii) (“subdistribution hazard”[All Fields] OR “Fine‐Gray”[All Fields]) AND (“2019/01/01”[PDAT]: “2019/05/31”[PDAT]). The first strategy identified papers published in 2015, while the second identified papers published in the first 5 months of 2019. The first search strategy was identical to that which we used in an earlier review.11 The second search strategy was conducted on June 6, 2019.

The first search process identified 64 papers published in 2015. Of these, eight methodologically‐oriented publications were excluded. One paper was excluded because it did not use the Fine‐Gray subdistribution hazard model. We examined the remaining 55 papers to see whether the authors incorporated time‐varying covariates into the Fine‐Gray subdistribution hazard model. The second search process identified 108 papers published in the first 5 months of 2019. Of these, five methodologically‐oriented publications were excluded, and one paper was excluded because an English‐language version was not available. We examined the remaining 102 papers to see whether the authors incorporated time‐varying covariates into the Fine‐Gray subdistribution hazard model.

Of the 55 papers published in 2015, six (11%) included time‐varying covariates in Fine‐Gray subdistribution hazard model.12, 13, 14, 15, 16, 17 Of the 102 papers published in the first 5 months of 2019, five (5%) included time‐varying covariates in the Fine‐Gray subdistribution hazard model.18, 19, 20, 21, 22 These 11 papers are summarized in Table 1.

Table 1

Papers reporting time‐varying covariates in a Fine‐Gray subdistribution hazard model

Author	Primary outcome	Competing risk	Binary time‐	Continuous time‐
			varying covariate	varying covariate
Han et al12	Cervical cancer‐	Death from other		Cumulative dose
	specific mortality	causes		of metformin
				after cervical
				cancer diagnosis
Ong et al13	Death	Discharge from	The occurrence
		the ICU	of a specific
			infection
Dautzenberg et al15	ICU mortality	ICU discharge	Colonization by
			a specific organism
Vajdic et al16	Cause‐specific	Deaths from	Specific	Age, calendar
	mortality	other causes	infections and	year of follow‐
			undergoing	up, and number
			opioid	of episodes of
			substitution	opioid
			therapy	substitution
				therapy
Ong et al14	Mortality	Liberation from	Viral reactivation
		mechanical	status
		ventilator
		support
Klein Hesselink et al17	Atrial fibrillation	All‐cause	Use of neck	Cumulative
		mortality	radiotherapy	thyroid
				stimulating
				hormone level
				and radioiodine dose
Aksnessaether et al18	Second cancer	Death	Treatment
Li et al19	Acute coronary	Death	Antiviral
	syndrome/end‐		treatment
	stage renal
	disease/ischaemic
	stroke/retinopathy
Deka et al20	Development of	Death	Androgen
	depression		deprivation
			therapy
Mori et al21	TCZ	TCZ	Clinical Disease
	discontinuation	discontinuation	Activity Index
	due to secondary	due to remission	(as binary or
	loss of efficacy	or other reasons	categorical
			variable)
Yaffe et al22	Dementia	Death		Number of
				follow‐up visits
				per year

Across the 11 studies, all but two of the time‐varying covariates were internal covariates. Age and calendar year of follow‐up can be thought of as external time‐varying covariates, as they can be fully specified at all time points after baseline, regardless of whether the subject had experienced a competing event. Several of the covariates (eg, cumulative dose of metformin, undergoing opioid substitution therapy, number of episodes of opioid substitution therapy, use of neck radiotherapy, cumulative thyroid stimulating hormone level, and radioiodine dose) could be external time‐varying covariates in a controlled experimental study in which these variables were dictated by a protocol that was pre‐specified and under the control of the investigator. However, the studies in which these time‐varying covariates were used were retrospective observational studies in which the study investigators passively recorded the treatments and exposures that had been applied. Accordingly, these variables are all internal time‐varying covariates. Nine of the 11 studies had at least one binary (or categorical) internal time‐varying covariate. It was difficult to determine how the internal time‐dependent covariates were coded for subjects who experienced a competing event.

Of the 11 studies that incorporated time‐varying covariates, six (55%) provided, in at least one instance, an interpretation that suggested that the time‐varying covariate was associated with the risk of the event.12, 13, 17, 18, 19, 21 As risk has a distinct probabilistic interpretation, it appears that these authors were suggesting that the time‐varying covariate was associated with the CIF. In the remaining papers, one set of authors explicitly noted that “cumulative incidence functions cannot be interpreted for time‐varying risk factors”15 (page 1172). In the papers that did not infer an association between the time‐varying covariate and the risk of the event, the authors restricted themselves to the reporting of subdistribution hazard ratios or describing the effect of the time‐varying covariate in terms of its association with the rate or hazard of the outcome. These results are important, as they provide information as to how authors think their results should be interpreted, regardless of whether the interpretation is valid.

DISCUSSION

The Fine‐Gray subdistribution hazard model is increasingly being used for the analysis of time‐to‐event outcomes in the presence of competing events. However, many applied analysts and clinical researchers appear to be unaware that incorporating time‐dependent covariates in the Fine‐Gray model requires considerable care and that their inclusion may have undesirable effects on the interpretability of the resultant model.

The focus of the current paper is on the use of time‐varying covariates with the Fine‐Gray subdistribution hazard model. However, similar issues arise with the conventional Cox proportional hazards model in the absence of competing risks. When including time‐varying covariates in a conventional Cox regression model, one can still estimate the effect of covariates on the hazard function. However, one may lose the ability to estimate the effect of covariates on the survival function and to make inferences about the direction of the effect of the covariate on the incidence of the outcome (see Section 2.3).

Several sets of authors have highlighted how the unconventionality of the risk set for the Fine‐Gray subdistribution hazard model makes the inclusion of internal or subject‐level time‐varying covariates problematic.7, 8, 9 The risk set contains those subjects who have experienced a competing event. However, their inclusion in the risk set requires that the value of the time‐dependent covariate be known for these subjects after they have experienced a competing event. If the competing event includes death or cause‐specific death, then it will typically not be possible to know the value of the time‐varying covariate after the subject has died and ad hoc approaches for extrapolation are needed, complicating the interpretation of the model. Our review of the use of the Fine‐Gray subdistribution hazard model in the medical literature highlights that, despite this fact being known to those who are experienced with competing risks regression models, it appears to be unfamiliar to many applied analysts and clinical researchers. We observed that the inclusion of internal time‐varying covariates in a Fine‐Gray subdistribution hazard model was not uncommon in the medical literature, with the presentation being unclear and potentially inappropriate.

Beyersmann and Schumacher used the relationship between discrete covariates and multistate models to develop an approach to incorporate time‐dependent covariates in a subdistribution hazards framework.8 However, in a recent study, Poguntke et al used simulations to examine the impact of including time‐varying covariates in subdistribution hazard models.23 When they simulated data under a setting in which the time‐varying covariate had no association with the incidence of the primary outcome, they detected a non‐null subdistribution hazard ratio. Consequently, they strongly recommended “avoiding the use of the subdistribution approach for assessing the effect of a time‐dependent covariate”23 (page 9). However, this latter paper did not distinguish between internal and external time‐dependent covariates and their conclusions should be considered to apply to the setting with internal time‐dependent covariates. Importantly, for external covariates, valid inferences about the CIF may be obtained from the subdistribution hazard model, similar to data without competing risks where the survival function may be obtained from the proportional hazards model.

When the focus is on prediction or estimation of the CIF in the presence of internal time‐varying covariates, Cortese et al describe how the Fine‐Gray model can be combined with a landmark approach to estimate cumulative incidence.24 This approach is based on a similar approach described by van Houwelingen for use with the all‐cause hazard model and the cause‐specific hazard model.25 This approach entails defining a set of k landmark times: {. One then conducts a sequence of conventional competing risks survival analyses with time‐invariant covariates at each of the k landmark times. The ith landmark analysis is restricted to those subjects who are at risk of the event at the ith landmark time (ie, those subjects who have not experienced the primary outcome of interest or a competing event prior to time s ). The covariates for the ith landmark analysis are fixed at their values at time s . Thus, in the ith landmark analysis, one is estimating Pr(. The vector X(s) denotes the values of the covariates at time s and is treated as a time‐invariant covariate in this analysis. Thus, a conventional Fine‐Gray subdistribution hazard model can be used to model the incidence of the outcome conditional on survival to time s and on the values of the covariates at time s. When prediction of incidence or survival probabilities is of interest, this approach may provide easier and more transparent interpretations with internal time‐dependent covariates for survival data both with and without competing risks. It should be noted that the use of landmarking comes at the cost of increasing variability (or decreasing precision) of the estimates with increasing landmark time (or as the size of the risk set decreases).

Even when a well‐defined time‐dependent covariate is included in a Fine‐Gray subdistribution hazard model, unanticipated consequences may arise of which the applied analyst or clinical researcher may be unaware. An appealing feature of the Fine‐Gray subdistribution hazard model, and indeed the motivation for its development, is the ability to make inferences about the effect of covariates on the CIF in the presence of competing risks (which cannot be done with the conventional Cox regression model for the cause‐specific hazard function). However, this capacity is lost when incorporating internal time‐varying covariates, similarly to the proportional hazards model without competing risks. Thus, while the inclusion of such time‐varying covariates in the Fine‐Gray model permits one to interpret the covariates as having an effect on the subdistribution hazard of the primary event, that is, on the rate of the event in those subjects who have not yet experienced that event, it no longer permits inferences about the association of a covariate with the CIF or the incidence of the event. As we have noted elsewhere, this interpretation can be difficult for many to conceptualize11 and is less appealing than an interpretation about the effect of the covariate on the CIF. Furthermore, the study by Poguntke et al suggests that the estimated subdistribution hazard ratio may be misleading when the time‐dependent covariate is an internal time‐dependent covariate.

van Walraven and McAlister examined 100 studies published in 2013 in prominent medical journals.26 Of these 100 studies, 46% contained Kaplan‐Meier analyses that were susceptible to competing risks bias. Similarly, Schumacher et al conducted an examination of the presence of competing risk bias in 2015 in a prominent medical journal.27 In 51 articles in which competing risks were present, 25 (49%) were susceptible to competing risk bias (eg, a Kaplan‐Meier survival function was estimated rather than a CIF). We recently published two reviews addressing issues related to the handling of competing risks in the medical literature. The first review examined how competing risks were addressed in reports of randomized controlled trials (RCTs) published in four leading general medical journals.28 We estimated that 77.5% of RCTs with a time‐to‐event outcome were potentially susceptible to competing risks. Of those studies that were potentially susceptible to competing risks, 77.4% reported the results of a Kaplan‐Meier survival analysis, while only 16.1% reported using CIFs to estimate the incidence of the outcome over time in the presence of competing risks. The second review examined how authors interpreted the hazard ratios arising from a Fine‐Gray subdistribution hazard model.11 We found that many authors had an incorrect or unclear interpretation of the estimated subdistribution hazard ratios. In that review, we provided guidelines for interpreting the subdistribution hazard ratio resulting from a Fine‐Gray competing risk regression model. The focus of the current review was on how authors in the biomedical literature incorporated time‐varying covariates in these models.

In our current review, we found that 7% of studies that used the Fine‐Gray subdistribution hazard model included time‐varying covariates in the model. In most instances, these were internal time‐varying covariates. As noted above, the inclusion of internal time‐varying covariates can be problematic due to the unconventional definition of the risk set that contains those who have experienced a competing event. In half of these studies, the competing event was either all‐cause mortality or cause‐specific mortality. It is unclear how the time‐varying covariates were defined following the subject's death from a competing cause. Putting aside the concern about defining the time‐varying covariate in subjects who have experienced a competing event, as we have previously discussed, the inclusion of internal time‐dependent covariates implies that the resultant model is no longer a model for the CIF. Thus, one cannot make claims about the association of the covariate with increased risk of the event. Despite this limitation, half of the studies that included internal time‐varying covariates interpreted a time‐varying covariate as having an effect on the risk of the occurrence of the event of interest.

Statistical software packages differ in their ability to incorporate time‐varying covariates in the Fine‐Gray subdistribution hazard regression model. The crr function in the cmprsk package for R permits the inclusion of continuous time‐varying covariates that are interactions of time‐invariant covariates and known polynomial functions of time. Thus, it does not permit the inclusion of binary or categorical time‐varying covariates (these models can also be fit using functions from the riskRegression package, which provide a formula interface similar to that in the survival package). In contrast to this, both the stcrreg function in Stata and the PHREG procedure in SAS permit the inclusion of both continuous and categorical time‐varying covariates. Regardless of the capabilities of a given statistical software package, one must exercise extreme caution in the interpretation of the resultant model.

We summarize our recommendations as follows. First, external time‐dependent covariates may be included in the Fine‐Gray subdistribution hazard model. However, analysts must interpret the resultant model carefully, as the estimated regression coefficients no longer describe the association of the covariate with the CIF. Thus, the fitted model loses its original attractive feature: the ability to directly estimate the association between covariates and the CIF. The interpretation of the covariate effects on the CIF requires careful consideration, similarly to the proportional hazards model with external time‐dependent covariates in the absence of competing risks. For a fixed realization of the external time‐varying covariate, the CIF may be calculated via either analytic or numerical integration, permitting a comparison of the risks for subjects with different time‐varying covariates. For the case of ordered covariates, there is a simple ordering of the CIFs based on the sign of the regression coefficient in the proportional subdistribution hazard regression model.

Second, analysts should exercise extreme caution in including internal time‐varying covariates in the Fine‐Gray subdistribution hazard model. The nature of the subdistribution hazard risk set can make defining the covariates difficult for subjects who have experienced a competing event. Furthermore, one can no longer make inferences (either directly or via integration) about the effect of the covariate on the CIF. Finally, based on the findings of Poguntke et al, there is the risk that the estimated subdistribution hazard ratio based on the assumed covariate definition may be misleading, yielding unclear inferences. While we are reluctant to suggest that internal time‐varying should never be included in a Fine‐Gray subdistribution hazard model, these arguments suggest that their inclusion should be rare and that analysts should proceed with utmost caution.