search

harrelfe / rms

Regression Modeling Strategies
https://hbiostat.org/R/rms
Other
177 stars 49 forks source link

Time-dependent coefficients for cph #6

Open gforge opened 10 years ago

gforge commented 10 years ago

Dear prof. Harrell,

I've recently learned of a way to work with time-dependent coefficients when the proportional hazards assumption is violated in the cox-regression: By splitting the dataset using the Epi-package I use the new start-times as an interaction term. Unfortunately this use was not anticipated and causes some issues with the cph() and its related functions.

I've forked your package and worked on a fix but while I get the same coefficients as coxph() gives (although there are some minimal differences) my favorite function, the contrast(), still doesn't work with this fix. I'm not entirely sure how to best implement the needed functionality. You can find my test-file here and the fixes are on line 143 and line 152-155 .

harrelfe commented 10 years ago

Please don't fork the package until I get a chance to correct simple errors.

I tried running your script but received an error before getting to time-dependent covariates with cph:

lxs_obj <- Lexis(entry = list(Timeband = Start),

  • exit = list(Timeband = dt, Age = age + dt),
  • exit.status = e,
  • data = test_df) Error in Lexis(entry = list(Timeband = Start), exit = list(Timeband = dt, : Duration is not the same on all time scales

And best not to add "obj" or "df" to the end of object names and not to include underscores. These slow down typing for me and "obj" and "df" don't provide new information.

Frank

On 08/24/2014 07:56 AM, Max Gordon wrote:

Dear prof. Harrell,

I've recently learned of a way to work with time-dependent coefficients when the proportional hazards assumption is violated in the cox-regression: By splitting the dataset using the Epi-package I use the new start-times as an interaction term. Unfortunately this use was not anticipated and causes some issues with the cph() and its related functions.

I've forked your package and worked on a fix but while I get the same coefficients as coxph() gives (although there are some minimal differences) my favorite function, the contrast(), still doesn't work with this fix. I'm not entirely sure how to best implement the needed functionality. You can find my test-file here https://github.com/gforge/rms/blob/master/tests/testthat/test-time_interaction.R and the fixes are in this commit https://github.com/gforge/rms/commit/87b9a367c2aca98f79a432541df75b367de520d8.

— Reply to this email directly or view it on GitHub https://github.com/harrelfe/rms/issues/6.

Frank E Harrell Jr Professor and Chairman School of Medicine

Department of *Biostatistics*   *Vanderbilt University*
gforge commented 10 years ago

Sorry, the label() was causing errors. I thought it was the units() but forgot to properly check after removing them. I've also cleaned the variable names.

The fix uses the assume.code == 1 as I assumed that it indicates an asis-object but I think that the problem may actually originate from the assume.code not being 9 (the interaction indicator) but I must admit that I don't really grasp this section of the code.

I regularly retrieve and incorporate you commits into the fork, especially before I try to fix anything. I think it's more convenient to investigate and report issues using the very latest Github-version. It should limit the risk of fixing the same bug twice. I keep my test-code within the testthat folder in order to separate my test-cases from the originals (I frequently use the testthat-package, after some debugging I found that with the <<- it plays rather nicely with the datadist-setup).

harrelfe commented 10 years ago

rms was not designed to work with less-than-full-rank models so don't change anything related to how the different terms are recognized. This will likely break something else.

Frank

On 08/24/2014 09:16 AM, Max Gordon wrote:

Sorry, the label() was causing errors. I thought it was the units() but forgot to properly check after removing them. I've also cleaned the variable names.

The fix uses the assume.code == 1 as I assumed that it indicates an asis-object but I think that the problem may actually originate from the assume.code not being 9 (the interaction indicator) but I must admit that I don't really grasp this section of the code.

I regularly retrieve and incorporate you commits into the fork, especially before I try to fix anything. I think it's more convenient to investigate and report issues using the very latest Github-version. It should limit the risk of fixing the same bug twice. I keep my test-code within the testthat folder in order to separate my test-cases from the originals (I frequently use the testthat-package, after some debugging I found that with the <<- it plays rather nicely with the datadist-setup).

— Reply to this email directly or view it on GitHub https://github.com/harrelfe/rms/issues/6#issuecomment-53194655.

Frank E Harrell Jr Professor and Chairman School of Medicine

Department of *Biostatistics*   *Vanderbilt University*
gforge commented 10 years ago

Ok, thank you. I can get around the problem by manually generating an interaction variable - seems to work satisfactory:

spl_alt <- 
  within(spl, {
    Male_time_int = (sex == "Male")*Timeband
  })

spl_alt$lex.Cst <- NULL
spl_alt$Start <- NULL
dd <- datadist(spl_alt)
options(datadist = "dd")

model <-
  cph(Surv(time = Timeband, time2 = Stop, event = lex.Xst) ~ 
        Age + sex + Male_time_int,
      data = spl_alt)

contrast(model, 
         a = list(sex = "Male", Male_time_int = 0:5),
         b = list(sex = "Female", Male_time_int = 0))
harrelfe commented 10 years ago

It's too bad that rms needs you to do it that way but that is probably the most expedient solution.

On 08/27/2014 10:32 AM, Max Gordon wrote:

Ok, thank you. I can get around the problem by manually generating an interaction variable - seems to work satisfactory:

spl_alt<-
within(spl, { Male_time_int= (sex== "Male")*Timeband })

spl_alt$lex.Cst<- NULL spl_alt$Start<- NULL dd<- datadist(spl_alt) options(datadist= "dd")

model<- cph(Surv(time= Timeband, time2= Stop, event= lex.Xst) ~
Age+ sex+ Male_time_int, data= spl_alt)

contrast(model,
a= list(sex= "Male", Male_time_int= 0:5), b= list(sex= "Female", Male_time_int= 0))

— Reply to this email directly or view it on GitHub https://github.com/harrelfe/rms/issues/6#issuecomment-53590401.

Frank E Harrell Jr Professor and Chairman School of Medicine

Department of *Biostatistics*   *Vanderbilt University*