Understanding Landmarking and Its Relation with Time-Dependent Cox Regression

Time-dependent Cox regression and landmarking are the two most commonly used approaches for the analysis of time-dependent covariates in time-to-event data. The estimated effect of the time-dependent covariate in a landmarking analysis is based on the value of the time-dependent covariate at the landmark time point, after which the time-dependent covariate may change value. In this note we derive expressions for the (time-varying) regression coefficient of the time-dependent covariate in the landmark analysis, in terms of the regression coefficient and baseline hazard of the time-dependent Cox regression. These relations are illustrated using simulation studies and using the Stanford heart transplant data.


Introduction
Time-dependent covariates play an important role in the analysis of censored time-toevent data. Prominent examples include the effect of heart transplant on survival for heart patients [6] and the effect of CD4+ T-cell counts on the occurrence of AIDS or death for HIV-infected patients [15]. In the first example the time-dependent covariate is a binary covariate, in the second example a numerical covariate, typically mea-sured longitudinally. These two examples constitute the most common instances of time-dependent covariates in survival analysis. Examples in the context of organ transplantation include the aforementioned heart transplant example, but also the effect of kidney transplantation or changing dialysis modality in end-stage renal disease, or the effect of graft failure on survival for patients with a liver transplant.
Broadly speaking, two approaches have come in general use in estimating the effect of time-dependent covariates. The first is time-dependent Cox regression, already mentioned by [5]. In this approach the hazard at time t is assumed to depend on the current value at time t of the time-dependent covariate, X (t), through the product of a baseline hazard and exp(β X (t)). This approach yields valid inference if the value of the time-dependent covariate is known for all subjects at all event time points without error, and the regression model is correctly specified. Especially in the longitudinal setting, this is typically not the case and many researchers have studied the amount of bias when measurement error and ageing of covariates are present [2,11,15]. A second approach is landmarking [3], which involves setting a landmark time point s, and using the value of the time-dependent covariate at s (or using some other appropriate summary of the history of the time-dependent covariate up to s) as a timefixed covariate in an analysis of survival from s onwards, in a subset of subjects at risk at s. For overviews we refer to [7,12].
In principle, both approaches can be used for estimating the effects of timedependent covariates, and there are no clear settings where the choice between time-dependent Cox and landmarking is obvious. Each of the two methods has its advantages and disadvantages. To get some feeling of the relative merits of the two approaches, let us consider the situation of the Stanford heart transplant example [6]. Patients are admitted to a waiting list, the time to event is time from admittance to the waiting list until death, which may be subject to censoring, and interest is in the effect of a heart transplant on survival. The time-dependent covariate heart transplant, X (t), is thus initially equal to 0, and attains the value 1 as soon as the patient receives a heart transplant. (If the patient never receives a heart transplant, the value remains 0 throughout his/her follow-up.) The most important reason for the popularity of landmarking, especially in the present context of a binary time-dependent covariate, is its transparency. It is clear what is being compared: at the landmark time s, two groups are compared with regard to their survival from time s onwards, one group without, the other group with a heart transplant received before or at time s. Differences between these two groups can be visualized by plotting the Kaplan-Meier survival curves for the two groups. In contrast, such a visualization is much less obvious for the time-dependent Cox model. Survival curves can only be shown for patients with X (t) ≡ 0 and X (t) ≡ 1, respectively. Model-free curves have also been proposed in this context, sometimes referred to as Simon-Makuch curves [13]. Both model-based and Simon-Makuch curves show the survival for fictional patients who either never receive a heart transplant (X (t) ≡ 0) or have received a heart transplant at t = 0 (X (t) ≡ 1). In both cases one may question whether these curves reflect clinically realistic quantities. The landmarking curves have a much clearer interpretation; see also [4] for a detailed discussion on this topic. Disadvantages of landmarking are first of all the need for a (to some extent arbitrary) choice of landmark time point, and also a loss of power because, especially for later landmark time points, subjects with an event before the landmark time point are excluded from analysis. In many cases, early landmark time points also lead to a loss of power, because in the beginning of follow-up, there will be few subjects with a heart transplant, leading to highly unbalanced groups. For estimation of the effect of a heart transplant on survival, both methods are equally applicable. The topic of this paper is the following: suppose that a time-dependent Cox regression model is valid. In a landmark analysis at landmark time s two groups, one with (X (s) = 1) and one without (X (s) = 0) a heart transplant are compared. Those subjects without a heart transplant at time s might in the future receive one, so that X (t) = 1 for some future t, whereas for subjects with a heart transplant at time s, X (t) will keep the value 1. This will make the groups with X (s) = 0 and X (s) = 1 more similar in the future and will therefore attenuate the effect of the time-dependent covariate, compared to the time-dependent Cox regression.
The purpose of this paper is to understand, quantify and illustrate the differences between the regression coefficients of a time-dependent Cox model and those obtained in a landmark analysis. Starting from a time-dependent Cox regression model, which we assume to be correctly specified, we will derive formulas for the (time-varying) regression coefficient corresponding to the landmark model. We study two special cases of time-dependent covariates in more detail and show that if the time-dependent Cox model satisfies the proportional hazards assumption, there will be attenuation in the sense that the landmark regression coefficient is between the time-dependent Cox regression coefficient and 0. We show that the degree of attenuation depends on the rate of change of the time-dependent covariate. An illustration using the Stanford heart transplant data is provided.

