<p>This is essentially a 3 part question: </p>

<ol>
<li>How to estimate time-varying effects?</li>
<li>What is the difference between different specifications of time-varying effects using <code>survival::coxph</code> function</li>
<li>How to decide what shape the time-variation has, i.e., linear, logarithmic, ...</li>
</ol>

<p>I will try to answer these questions in the following using the veteran data example, which is featured in section 4.2 of the <a href="https://cran.microsoft.com/web/packages/survival/vignettes/timedep.pdf" rel="noreferrer">vignette on time-dependent covariates and time-dependent coefficients</a> (also known as time-varying effects) in the <strong><code>survival</code></strong> package: </p>

<pre><code>library(dplyr)
library(survival)

data("veteran", package = "survival")
veteran &lt;- veteran %&gt;%
  mutate(
    trt   = 1L * (trt == 2),
    prior = 1L * (prior == 10))
head(veteran)
#&gt;   trt celltype time status karno diagtime age prior
#&gt; 1   0 squamous   72      1    60        7  69     0
#&gt; 2   0 squamous  411      1    70        5  64     1
#&gt; 3   0 squamous  228      1    60        3  38     0
#&gt; 4   0 squamous  126      1    60        9  63     1
#&gt; 5   0 squamous  118      1    70       11  65     1
#&gt; 6   0 squamous   10      1    20        5  49     0
</code></pre>

<h2>1. How to estimate time-varying effects</h2>

<p>There are different popular methods and implementations, e.g. <code>survival::coxph</code>, <code>timereg::aalen</code> or using GAMs after appropriate data transformation (see below).</p>

<p>Although the specific methods and their implementaitons differ, a general idea ist to create a long form data set where </p>

<ul>
<li>the follow-up is partitioned into intervals</li>
<li>for each subject, the status is 0 in all intervals except the last (if an event) </li>
<li>the time variable is updated in each interval</li>
</ul>

<p>Then, the time (or a transformation of time, e.g. log(t)) is simply a covariate and time-varying effects can be estimated by an interaction between the covariate of interest and the (transformed) covariate of time. </p>

<p>If the functional form of the time-variation is known, you can use the <code>tt()</code> aproach: </p>

<pre><code>cph_tt &lt;- coxph(
      formula = Surv(time, status) ~ trt + prior + karno + tt(karno),
      data    = veteran,
      tt      = function(x, t, ...) x * log(t + 20))
</code></pre>

<h2>2. What is the difference between different specifications of time-varying effects using <code>survival::coxph</code> function</h2>

<p>There is no difference. I assume the <code>tt()</code> function is simply a short-cut for the estimation via transformation to the long-format. You can verify that the two approaches are equivalent using the code below: </p>

<h2>transform to long format</h2>

<pre><code>veteran_long &lt;- survSplit(Surv(time, status)~., data = veteran, id = "id",
  cut = unique(veteran$time)) %&gt;%
  mutate(log_time = log(time + 20))
head(veteran_long) %&gt;% select(id, trt, age, tstart, time, log_time, status)
#&gt;   id trt age tstart time log_time status
#&gt; 1  1   0  69      0    1 3.044522      0
#&gt; 2  1   0  69      1    2 3.091042      0
#&gt; 3  1   0  69      2    3 3.135494      0
#&gt; 4  1   0  69      3    4 3.178054      0
#&gt; 5  1   0  69      4    7 3.295837      0
#&gt; 6  1   0  69      7    8 3.332205      0

cph_long &lt;- coxph(formula = Surv(tstart, time, status)~
    trt + prior + karno + karno:log_time, data = veteran_long)

## models are equivalent, just different specification
cbind(coef(cph_long), coef(cph_tt))
#&gt;                       [,1]        [,2]
#&gt; trt             0.01647766  0.01647766
#&gt; prior          -0.09317362 -0.09317362
#&gt; karno          -0.12466229 -0.12466229
#&gt; karno:log_time  0.02130957  0.02130957
</code></pre>

