\documentclass{article} \usepackage{amsmath} \usepackage{Sweave} \addtolength{\textwidth}{1in} \addtolength{\oddsidemargin}{-.5in} \setlength{\evensidemargin}{\oddsidemargin} %\VignetteIndexEntry{Using Time Dependent Covariates} \title{Using Time Dependent Covariates and Time Dependent Coefficients in the Cox Model} \author{Terry Therneau \and Cindy Crowson\\ Mayo Clinic} \begin{document} \maketitle \SweaveOpts{keep.source=TRUE} <>= options(width=60, continue=" ") makefig <- function(file, top=1, right=1, left=4) { pdf(file, width=9.5, height=7, pointsize=18) par(mar=c(4, left, top, right) +.1) } library(survival) @ \section{Introduction} One of the strengths of the Cox model is its ability to encompass covariates that change over time, due to the theoretical foundation in martigales. A \emph{martingale} (original definition) is a betting strategy in games of chance. One of the simplest and best known is doubling the bet each time you lose. For instance consider the following game of roulette: \begin{equation*} \begin{tabular}{rcc|c} Bet & Outcome & Win & Running total \\ \hline R \$1 & Red & 2 & 1 \\ &&&\\ R \$1 & Black & 0 & 0 \\ R \$2 & Black & 0 & -2 \\ B \$4 & Red & 0 & -6 \\ R \$8 & Black & 0 & -14 \\ B \$16& Black & 32 & 2 \\ &&&\\ B \$1 & Red & 0 & 1 \\ B \$2 & Black & 4 & 3 \\&&&\\ \vdots & \vdots & \vdots &\vdots \end{tabular} \end{equation*} At the end of each cycle of bets the player is another \$1 ahead. The problem is that a modest sequence of losses will exhaust their stake. The rule for time dependent covariates in a Cox model is simple and essentially the same as that for gambling: you cannot look into the future. A covariate may change in any way based on past data or outcomes, but it may not reach ``forward'' in time. One of the more well known examples of this error is analysis by response: at the end of a trial a survival curve is made comparing those who had an early response to treatment (shrinkage of tumor, lowering of cholesterol, or whatever), and it discovered that response predicts survival. The problem arises because subjects are classified as responders or non-responders from the beginning of the study, i.e., they are placed into group A or B before the response has occurred. As a consequence, any early deaths that occur before response can be assessed will be assigned to the non-responder group, even deaths that have nothing to do with the condition under study. There are many variations on the error: interpolation of the values of a laboratory test linearly between observation times, removing subjects who do not finish the treatment plan, imputing the date of an adverse event midway between observation times, etc. All of these are similar to running a red light in your car: disaster is not guarranteed --- but it is likely. The most common way to encode time-dependent covariates is to use the (start, stop] form of the model. <>= fit <- coxph(Surv(time1, time2, status) ~ age + creatinine, data=mydata) @ In this data set a patient might have the following observations \begin{center} \begin{tabular}{ccccccc} subject & time1 & time2 & status & age & creatinine & \ldots \\ \hline 1 & 0 & 15 & 0 & 25 & 1.3 \\ 1 & 15& 46 & 0 & 25 & 1.5 \\ 1 & 46& 73 & 0 & 25 & 1.4 \\ 1 & 73& 100& 1 & 25 & 1.6 \\ \end{tabular} \end{center} In this case the variable \emph{age} = age at entry to the study stays the same from line to line, while the value of creatinine varies and is treated as 1.3 over the interval $(0, 15]$, 1.5 over $(15, 46]$, etc. The intervals are open on the left and closed on the right, which means that the creatinine is taken to be 1.3 on day 15. The status variable describes whether or not each interval ended in an event. One commmon question with this data setup is whether we need to worry about correlated data, since a given subject has multiple observations. The answer is no, we do not. The reason is that this representation is simply a bookkeeping trick; the likelihood equations at any time point use only one copy of any subject, the program picks out the correct row of data at any given time. There are two exceptions to this rule, in which case the cluster variance is necessary: \begin{itemize} \item When subjects have multiple events. \item When a subject appears in overlapping intervals. This however is almost always a data error, since it corresponds to two copies of the subject being present at the same time, e.g., they could meet themselves on the sidewalk. \end{itemize} \section{Examples} \subsection{Multiple events} Chronic granulomatous disease (CGD) is a heterogenous group of uncommon inherited disorders characterized by recurrent pyogenic infections that usually begin early in life and may lead to death in childhood. Interferon gamma is a principal macrophage-activating factor shown to partially correct the metabolic defect in phagocytes. It was hypothesized that treatment with interferon might reduce the frequency of serious infections in patients with CGD. In 1986, Genentech, Inc. conducted a randomized, double-blind, placebo-controlled trial in 128 CGD patients who received Genentech's humanized interferon gamma (rIFN-g) or placebo three %' times daily for a year. The primary endpoint of the study was the time to the first serious infection. However, data were collected on all serious infections until the end of followup, which occurred before day 400 for most patients. Thirty of the 65 patients in the placebo group and 14 of the 63 patients in the rIFN-g group had at least one serious infection. The total number of infections was 56 and 20 in the placebo and treatment groups, respectively. One patient was taken off on the day of his last infection; all others have some followup after their last episode. Below are the first 10 observations, but with the listing truncatated beyond the fourth infection. Subject 2 has 7 infections, no one in the study has more. \small\begin{verbatim} 1 204 082888 1 2 12 147.0 62.0 2 2 2 2 414 219 373 2 204 082888 0 1 15 159.0 47.5 2 2 1 2 439 8 26 152 241 249 322 350 3 204 082988 1 1 19 171.0 72.7 1 2 1 2 382 4 204 091388 1 1 12 142.0 34.0 1 2 1 2 388 5 238 092888 0 1 17 162.5 52.7 1 2 1 1 383 246 253 6 245 093088 1 2 44 153.3 45.0 2 2 2 2 364 7 245 093088 0 1 22 175.0 59.7 1 2 1 2 364 292 8 245 093088 1 1 7 111.0 17.4 1 2 1 2 363 9 238 100488 0 1 27 176.0 82.8 2 2 1 1 349 294 10 238 100488 1 1 5 113.0 19.5 1 2 1 1 371 \end{verbatim} \normalsize The data set \texttt{cgd} is included in the survival package. Below we list the first few observations of that data set: <<>>= cgd[1:10, c("id", "tstart", "tstop", "status", "enum", "treat")] cfit <- coxph(Surv(tstart, tstop, status) ~ treat + sex + age + inherit + cluster(id), data=cgd) @ The logic for creating the time variables is \begin{itemize} \item If a subject has no events there is a single interval from 0 to last follow-up, with a status of 0. \item If there are events \begin{itemize} \item There is an interval for each event (0, event1], (event1, event2], etc. each with status =1. \item If the follow-up time exceeds the last event there will be a final interval with status =0. \end{itemize} \item The enum variable is simply 1,2,\ldots for each subject. It often proves useful in analysis, for instance \texttt{subset=(enum==1)} within a model restricts analysis to the first event for each subject. \end{itemize} \subsection{Changing lab tests} To be filled in. \subsection{Predictable time-dependent covariates} Occasionaly one has a time-dependent covariate whose values in the future are predictable. The most obvious of these is patient age, occasionally this may also be true for the cumulative dose of a drug. If age is entered as a linear term in the model, then the effect of changing age can be ignored in a Cox model, due to the structure of the partial likelihood. Assume that subject $i$ has an event at time $t_i$, with other subject $j \in R_i$ at risk at that time, with $a$ denoting age. The partial likelihood term is \begin{equation*} \frac{e^{\beta * a_i}}{\sum_{j \in R_i} e^{\beta* a_j}} = \frac{e^{\beta * (a_i + t_i)}}{\sum_{j \in R_i} e^{\beta* (a_j + t_i)}} \end{equation*} We see that using time-dependent age (the right hand version) or age at baseline (left hand), the partial likelihood term is identical since $\exp(\beta t_i)$ cancels out of the fraction. Howevever, if the effect of age on risk is \emph{non-linear}, this cancellation does not occur. Since age changes continuously, we would in theory need a very large (start, stop] data set to completely capture the effect --- an interval per day to match the usual resolution for death times. In practice this level of resolution is not necessary; though we all grow older, risk does not increase so rapidly that we need to know out age to the day! For most medical applications year is sufficient, but this still leads to a large data set. One useful way to generate this data set is through use of the \texttt{pyears} function. The following example uses data on rehospitalization for cohort of rheumatoid arthritis patients who also have conjestive heart failure (CHF)\cite{Nicola05}. The variables are \begin{itemize} \item patid: patient identifier \item agechf: age at onset of CHF \item yearchf calendar year of CHF, relative to the start of the study \item startday, stopday: an interval of risk \item hospevt: =1 if the interval ends with a hospitalization \item prevhosp: number of prior hosptializations \item duration: duration of RA prior to CHF \item male sex: 1 for male, 0 for female \end{itemize} The age and duration variables have been rounded to .1 year to maintain patient privacy. <<>>= load('raheart.rda') age2 <- tcut(raheart$agechf*365.25, 0:110* 365.25, labels=0:109) rowid <- 1:nrow(raheart) pfit <- pyears(Surv(startday, stopday, hospevt) ~ age2 + rowid, data=raheart, data.frame=TRUE, scale=1) print(pfit$offtable) pdata <- pfit$data print(pdata[1:6,]) @ The \texttt{tcut} function attaches a set of cutpoints to the starting age for each subject, it's primary job is to mark the variable as time-increasing%' for the \texttt{pyears} function. In the \texttt{pyears} call we set scale=1 to prevent the age intervals from being rescaled to years, this is not critical. Printing out the value of offtable is important, however. One of them most common mistakes in using \texttt{pyears} is mismatched scales, for instance if the age were in years and the follow-up time in days, and a common result of that error is to have follow-up time that fits into none of the categories. This will give a large amount of time that is outside the boundaries of the table. In the resulting data frame the first observation has 7 days of follow-up, exactly as in the starting data. The second observation has been broken into 4 rows, 109.6 days at age 92, then a year each at age 93 and 94, and a final 234.9 days at age 95 ending with a hospitalization. (The odd fractions of a day like .575 are a consequence of rounding the age values to 1 digit.) Now we combine this with the original data set using the same indexing trick found in the first example. We also need a variable containing the end time for the prior row of each subject and zero for the first row of each subject, which is \texttt{lagtime} below. We then fit two models. The first looks at the effect of age at diagnosis of CHF, the second at the effect of current age. <<>>= index <- as.integer(pdata$rowid) lagtime <- c(0, pdata$pyears[-nrow(pdata)]) lagtime[1+ which(diff(index)==0)] <- 0 #starts at 0 for each subject temp <- raheart$startday[index] + lagtime #start of each new interval data2 <- data.frame(raheart[index,], time1= temp, time2= temp + pdata$pyears, event= pdata$event, age2= 1+ as.numeric(pdata$age2) ) afit1 <- coxph(Surv(startday, stopday, hospevt) ~ male + pspline(agechf), data=raheart) afit2 <- coxph(Surv(time1, time2, event) ~ male + pspline(age2), data2) #termplot(afit1, terms=2, se=TRUE, xlab="Age at Diagnosis of CHF") #termplot(afit2, terms=2, se=TRUE, xlab="Current Age") table(with(raheart, tapply(hospevt, patid, sum))) @ In this particular case the two fits are quite similar. In retrospect this perhaps should have been expected: the mean age for onset of CHF in this group is 75 years, which does leave a lot of time for aging. (This analysis is very preliminary, however. As shown by the last line above there are a few patients with a very large number of admissions, sometimes referred to as ``entering a revolving door'' near the end of their disease. These have an untoward influence on the fit.) \subsection{Predictable covariates, method 2} Another method to create a time-changing covariate is to use the \emph{time-transform} feature of coxph. <<>>= afit2b <- coxph(Surv(startday, stopday, hospevt) ~ male + tt(agechf), data=raheart, tt=function(x, t, ...) pspline(x + t/365.25)) afit2b @ If there are one or more terms on the right hand side of the equation marked with the tt() operator, the program will pre-compute the values of that variable for each unique event time and strata combination. A user-defined function is called with arguments of \begin{itemize} \item the covariate: whatever is inside the tt() call \item the event time \item the event number: if there are multiple strata and the same event time occurs in two of them, they are treated separately \item the weight for the observation, if the call used weights \end{itemize} There is a single call to the function with a very large $x$ vector, it contains an element for each subject at risk at each event time. If there are multiple tt() terms in the formula, then the tt argument should be a list of functions with the requisite number of elements. There are other interesting uses for the time-transform capability. One example is O'Brien's logit-rank test procedure \cite{obrien78}. He proposed replacing the covariate at each event time with a logit transform of its ranks. This removes the influence of any outliers in the predictor $x$. For this case we ingore the event time argument and concentrate on the groupings. <<>>= function(x, t, riskset, weights){ obrien <- function(x) { r <- rank(x) (r-.5)/(.5+length(r)-r) } unlist(tapply(x, riskset, obrien)) } @ This relies on the fact that the input argments to tt() are ordered by the event number or riskset. This function is used as a default if no tt argument is present in the coxph call, but there are tt terms in the model formula. (Doing so allowed me to depreciate the survobrien function). Another interesting useage is to replace the data by simple ranks, not rescaled to 0--1. <<>>= function(x, t, riskset, weights) unlist(tapply(x, riskset, rank)) @ The score statistic for this model is $(C-D)/2$, where $C$ and $D$ are the number of concordant and discordant pairs, see the survConcordance function. The score statistic from this fit is then a test for significance of the concordance statistics, and is in fact the basis for the standard error reported by survConcordance. The O'Brien test can be viewed as concordance statistic that gives equal %' weight to each event time, whereas the standard concordance weights each event proportionally to the size of the risk set. (The Cox score statistic depends on the mean $x$ at each event time; since ranks go from 1 to number at risk the mean also scales.) Although handy, the computational impact of the tt argument should be considered before using it. The Cox model requires computation of a weighted mean and variance of the covariates at each event time, a process that is inherently $O(ndp^2)$ where $n$ = the sample size, $d$ = the number of events and $p$= the number of covariates. Much of the algorithmic effort in coxph() is to use updating methods for the mean and variance matrices, reducing the compute time to $O((n+d) p^2)$. When a tt term appears updating is not possible; for even moderate size data sets the impact of $nd$ versus $n+d$ can be surprising. The time-transform is a new addition and still has some rough edges. At this moment the $x=T$ argment is needed to get proper residuals and predicted values, and termplot is unable to properly reconstruct the data to plot the spline fit. Please communicate any concerns or interesting examples to the author. \begin{thebibliography}{9} \bibitem{Nicola05} Nicola PJ, Maradit-Kremers H, Roger VL, Jacobsen SJ, Crowson CS, Ballman KV, Gabriel SE. ``The risk of Congestive Heart Failure in Rheumatoid Arthritis: a Population-Based Study Over 46 Years.'' \emph{Arthritis Rheum} 52: 412--20, 2005. \bibitem{obrien78} O'Brien, Peter, ``A non-parametric test for %' association with censored data'', \emph{Biometrics} 34:243--250, 1978. \end{thebibliography} \end{document}