Theory
Let T denote the time-to-event random variable of interest. Let X (t) denote a timedependent covariate, which for simplicity we take to be one-dimensional, and denote its complete history until time t by X (t) = {X (s); 0 ≤ s ≤ t}. We assume that the hazard of T , conditional on X (t), is given by i.e. the hazard at time t is assumed to depend on the whole history X (t) only through the present value, and it follows a Cox model with possibly time-varying effect given by β(t). The relations to be derived in this paper are valid irrespective of censoring, but for consistent estimation of the parameters in the model and survival predictions it is assumed that censoring is independent of T and X (·), possibly given other time-fixed covariates in the model, which have been omitted here for the sake of simplicity.
Fix a time point s. A landmark analysis at time s will fix the value of X (·) at X (s), and assume a model for T with X (s) as time-fixed covariate, based on the survivors at risk at time s [3]. Independent censoring, and, if appropriate, truncation, is assumed, which implies that the survivors at risk at time s are representative for survivors at time s. As a result the hazard considered for the landmark model at landmark time s is given by where the additional conditioning on T ≥ s is in fact superfluous. It is only retained in the notation of h LM (t | s, X (s)) to emphasize that the landmark analysis is based on survivors at s. The postulated model for this landmark hazard is typically taken to be a proportional hazards model as well: Often β LM (t | s) is taken to be time-fixed in the analysis, i.e. β LM (t | s) ≡ β LM (s), but it may (and typically will) depend on s. The question addressed in this note is as follows: what is the relation between β LM (t | s) and β(t), and how does it depend on h 0 (t) and on the development of X (t)?
Intuitively, the landmark model employs an old value, X (s), instead of the current value X (t) to describe the hazard at time t. As a result, if X (t) changes rapidly between X (s) and X (t) and if X (t) is strongly related to T , one may expect a large discrepancy between β(t) and β LM (t | s). The aim is to quantify how quickly β LM (t | s) changes for t ≥ s, depending on the rate of change of X (t), on β(t) and on h 0 (t).
Our development will start by considering the conditional survival function at time t > s, given survival until time s and given X (s), The landmark hazard h LM (t | s, X (s)) is the derivative with respect to t of the negative logarithm of S(t | s, X (s)) at t = s + w, provided that, conditional on T ≥ t and X (t), T is independent of X (s). This expression can also be found in [15]. It is possible to further evaluate h LM (t | s, X (s)) for a number of specific models for the development of X (t). In what follows we consider two of such models. The first of these is the common situation where X (t) is dichotomous, the other is a specific case where X (t) is continuous.

