17
$\begingroup$

Update: Sorry for another update but I've found some possible solutions with fractional polynomials and the competing risk-package that I need some help with.


The problem

I can't find an easy way to do a time dependent coefficient analysis is in R. I want to be able to take my variables coefficient and do it into a time dependent coefficient (not variable) and then plot the variation against time:

$\beta_{my\_variable}=\beta_0+\beta_1*t+\beta_2*t^2...$

Possible solutions

1) Splitting the dataset

I've looked at this example (Se part 2 of the lab session) but the creation of a separate dataset seems complicated, computationally costly and not very intuitive...

2) Reduced Rank models - The coxvc package

The coxvc package provides an elegant way of dealing with the problem - here's a manual. The problem is that the author is no longer developing the package (last version is since 05/23/2007), after some e-mail conversation I've gotten the package to work but one run took 5 hours on my dataset (140 000 entries) and gives extreme estimates at the end of the period. You can find a slightly updated package here - I've mostly just updated the plot function.

It might be just a question of tweaking but since the software doesn't easily provide confidence intervals and the process is so time consuming I'm looking right now at other solutions.

3) The timereg package

The impressive timereg package also addresses the problem but I'm not certain of how to use it and it doesn't give me a smooth plot.

4) Fractional Polynomial Time (FPT) model

I found Anika Buchholz' excellent dissertation on "Assessment of time–varying long–term effects of therapies and prognostic factors" that does an excellent job covering different models. She concludes that Sauerbrei et al's proposed FPT seems to be the most appropriate for time-dependent coefficients:

FPT is very good at detecting time-varying effects, while the Reduced Rank approach results in far too complex models, as it does not include selection of time-varying effects.

The research seems very complete but it's slightly out of reach for me. I'm also a little wondering since she happens to work with Sauerbrei. It seems sound though and I guess the analysis could be done with the mfp package but I'm not sure how.

5) The cmprsk package

I've been thinking of doing my competing risk analysis but the calculations have been to time-consuming so I switched to the regular cox regression. The crr has thoug an option for time dependent covariates:

....
cov2        matrix of covariates that will be multiplied 
            by functions of time; if used, often these 
            covariates would also appear in cov1 to give 
            a prop hazards effect plus a time interaction
....

There is the quadratic example but I'm don't quite follow where the time actually appears and I'm not sure of how to display it. I've also looked at the test.R file but the example there is basically the same...

My example code

Here's an example that I use to test the different possibilities

library("survival")
library("timereg")
data(sTRACE)

# Basic cox regression    
surv <- with(sTRACE, Surv(time/365,status==9))
fit1 <- coxph(surv~age+sex+diabetes+chf+vf, data=sTRACE)
check <- cox.zph(fit1)
print(check)
plot(check, resid=F)
# vf seems to be the most time varying

######################################
# Do the analysis with the code from #
# the example that I've found        #
######################################

# Split the dataset according to the splitSurv() from prof. Wesley O. Johnson
# http://anson.ucdavis.edu/~johnson/st222/lab8/splitSurv.ssc
new_split_dataset = splitSuv(sTRACE$time/365, sTRACE$status==9, sTRACE[, grep("(age|sex|diabetes|chf|vf)", names(sTRACE))])

surv2 <- with(new_split_dataset, Surv(start, stop, event))
fit2 <- coxph(surv2~age+sex+diabetes+chf+I(pspline(stop)*vf), data=new_split_dataset)
print(fit2)

######################################
# Do the analysis by just straifying #
######################################
fit3 <- coxph(surv~age+sex+diabetes+chf+strata(vf), data=sTRACE)
print(fit3)

# High computational cost!
# The price for 259 events
sum((sTRACE$status==9)*1)
# ~240 times larger dataset!
NROW(new_split_dataset)/NROW(sTRACE)

########################################
# Do the analysis with the coxvc and   #
# the timecox from the timereg library #
########################################
Ft_1 <- cbind(rep(1,nrow(sTRACE)),bs(sTRACE$time/365,df=3))
fit_coxvc1 <- coxvc(surv~vf+sex, Ft_1, rank=2, data=sTRACE)

fit_coxvc2 <- coxvc(surv~vf+sex, Ft_1, rank=1, data=sTRACE)

Ft_3 <- cbind(rep(1,nrow(sTRACE)),bs(sTRACE$time/365,df=5))
fit_coxvc3 <- coxvc(surv~vf+sex, Ft_3, rank=2, data=sTRACE)