<h2>3. How to decide what shape the time-variation has?</h2>

<p>As mentioned before, time-varying effects are simply interactions of a covariate <code>x</code> and time <code>t</code>, thus time-varying effects can have different specifications, equivalent to interactions in standard regression models, e.g.</p>

<ul>
<li><code>x*t</code>:  linear covariate effect, linearly time-varying effect </li>
<li><code>f(x)*t</code>: non-linear covariate effect, linearly time-varying effect</li>
<li><code>f(t)*x</code>: linear covariate effect, non-linearly time-varying (for categorical x) this essentially represents a stratified baseline hazard </li>
<li><code>f(x, t)</code>: non-linear, non-linearly time-varying effect</li>
</ul>

<p>In each case, the functional form of the effect <code>f</code> can either be estimated from the data or prespecified (e.g. <code>f(t)*x = karno * log(t + 20)</code> above). </p>

<p>In most cases you would prefer to estimate <code>f</code> from the data. The support for the (penalized) estimation of such effects is to my knowledge limited in the <strong><code>survival</code></strong> package. However, you can use <code>mgcv::gam</code> to estimate any of the effects specified above (after appropriate data transformation). An example is given below and shows that the effect of <code>karno</code> goes towards 0 as time progresses, regardless of the Karnofsky score at the beginning of the follow-up (see <a href="https://adibender.github.io/pammtools/articles/tveffects.html" rel="noreferrer">here</a> for details and also Section 4.2 <a href="https://arxiv.org/pdf/1806.01042.pdf" rel="noreferrer">here</a>): </p>

<pre><code>library(pammtools)
# data transformation
ped &lt;- as_ped(veteran, Surv(time, status)~., max_time = 400)
# model
pam &lt;- mgcv::gam(ped_status ~ s(tend) + trt + prior + te(tend, karno, k = 10),
  data = ped, family = poisson(), offset = offset, method = "REML")
p_2d &lt;- gg_tensor(pam)
p_slice &lt;- gg_slice(ped, pam, "karno", tend = unique(tend), karno = c(20, 50, 80), reference = list(karno = 60))
gridExtra::grid.arrange(p_2d, p_slice, nrow = 1)
</code></pre>

<p><a href="https://i.stack.imgur.com/llPOl.png" rel="noreferrer"><img src="https://i.stack.imgur.com/llPOl.png" alt=""></a></p>
How to model interaction of covariate with time when proportionality assumption is violated in survival analysis in R
<p>In R, what is the best way to incorporate the interaction term between a covariate and time, when the proportionality test (with coxph) shows that the proportionality assumption in the Cox model is violated? I know that you can either use strata or an interaction with time term, I'm interested in the latter. I haven't been able to find a definitive clear explanation with examples on how to do this on the internet. In the most common example using the Rossi dataset, Fox suggested to do,</p>

<pre><code>coxph(formula = Surv(start, stop, arrest.time) ~ fin + age + age:stop + prio, data = Rossi.2)
</code></pre>

<p>Is there a difference between modeling with age:stop versus age:start? Does the formula have to use this format? If I use the Surv with the two parameter format, would the following also make sense?</p>

<pre><code>coxph(formula = Surv(week, arrest) ~ fin + age + age:week + prio, data = Rossi)
</code></pre>

<p>Or you have to split the dataset and use the Surv(start,stop,event) method?
 Also, there is the time-transform method, so,</p>

<pre><code>coxph(formula = Surv(week, arrest) ~ fin + age + tt(age) + prio, data = Rossi, tt=function(x,t,...) x*t)
</code></pre>

<p>I know that some people would prefer model with the <code>log(t)</code> instead of <code>t</code> here. But which one of these is the correct method to model interaction with time? Do these all refer to the same/different underlying statistical model? And the end, are all modeling (for the interaction term): <code>h(t) = h0(t)exp(b*X*t)</code>?</p>