Dichotomous Time-Dependent Covariates
Let X (t) be a dichotomous time-dependent covariate. This is the type of situation for which the landmarking approach was originally proposed [3]. In that paper the endpoint was overall survival, and interest was in the effect of response to chemotherapy (a timedependent covariate, coded as 0 = no response, 1 = response). Many more examples can be given, including the effect of disease recurrence on survival [16], the effect of treatment adherence on disease recurrence [8] or the effect of adverse events on recurrence rates [9].

Theory
In this context our aim is to obtain an expression of The conditional expectation can be evaluated by considering the multi-state model given in Fig. 1.
States 0 and 1 correspond to the values of the time-dependent covariate (response) being 0 and 1, respectively, and state 2 is the death state. The original landmarking paper did not consider a possible transition back from state 1 to 0, but we will allow this here. The relation between the event time T and the multi-state model X (t) is given by the equivalence In the present context, Eq. (3) is valid if for those alive at time t, conditional on the current state X (t), the transition to death is independent of X (s), the state at time s. This assumption is fulfilled if the multi-state model is Markov. If that is the case, then the transition probabilities P gh (s, t) = P(X (t) = h | X (s) = g) can be calculated if the hazards of making a transition from 0 to 1 or backwards and the transitions from 0 or 1 to state 2 (death) are known, using the Kolmogorov-Chapman forward equations. The transition hazard from g to h will be denoted by λ gh (t). The conditional prevalence probabilities defined for g = 0, 1, describe the conditional probability of being in response (in state 1) at time t, given in state g at the earlier landmark time s, and given alive at time t. Finally, we define π g0 (s, t) = 1 − π g1 (s, t). The conditional expectation in (3) can then be written as This implies that the hazard ratio in the landmark model can be expressed as A number of remarks can be made regarding this expression. First, if β(t) ≡ 0, that is, when the time-dependent covariate has no effect at all on survival, then we have β LM (t | s) ≡ 0, independently of π g0 (s, t) and π g1 (s, t). Second, the landmark regression coefficient β LM (t | s) is always between −β(t) and β(t). A further simplification can be made when considering the most common situation, the irreversible case, where the 1 → 0 transition is not possible, and hence π 11 (s, t) ≡ 1. This gives and has 0 ≤ β LM (t | s) ≤ β(t). The intuitive explanation of the formula is as follows: if those with X (s) = 0 quickly jump to 1 (so if π 01 (s, t) is high), then the effect is quickly attenuated, so β LM (t | s) is close to 0, while if those with X (s) = 0 remain in 0, then the effect is not attenuated, so β LM (t | s) will be close to β(t). The derivative with respect to t of β LM (t | s), evaluated at the landmark s, is given by If the time-dependent Cox model is correct and satisfies the proportional hazards assumption, i.e. if β(t) ≡ β, then β LM (t | s) will usually vary over t, unless for instance β = 0. As mentioned earlier, usually the effect of the time-dependent covariate is estimated through a proportional hazards model like (1), where the proportional hazards assumption would ignore the possibly time-varying nature of β LM (t | s). Equations (4) and (5) show that a proportional hazards landmark model would typically be misspecified. If such a misspecified Cox regression landmark model is fitted, then [17][18][19] the estimate obtained in that model will be approximately equal to Figure 2 shows a plot of β LM (t | s) of Eq. (5) on the y-axis in the irreversible case, where both time to response and time to death follow exponential distributions with different values of the rates λ 01 (t) ≡ ρ for response, and β; the death rate without response was set to λ 02 (t) ≡ λ = 0.1. The landmark regression coefficients β LM (t | s) have been recalculated for different values of s. The attenuation as time between s and t increases can be clearly seen. The degree of attenuation increases with larger values of ρ.

