Time dependent coefficients in R - how to do it?

Question

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)?

Answer 1

@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."

Answer 2

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:

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:

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.

Answer 3

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.