layout(matrix(1:3, ncol=1))
my_plotcoxvc <- function(fit, fun="effects"){
    plotcoxvc(fit,fun=fun,xlab='time in years', ylim=c(-1,1), legend_x=.010)
    abline(0,0, lty=2, col=rgb(.5,.5,.5,.5))
    title(paste("B-spline =", NCOL(fit$Ftime)-1, "df and rank =", fit$rank))
}
my_plotcoxvc(fit_coxvc1)
my_plotcoxvc(fit_coxvc2)
my_plotcoxvc(fit_coxvc3)

# Next group
my_plotcoxvc(fit_coxvc1)

fit_timecox1<-timecox(surv~sex + vf, data=sTRACE)
plot(fit_timecox1, xlab="time in years", specific.comps=c(2,3))

The code results in these graphs: Comparison of different settings for coxvc and of the coxvc and the timecox plots. I guess the results are ok but I don't think I'll be able to explain the timecox graph - it seems to complex...

My (current) questions

  • How do I do the FPT analysis in R?
  • How do I use the time covariate in cmprsk?
  • How do I plot the result (preferably with confidence intervals)?
$\endgroup$
  • 3
    $\begingroup$ The example in the link is about time-varying covariates, not time varying coefficients. These are too different things. To get time varying parameters the way you described use interactions, i.e. instead of model $y=x\beta_{my}$, fit model $y=x\beta_0+x\cdot t\beta_1+x\cdot t^2\beta_2$. In R the formulas would look like y~x and y~x*(t+t^2)-t (also possible to write y~x+x:t+x:t^2 for the last one). $\endgroup$ – mpiktas Nov 18 '11 at 14:45
  • $\begingroup$ I thought the second part: "2. Model Deterministic Time-dependent Covariates to Check the PH Assumption" would be the part dealing with my question. I was hoping to do something of the formula that you describe but when I tried it I either got an error or a separate time variable... I got a low p-value for time though :-D $\endgroup$ – Max Gordon Nov 18 '11 at 21:08
  • $\begingroup$ @max-gordon, is your response variable a quantity, or the time elapsed until an even occurs? Because most of the methods you cite are specifically for time-to-event data. $\endgroup$ – f1r3br4nd Jul 10 '13 at 20:02
  • $\begingroup$ @f1r3br4nd: It is a quantity (Age in my study) where the hazard is non-proportional, i.e. it varies over time in my time-to-event model. In the end I decided to split into two different time-frames as I was not thrilled by doing a 3-D graph - that would have never passed the reviewers... $\endgroup$ – Max Gordon Jul 11 '13 at 20:12
  • $\begingroup$ There's a difference between time dependent/varying predictors and time interaction. Most variables are time dependent (sex is an exception). If you have one observation per person, then you'll have little or no chance to perform a time dependent/varying analysis. Anderson-Gill's method is the most frequently used for time- dependent survival analysis. The advantage of time dependent methods is that values during follow-up might be more predicitve of survival experience than baseline values. The second situation, time interacting predictors are simply tests of the PH assumption. $\endgroup$ – Adam Robinsson Aug 15 '15 at 6:31
8
+50
$\begingroup$