Illustration
The time-dependent effect of β LM (t | s) is further illustrated by generating a single large dataset (n = 10, 000) based on the same exponential distributions as before with ρ = 0.2, λ = 0.1, β = 1. A landmark analysis was performed at s = 2, with death as endpoint and X (s), response at 2 years, as time-fixed covariate. Subsequently the method of [10] based on Schoenfeld residuals, as implemented in cox.zph in the survival package [14] in R was applied. This method gives an approximation of β LM (t | s) as function of t. Figure 3 shows the result, on the left the plot including the residuals, on the right only the (default) lowess curve in black showing the approximation ofβ LM (t | s), with 95 % confidence intervals in black dashed lines. In dotted lines the theoretical formula of (4) is shown. Time  The conclusion is that the theoretical time-varying effect of β LM (t | s) can be detected and retrieved (the lowess curve being close to the theoretical curve) in a large dataset.
After focusing on the time-varying nature of β LM (t | s) of the landmark model, we now turn to what is being estimated when (as is usually done in practice) a proportional hazards landmark model is fitted to the data. Figure 4 shows box-plots of the estimates of β * LM (s) of (6), obtained from 1000 simulations from datasets of 10,000 subjects, with ρ = 0.2, λ = 0.1, β = 1, for values of s = 2, 4, 6, 8, 10, together with the time-dependent Cox regression estimate. In Fig. 4a no censoring was applied, while in Fig. 4b independent uniform censoring between 7.5 and 12.5 was applied.
The mean of the estimates of β * LM (s) hardly changed with s (a minimal decrease from 0.504 at s = 2 to 0.476 at s = 10). The mean for the time-dependent Cox regression was 1.001. With censoring, the mean of the estimates of β * LM (s) increased from 0.558 at s = 2 to 0.808 at s = 10. The reason for this increase compared to the case of no censoring is that for later landmark times, contributions from β LM (t | s) for

Theory
Explicit calculations as for the case of dichotomous time-dependent covariates are much harder to obtain for continuous time-dependent covariates. Recall from (2) We will follow [2,15] in assuming that conditioning on T ≥ t is negligible [11] (see also), and that the joint distribution of X (s) and X (t) is Gaussian. The results given below should thus be understood as approximations that should be reasonable in low-risk situations. We will derive such approximations for the case where the timedependent covariate X (t) follows a Gaussian process with mean μ(t) = EX (t) and covariance function K (s, t) = cov(X (s), X (t)). Then the conditional distribution of 1 (s, s). The conditional expectation on the right-hand side of (7) is in fact the moment generating function of the conditional distribution of X (t), given X (s), evaluated at β(t) (recalling that the additional condition that T ≥ t is ignored). This leads to the approximation X (s)) . (8) Note again that the approximation in (8) is due to having ignored the selection on having survived to time t. A quite useful special case, also considered in [2], is the case where the timedependent covariate of subject i at time t follows with the mean μ(t) a deterministic time trend, b i a random person effect, following a zero mean normal distribution with variance ω 2 , and X * i (t) a mean zero Ornstein-Uhlenbeck (OU) process, starting at X * i (0) = 0, and defined further by where W i (t) is a Wiener process and θ and σ are parameters describing the degree of mean reversal (to zero) and influence of the random fluctuations of the Wiener process, respectively. The solution of (10) is given by [1] (see, e.g.) [A.4] This is a stationary Wiener process with covariances Adding the random person effect b, the result is a stationary Wiener process with where σ 2 tot = ω 2 + σ 2 /(2θ) is the total variance of X (t) and ρ = ω 2 /σ 2 tot is the intraclass correlation, the proportion of the total variance represented by the random person effect variance. The correlation drives the conditional distribution of X i (t), given X i (s), which can be seen to be normal with mean μ(t | X i (s)) and variance σ 2 (t | X i (s)), given by Using (8), this leads to the approximation t) . (12) Taking logarithms in (12) and differentiating with respect to X (s) yields So, when the time-dependent covariate follows (9)-(10), the regression coefficient β LM (t | s) in the landmark model is approximately equal to the regression coefficient β(t) in the time-dependent Cox model, attenuated by the correlation function ρ(s, t), defined in (11). The correlation function ρ(s, t), defined in (11), decreases exponentially from 1 to ρ as t increases from s to ∞, the speed of decay depending on θ .

