I&#39;m trying to build a survival model in JAGS that allows for time-varying covariates. I&#39;d like it to be a parametric model — for example, assuming survival follows the Weibull distribution (but I&#39;d like to allow the hazard to vary, so exponential is too simple). So, this is essentially a Bayesian version of what can be done in the `flexsurv` package, which allows for time-varying covariates in parametric models.

Therefore, I want to be able to enter the data in a &#39;counting-process&#39; form, where each subject has multiple rows, each corresponding to a time interval in which their covariates remained constant (as described [in this pdf](https://cran.r-project.org/web/packages/survival/vignettes/timedep.pdf) or [here](http://pages.uoregon.edu/aarong/teaching/G4075_Outline/node31.html). This is the `(start, stop]` formulation that the `survival` or `flexurv` packages allow.

Unfortunately, every explanation of how to perform survival analysis in JAGS seems to assume one row per subject.

I attempted to take this simpler approach and extend it to the counting process format, but the model does not correctly estimate the distribution.

A Failed Attempt:

Here&#39;s an example. First we generate some data:

    library(&#39;dplyr&#39;)
    library(&#39;survival&#39;)
    
    ## Make the Data: -----
    set.seed(3)
    n_sub &lt;- 1000
    current_date &lt;- 365*2
    
    true_shape &lt;- 2
    true_scale &lt;- 365
    
    dat &lt;- data_frame(person = 1:n_sub,
                      true_duration = rweibull(n = n_sub, shape = true_shape, scale = true_scale),
                      person_start_time = runif(n_sub, min= 0, max= true_scale*2),
                      person_censored = (person_start_time + true_duration) &gt; current_date,
                      person_duration = ifelse(person_censored, current_date - person_start_time, true_duration)
    )

      person person_start_time person_censored person_duration
       (int)             (dbl)           (lgl)           (dbl)
    1      1          11.81416           FALSE        487.4553
    2      2         114.20900           FALSE        168.7674
    3      3          75.34220           FALSE        356.6298
    4      4         339.98225           FALSE        385.5119
    5      5         389.23357           FALSE        259.9791
    6      6         253.71067           FALSE        259.0032
    7      7         419.52305            TRUE        310.4770

Then we split the data into 2 observations per subject. I&#39;m just splitting each subject at time = 300 (unless they didn&#39;t make it to time=300, in which they get just one observation).
        
    ## Split into multiple observations per person: --------
    cens_point &lt;- 300 # &lt;----- try changing to 0 for no split; if so, model correctly estimates
    dat_split &lt;- dat %&gt;%
      group_by(person) %&gt;%
      do(data_frame(
        split = ifelse(.$person_duration &gt; cens_point, cens_point, .$person_duration),
        START = c(0, split[1]),
        END = c(split[1], .$person_duration),
        TINTERVAL = c(split[1], .$person_duration - split[1]), 
        CENS = c(ifelse(.$person_duration &gt; cens_point, 1, .$person_censored), .$person_censored), # &lt;— edited original post here due to bug; but problem still present when fixing bug
        TINTERVAL_CENS = ifelse(CENS, NA, TINTERVAL),
        END_CENS = ifelse(CENS, NA, END)
      )) %&gt;%
      filter(TINTERVAL != 0)

      person    split START      END TINTERVAL CENS TINTERVAL_CENS
       (int)    (dbl) (dbl)    (dbl)     (dbl) (dbl)        (dbl)
    1      1 300.0000     0 300.0000 300.00000     1           NA
    2      1 300.0000   300 487.4553 187.45530     0    187.45530
    3      2 168.7674     0 168.7674 168.76738     1           NA
    4      3 300.0000     0 300.0000 300.00000     1           NA
    5      3 300.0000   300 356.6298  56.62979     0     56.62979
    6      4 300.0000     0 300.0000 300.00000     1           NA

Now we can set up the JAGS model.
    
    ## Set-Up JAGS Model -------
    dat_jags &lt;- as.list(dat_split)
    dat_jags$N &lt;- length(dat_jags$TINTERVAL)
    inits &lt;- replicate(n = 2, simplify = FALSE, expr = {
           list(TINTERVAL_CENS = with(dat_jags, ifelse(CENS, TINTERVAL + 1, NA)),
                END_CENS = with(dat_jags, ifelse(CENS, END + 1, NA)) )
    })
    
    model_string &lt;- 
    &quot;
      model {
        # set priors on reparameterized version, as suggested
        # here: https://sourceforge.net/p/mcmc-jags/discussion/610036/thread/d5249e71/?limit=25#8c3b
        log_a ~ dnorm(0, .001) 
        log(a) &lt;- log_a
        log_b ~ dnorm(0, .001)
        log(b) &lt;- log_b
        nu &lt;- a
        lambda &lt;- (1/b)^a
        
        for (i in 1:N) {
          # Estimate Subject-Durations:
          CENS[i] ~ dinterval(TINTERVAL_CENS[i], TINTERVAL[i])
          TINTERVAL_CENS[i] ~ dweibull( nu, lambda )
        }
      }
    &quot;
    
    library(&#39;runjags&#39;)
    param_monitors &lt;- c(&#39;a&#39;, &#39;b&#39;, &#39;nu&#39;, &#39;lambda&#39;) 
    fit_jags &lt;- run.jags(model = model_string,
                         burnin = 1000, sample = 1000, 
                         monitor = param_monitors,
                         n.chains = 2, data = dat_jags, inits = inits)
    # estimates:
    fit_jags
    # actual:
    c(a=true_shape, b=true_scale)

Depending on where the split point is, the model estimates very different parameters for the underlying distribution. It only gets the parameters right if the data isn&#39;t split into the counting process form. It seems like this is not the way to format the data for this kind of problem. 

If I am missing an assumption and my problem is less related to JAGS and more related to how I&#39;m formulating the problem, suggestions are very welcome. I might be despairing that time-varying covariates can&#39;t be used in parametric survival models (and can only be used in models like the Cox model, which assumes constant hazards and which doesn&#39;t actually estimate the underlying distribution)— however, as I mentioned above, the `flexsurvreg` package in R does accommodate the `(start, stop]` formulation in parametric models.

If anyone knows how to build a model like this in another language (e.g. STAN instead of JAGS) that would be appreciated too.

### Edit:

Chris Jackson provides some helpful advice via email:

&gt; I think the T() construct for truncation in JAGS is needed here.   Essentially for each period (t[i], t[i+1]) where a person is alive but the covariate is constant, the survival time is left-truncated at the start of the period, and possibly also right-censored at the end. So you&#39;d write something like `y[i] ~ dweib(shape, scale[i])T(t[i], )`

I tried implementing this suggestion as follows:

    model {
      # same as before
      log_a ~ dnorm(0, .01) 
      log(a) &lt;- log_a
      log_b ~ dnorm(0, .01)
      log(b) &lt;- log_b
      nu &lt;- a
      lambda &lt;- (1/b)^a
      
      for (i in 1:N) {
        # modified to include left-truncation
        CENS[i] ~ dinterval(END_CENS[i], END[i])
        END_CENS[i] ~ dweibull( nu, lambda )T(START[i],)
      }
    }

Unfortunately this doesn&#39;t quite do the trick. With the old code, the model was mostly getting the scale parameter right, but doing a very bad job on the shape parameter. With this new code, it gets very close to the correct shape parameter, but consistently over-estimates the scale parameter. I have noticed that the degree of over-estimation is correlated with how late the split point comes. If the split-point is early (`cens_point = 50`), there&#39;s not really any over-estimation; if it&#39;s late (`cens_point  = 350`), there is a lot.

I thought maybe the problem could be related to &#39;double-counting&#39; the observations: if we see a censored observation at t=300, then from that same person, an uncensored observation at t=400, it seems intuitive to me that this person is contributing two data-points to our inference about the Weibull parameters when really they should just be contributing one point. I, therefore, tried incorporating a random-effect for each person; however, this completely failed, with huge estimates (in the 50-90 range) for the `nu` parameter. I&#39;m not sure why that is, but perhaps that&#39;s a question for a separate post. Since I&#39;m not whether the problems are related, you can find the code for this whole post, including the JAGS code for that model, [here](https://gist.github.com/jwdink/c7db9c69ebfa4fd9b9e7063c0e828259).

All Survival Analysis Solutions on Gang of Coders

Total of 1 Survival Analysis Solutions

Representing Parametric Survival Model in 'Counting Process' form in JAGS