@mpiktas came close in offering a feasible model, however the term that needs to be used for the quadratic in time=t would be I(t^2)) . This is so because in R the formula interpretation of "^" creates interactions and does not perform exponentiation, so the interaction of "t" with "t" is just "t". (Shouldn't this be migrated to SO with an [r] tag?)

For alternatives to this process, which looks to me somewhat dubious for inference purposes, and one which probably fits your interest in using the supportive functions in Harrell's rms/Hmisc packages, see Harrell's "Regression Modeling Strategies". He mentions (but only in passing although he does cite some of his own papers) constructing spline fits to the time scale to model bathtub-shaped hazards. His chapter on parametric survival models describes a variety of plotting techniques for checking proportional hazards assumptions and for examining the linearity of estimated log-hazard effects on the time scale.

Edit: An additional option is to use coxph's 'tt' parameter described as an "optional list of time-transform functions."

$\endgroup$
  • $\begingroup$ I agree that this should probably be moved to the SO [r] tag. $\endgroup$ – Zach Dec 11 '11 at 19:14
  • $\begingroup$ +1 for for your answer, I wasn't aware that this would be so hard answering. It seems like a common problem, perhaps the question is more of a question of coding than and you might be right about SO being a better choice. I tried your formula it seems that vf+I(vflog(time)) has an excellent fit, I tried just vftime and vf*time^2 but the log gave by fare the lowest p-value. I tried to run it with the cph() function to get the AIC but it gave an error :( Do you have any idea of how to do a plot on the estimate? $\endgroup$ – Max Gordon Dec 11 '11 at 21:10
  • $\begingroup$ I thought that check <- cox.zph(fit1); print(check); plot(check, resid=F) as in your set up gave informative plots of the time "effect". Did you mean cph() which is from the rms package or coxph from survival? $\endgroup$ – DWin Dec 11 '11 at 21:53
  • $\begingroup$ Yes, the Schoenfeld residuals does give a nice idea of the time variation but I think people might have a hard time understanding it. The plot gives as I understand it the residual variation not explained by my model. I would like though a plot where I have the complete variable effect on the y-axis and time on the x-axis, I believe that this would be easier to interpret since you don't have to look at both the table and the plot to get the hazard at a specific point in time... Yes I meant cph() and not the coxph() since that one doesn't work with the AIC() $\endgroup$ – Max Gordon Dec 12 '11 at 5:37
  • $\begingroup$ I'm also a little confused to why I've found all the complex methods described in my question while this I(variable*time) seems very straight forward and intuitive - as a non-statistician I'm thinking - what have I missed? $\endgroup$ – Max Gordon Dec 12 '11 at 5:41
5
$\begingroup$

I've changed the answer to this as neither @DWin's or @Zach's answers fully answers how to model time-varying coefficients. I've recently wrote a post about this. Here's the gist of it.

The core concept in the Cox regression model is the hazard function, $h(t)$. It is defined as:

$$h(t) = \frac{f(t)}{S(t)}$$

Where the $f(t)$ is the risk of having an event at any given time while the $S(t)$ is the probability of surviving that fare. The number is thus a fraction with a theoretical range from $0$ to $\infty$.

An interesting feature of the hazard function is that we can include observations at other points in time than $time_0$, e.g. if "Peter" is operated with a hip arthroplasty in England, arrives after 1 year to Sweden, he has been alive for 1 year when we choose to include him. Note that any patient that would have died prior to that would never have come to our attention and we can therefore not include Peter in the $S(t)$ when looking at the hazard prior to 1 year. After 1 year we may include Peter.

When allowing subjects entering at other time points we must change the Surv from Surv(time, status) to Surv(start_time, end_time, status). While the end_time is strongly correlated with the outcome the start_time is now available as an interaction term (just as hinted in the original suggestions). In a regular setting the start_time is 0 except for a few subjects that appear later but is we split each observation into several periods we suddenly have plenty start times that are non-zero. The only difference is that each observation occurs multiple times where all but the last observation has the option of a non-censored outcome.

Time splitting in practice

I've just published on CRAN the Greg package that makes this time-split easy. First we start with some theoretical observations:

library(Greg)
test_data <- data.frame(
  id = 1:4,
  time = c(4, 3.5, 1, 5),
  event = c("censored", "dead", "alive", "dead"),
  age = c(62.2, 55.3, 73.7, 46.3),
  date = as.Date(
    c("2003-01-01", 
      "2010-04-01", 
      "2013-09-20",
      "2002-02-23"))
)

We can show this graphically with * being an indicator of event:

enter image description here

If we apply the timeSplitter as following:

library(dplyr)
split_data <- 
  test_data %>% 
  select(id, event, time, age, date) %>% 
  timeSplitter(by = 2, # The time that we want to split by
               event_var = "event",
               time_var = "time",
               event_start_status = "alive",
               time_related_vars = c("age", "date"))

We get the following:

enter image description here

As you can see each object has been split into multiple events where the last time span contains the actual event status. This allows us to now build models that have simple : interaction terms (do not use the * as this expands to time + var + time:var and we're not interested in the time per se). There is no need for using the I() function although if you want to check nonlinearity with time I often create a separate time-interaction variable that I add a spline to and then display using rms::contrast. Anyway, your regression call should now look like this:

coxp(Surv(start_time, end_time, event) ~ var1 + var2 + var2:time, 
     data = time_split_data)

Using the survival package's tt function

There is also a way to model time dependent coefficients directly in the survival package using the tt function. Prof. Therneau provides a thorough introduction to the subject in his vignette. Unfortunately in large datasets this is hard to do due to memory limitations. It seems that the tt function splits the time into very fine pieces generating in the process a huge matrix.

$\endgroup$
2
$\begingroup$

You can use the apply.rolling function in PerformanceAnalytics to run a linear regression through a rolling window, which will allow your coefficients to vary over time.

For example:

library(PerformanceAnalytics)
library(quantmod)
getSymbols(c('AAPL','SPY'), from='01-01-1900')
chart.RollingRegression(Cl(AAPL),Cl(SPY), width=252, attribute='Beta')
#Note: Alpha=y-intercept, Beta=regression coeffient

This works with other functions too.

$\endgroup$
  • $\begingroup$ Thank you for your answer, I guess a moving time-window should work just as well as my approaches. I can't get your example to run though, could you please give an example based on my sTRACE example so that I know exactly how to implement it? $\endgroup$ – Max Gordon Dec 10 '11 at 13:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.