Illustration
Four plots in Fig. 6 illustrate the approximation in (12) for a relatively simple situation with one time-dependent covariate. The baseline hazard h 0 (t) was taken to be constant, equal to 0.1, corresponding to an exponential distribution with mean 10. The hazard at time t, given the history X (t), was taken to be with β = 0.5. Figure 6a shows the reference situation where the trajectories of the time-dependent covariate were generated according to an OU process, with μ(t) ≡ 0, total variance σ 2 tot = 0.5, intraclass correlation ρ = 0.5 and θ = 1. For given σ 2 tot , ρ and θ , the random person effect variance ω 2 was chosen such that ω 2 /σ 2 tot = ρ, and σ 2 such that ω 2 + σ 2 /(2θ) = σ 2 tot . The landmark was fixed at s = 1. The approximation in (12) is contrasted with a Monte Carlo approximation where an OU process with the parameters specified above was generated, jointly with an event time following (14). For each of a very large number of subjects, i = 1, . . . , N , first a random effect b i was generated according to a mean zero normal distribution with variance ω 2 . The starting value of the realization of X * i (t) was set to X * i (0) = 0. Subsequently, for a series of very short intervals of length dt (we took 0.01), given the current value of X * i (t), the value of X i (t) was set at b i + X * i (t), the hazard was calculated according to (14), and a coin was tossed with probability h 0 (t) exp(β X i (t)) dt to decide if the subject would die in that interval. If not, a new value for X * i (t + dt) was calculated according to X * i (t) − θ X * i (t) dt + σ U dt, with U an independent standard normal random variable, and X i (t + dt) was set as b i + X * i (t + dt). The procedure for subject i stopped when the subject died or at a censoring time of 10 years. A very large database of processes and associated death or censoring times was stored. For a given landmark time point s, here 1 year, the number of evaluable subjects (i.e. surviving until s) was set at 500,000. Conditioning on X (s) = 0 was achieved by considering only those simulated processes for which X (s) was less than dX away from 0 (we took dX = 0.1). Finally, E e β X (t) | T ≥ t, X (s) = 0 for fixed t > s was approximated by considering the further subset of processes for which T ≥ t and calculating the average value of e β X (t) within those subsets. The procedure of taking appropriate subsets (from the same database) was repeated for X (s) = ±0.5. These latter approximations are shown as the logarithm in solid (wiggly, because of Monte Carlo error) lines, the approximations in (12) (also logarithm) as dotted lines. Figures 6b-d show results when the reference situation is changed by choosing different values of θ (b), ρ (c) and σ 2 tot (d). All curves in Fig. 6a-d start at β X (s) = 0.25, 0, −0.25, for X (s) = 0.5, 0, −0.5, respectively, at t = s. The parameters θ and ρ determine the speed of attenuation for t > s and the asymptotes for t → ∞ of the regression coefficients β LM (t | s) in the landmark analyses. Lower values of θ and ρ imply a higher degree of attenuation. The second term 1 2 β 2 σ 2 tot (1 − ρ 2 (s, t)) increases from 0 at t = s to 1 2 β 2 σ 2 tot at t → ∞ and is incorporated equally in each of the dotted curves with each plot, since it does not depend on X (s). Changing the value of σ 2 tot [comparing (a) and (d)] only leads to vertical shifts of the dotted curves and does not change the relative distance for different values of X (s). Finally, the difference between the theoretical approximations shown in the dotted lines and the Monte Carlo approximations in the solid lines reflect the effect of selective removal (because of death), which was ignored in (8) and onwards. Compared to a situation with no removal of subjects because of death, for a given X (s), subjects with higher X (t) have a higher probability of being removed. As a result, subjects with lower X (t) remain in the population, which results in the dotted curves being pulled downwards. The total variance σ 2 tot , which was seen not to influence the theoretical β LM (t | s) in (13), does influence this selective removal. This behaviour is quite similar to the selective removal of subjects with high frailty values in frailty models, the effect of which is also stronger with increasing frailty variance. Further simulation studies (not shown here) indicated that the effect of selective removal is stronger (i.e. the approximation of (12) less accurate) when the baseline death rate is increased, a phenomenon that is also present in frailty models.

