Discussion of Time Dependent Covariate Data Structure for "coxph" in R
======================================================================
E. Slud,  12/9/2021

This material is adapted and from the TimDep.RLog handout, topic (8) in the 
Statistical Computing Handouts.

> library(survival)
> data(cancer)   ## "NCCTG lung cancer study" = Mayo lung cancer study
  data(lung)     ## same dataset, 228 x 10
> names(lung)
 [1] "inst"      "time"      "status"    "age"       "sex"       "ph.ecog"  
 [7] "ph.karno"  "pat.karno" "meal.cal"  "wt.loss"  

  lung$status = lung$status-1    ### was originally coded 1,2

## Objective here is to take a dataset with only baseline covariates and
#   code it into a (start, stop) format to handle the covariate 
> sextim = (2-cancer$sex)*log(cancer$time)

## The function "tmerge" in the "survival" package can do this. We do it by hand 
##   below, for the particular case of deterministically time-varying covariate.

> ordtim = sort(unique(lung$time))             ## 186 unique times
  maxind = c(outer(lung$time, ordtim, function(x,y)  x>= y) %*% rep(1,186)) 
   ## these are the indices of the 228 patients' times among the ordered times

> sum(abs(ordtim[maxind]-lung$time))     ## this is the verification
[1] 0

  repind = rep(1:228,maxind) ## repeat each index i maxind[i] times
  length(repind)
[1] 20022       

## build up the vectors "start" and "stop" times step by step, in a for-loop

> start = NULL; stop = NULL
  for(i in 1:nrow(lung))    {
      stop = c(stop, ordtim[1:maxind[i]])
      start = c(start, c(0,ordtim)[1:maxind[i]])   }
  newfram = cbind.data.frame(start=start, stop=stop, status = 0,
      lung[repind, c(5:7,10)], sextim = 
      (2-lung$sex[repind])*log(stop))

> newfram$status[cumsum(maxind)] = lung$status
    0     1 
19857   165             ### agrees with 165=sum(lung$status)

### Now have a 10022 x 8 data-frame, with the following columns
> names(newfram)
[1] "start"    "stop"     "status"   "sex"      "ph.ecog"  "ph.karno" "wt.loss" 
[8] "sextim"
 
## Finally, we can do the time-dependent Cox-model fit, which we compare with 
#   the same Cox-model fit wuthout the time-dependent covariate

> lungA = coxph(Surv(time,status) ~ sex + ph.ecog +
      ph.karno + wt.loss, data=lung)
  lungA
Call:
coxph(formula = Surv(time, status) ~ sex + ph.ecog + ph.karno + 
    wt.loss, data = lung)

              coef exp(coef)  se(coef)      z        p
sex      -0.623193  0.536229  0.176389 -3.533 0.000411
ph.ecog   0.753428  2.124270  0.194168  3.880 0.000104
ph.karno  0.013814  1.013910  0.009908  1.394 0.163241
wt.loss  -0.009008  0.991032  0.006511 -1.383 0.166516

Likelihood ratio test=31.06  on 4 df, p=2.97e-06
n= 213, number of events= 151 
   (15 observations deleted due to missingness)

> lungB = coxph(Surv(start,stop,status) ~sex + ph.ecog +
      ph.karno + wt.loss, data=newfram)       
         ## uses time-dependent framework but has no time-dependent covariates
  lungB
Call:
coxph(formula = Surv(start, stop, status) ~ sex + ph.ecog + ph.karno + 
    wt.loss, data = newfram)

              coef exp(coef)  se(coef)      z        p
sex      -0.623193  0.536229  0.176389 -3.533 0.000411
ph.ecog   0.753428  2.124270  0.194168  3.880 0.000104
ph.karno  0.013814  1.013910  0.009908  1.394 0.163241
wt.loss  -0.009008  0.991032  0.006511 -1.383 0.166516

Likelihood ratio test=31.06  on 4 df, p=2.97e-06
n= 19303, number of events= 151 
   (719 observations deleted due to missingness)

> lungC = coxph(Surv(start,stop,status) ~sex + ph.ecog +
      ph.karno + wt.loss + sextim, data=newfram)   
  lungC
Call:
coxph(formula = Surv(start, stop, status) ~ sex + ph.ecog + ph.karno + 
    wt.loss + sextim, data = newfram)

              coef exp(coef)  se(coef)      z        p
sex      -2.006232  0.134495  1.127879 -1.779 0.075278
ph.ecog   0.735771  2.087091  0.194670  3.780 0.000157
ph.karno  0.012474  1.012553  0.009965  1.252 0.210615
wt.loss  -0.009612  0.990434  0.006577 -1.462 0.143877
sextim   -0.255592  0.774458  0.204005 -1.253 0.210254

Likelihood ratio test=32.76  on 5 df, p=4.204e-06
n= 19303, number of events= 151 
   (719 observations deleted due to missingness)

  ### To account for the patients with some missing data:

> sapply(lung[,c(5:7,10)], function(col) sum(is.na(col)))
     sex  ph.ecog ph.karno  wt.loss 
       0        1        1       14