Data Illustration
We further illustrate our results using the well-known Stanford heart transplant data [6], consisting of 103 patients admitted to a waiting list for a heart transplant. The event time is time from admittance to the waiting list until death; interest is in the effect of a heart transplant on survival. Of the 103 patients, 69 received a heart transplant, and a total of 75 patients died, 45 with a heart transplant and 30 without a heart transplant. Median follow-up calculated by reverse Kaplan-Meier was 2.51 years. An important covariate predictive of the effect of heart transplant is the mismatch score, measured for those patients with a heart transplant. It is a continuous score derived from antibody responses of pregnant women [6] and reflects the degree of incompatibility (based on tissue typing) between the donor and recipient. Median mismatch score was 1.08, with 0.75 and 1.58 as 25th and 75th percentiles. Because different effects of the heart transplant may be expected for patients with a high mismatch score and patients with a low mismatch score, we distinguish between patients with a mismatch score in the highest quartile and the rest. Four patients with a heart transplant and no mismatch score (because no tissue typing was performed) are included in the larger set of patients with a mismatch score ≤ 1.58. We define two time-dependent covariates of the type studied in Sect. 2.1: X 1 (t) = 1, if the patient has received a heart transplant before time t with a mismatch score higher than 1.58 and 0 otherwise, and X 2 (t) = 1, if the patient has received a heart transplant before time t with a mismatch score less than or equal to 1.58 and 0 otherwise. A time-dependent Cox regression with X 1 (t) and X 2 (t) as time-dependent covariates resulted in estimated regression coefficients of 0.605 with a standard error (SE) of 0.386 (hazard ratio; 95 % confidence interval 1.83; 0.86-3.90) for X 1 (t), and −0.031 with a standard error (SE) of 0.319 (hazard ratio; 95 % confidence interval 0.97; 0.52-1.81) for X 2 (t). The tentative conclusion, which has to be seen in the light of the fact that heart transplantation was in its infancy during data collection, seems to be that heart transplants with a high mismatch score do more harm than good and that no clear effect can be seen of heart transplants with a low/medium mismatch score. Since differences between time-dependent Cox regression and landmarking can be most clearly seen for time-dependent covariates with large effects, we will focus on the effect of X 1 (t) in a number of landmark models.
Most of the heart transplants in the data occur in the first couple of months, so we take s = 1, 1.5, 2 months as landmark time points for illustration. Earlier and later landmark time points result in numbers of transplant with (mismatch) score > 1.58 of less than ten. Table 1 gathers the results of the time-dependent Cox regression and of the three landmark analyses.
It is clear that, compared to the hazard ratio of 1.83 of the time-dependent Cox regression, the estimated hazard ratios are indeed attenuated towards 1, as predicted by our results. The degree of attenuation is determined by many factors: π 01 (s, t) in Eq. (4), possible time-varying effects β(t) in the time-dependent Cox model, and by the degree of censoring [see Eq. (6)], which makes it hard to compare the results of the different landmark models.

Discussion
In this paper we derived relations between the regression coefficients obtained in a landmark analysis and those of a time-dependent Cox regression, when interest is in the effect of a time-dependent covariate on survival. In case the time-dependent covariate has no effect on survival at all, i.e. when the time-dependent Cox regression coefficient is identically 0, the landmark regression coefficient is identically 0 as well. Otherwise the time-dependent Cox regression coefficient will be attenuated. Different formulas apply for dichotomous and continuous covariates, but the degree of attenuation is mainly determined by how quickly the value of the time-dependent covariate changes over time. For dichotomous time-dependent covariates this is expressed in the prevalence probabilities, while for the Ornstein-Uhlenbeck example it is expressed in terms of the intraclass correlation and the θ parameter describing the degree of mean reversal.
We did not study the effects of measurement error (misclassification error in the case of dichotomous time-dependent covariates) or ageing due to infrequent measurements. The approximations in Sect. 2.2.1 can be adapted to account for that, in the spirit of [2], and they will result in a further attenuation of the regression coefficient of the landmark analysis. The main reason for not considering this aspect in detail is that landmark analysis and time-dependent Cox regression analysis will both be affected by these complications.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.