# Survival models in rstanarm

**URL:** <https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998>\
**Category:** Developers\
**Created:** [April 26, 2018, 5:13am UTC](https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998 "2018-04-26T05:13:23Z")\
**Posts on this page:** 20\
**Page:** 4

<div class="post-metadata">

**Author:** ![ermeel](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ermeel/32/4180_2.png) [@ermeel](https://discourse.mc-stan.org/u/ermeel)\
**Post date:** [November 14, 2018, 6:26pm UTC](https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998/61 "2018-11-14T18:26:45Z")

</div>

> [@lcomm](#):
>
> Thus, reasonable priors for the log-baseline hazards are \mathcal{N}(log(E/PT),\left(\log(100)/2 \right)^2), where E is the number of observed events and PT is the total at-risk person-time.

I think log(E/PT) is just the log of the MLE of the baseline hazard (assuming a constant baseline hazard, i.e. survival times have exponential law). It’s a good heuristic, I also sometimes use.

---

<div class="post-metadata">

**Author:** ![ermeel](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ermeel/32/4180_2.png) [@ermeel](https://discourse.mc-stan.org/u/ermeel)\
**Post date:** [November 14, 2018, 6:40pm UTC](https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998/62 "2018-11-14T18:40:42Z")

</div>

@sambrilleman et al., what about posterior predictive checks, i.e. being able to generate (uncensored) survival times? I usually use this in my workflow and found the following algorithm to work well for the M-spline case. It generates _conditioned_ survival times for the “patients” without censoring:

 ![39](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/6/63769b256373d2aa0c1d81fa43076247ecd6ed53.png)

See also [this](https://ermeel86.github.io/case_studies/parametric_surv_models_20181111.html) collection of related snippets I created for a recent meetup.

I call this _conditioned_ posterior predictive checks, because posterior samples a generated conditioned to be less equal the length of study (or maximal event time observed). This is because I’ve seen cases where there was rather large support of the survival distribution on times larger then study length. Here _unconditioned_ or standard posterior predictive checks would indicate a bad model, see [this](https://ermeel86.github.io/case_studies/surv_stan_example_20181028.html) intro-level case study on Bayesian Survival Models in Stan, that I wrote, where I show this for a public dataset.

Lastly, what would we do in the case we would like to generate samples from the posterior predictive distribution without the conditioning. I am a bit concerned about the fact that the i-splines all are zero (or undefined?) beyond the right most knot (in time). One option would be to pick the right most knot large enough in advance (meaning larger than the observed samples), but this workaround seems suboptimal to me.

---

<div class="post-metadata">

**Author:** ![lcomm](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/lcomm/32/1238_2.png) [@lcomm](https://discourse.mc-stan.org/u/lcomm)\
**Post date:** [November 14, 2018, 8:18pm UTC](https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998/63 "2018-11-14T20:18:51Z")

</div>

This is exactly the rationale. You can either think of it as the MLE assuming an exponential distribution or just a simple naive estimator to get a rescaling factor. Truthfully, I even use MLEs for centering prior distributions on the Weibull scale _and_ shape parameters – it’s just more data “double dipping” than I would suggest as a default strategy.

---

<div class="post-metadata">

**Author:** ![lcomm](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/lcomm/32/1238_2.png) [@lcomm](https://discourse.mc-stan.org/u/lcomm)\
**Post date:** [November 14, 2018, 8:45pm UTC](https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998/64 "2018-11-14T20:45:07Z")

</div>

I wish I had seen this earlier. I just spent this afternoon implementing a posterior predictive checks based on Kaplan-Meier survival estimates and it looks like you’ve already done that.

Isn’t this just generating from a truncated event time distribution where the truncation takes the form of an upper bound? Is that what you mean by _conditioned_: you are conditioning on getting to observe the event, so you are forcing it to be less than T^\*?

---

<div class="post-metadata">

**Author:** ![ermeel](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ermeel/32/4180_2.png) [@ermeel](https://discourse.mc-stan.org/u/ermeel)\
**Post date:** [November 14, 2018, 8:49pm UTC](https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998/65 "2018-11-14T20:49:56Z")

</div>

Yes. I do this by rejection sampling, essentially. I’ve had some people suggesting to me to take the minimum of the generated sample and the maximum of observed event times, but this doesn’t seem to be right…

Here is some snippet:

```stan
generated quantities {
    vector[N_uncensored] times_uncensored_sampled;
    {
        real tmp;
        real max_time;
        real max_time_censored;
        max_time = max(times_uncensored);
        max_time_censored = max(times_censored);
        if(max_time_censored > max_time) max_time = max_time_censored;
        
        for(i in 1:N_uncensored) {
            tmp= max_time + 1; 
            while(tmp > max_time) {
                tmp = exponential_rng(exp(intercept+X_uncensored[i,]*betas));
            }
            times_uncensored_sampled[i] = tmp;
        }
    }

}

```

---

<div class="post-metadata">

**Author:** ![lcomm](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/lcomm/32/1238_2.png) [@lcomm](https://discourse.mc-stan.org/u/lcomm)\
**Post date:** [November 14, 2018, 9:39pm UTC](https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998/66 "2018-11-14T21:39:16Z")

</div>

I thought I was missing something because it was hard to invert due to the splines, but maybe not since you are using `exponential_rng` in your example code. To draw an exponential random number with an upper bound `ub` without rejection sampling, you could use the following (untested) code:

```
real exponential_ub_rng(real beta, real ub) {
    real p = exponential_cdf(ub, beta); // cdf for bounds
    real u = uniform_rng(0, p);
    return (-log1m(u) / beta); // inverse cdf for value
  }

```

The `beta` here matches the Stan manual parameterization and would be `exp(intercept+X_uncensored[i,]*betas)` in your notation.

This is somewhat tangential, though. You are presupposing that the same individuals will be censored (or not censored) in any replicated data, and I don’t see why that would be the case. Is it possible that when people are suggesting you take the minimum, they mean to suggest you apply artificial censoring at T^\*? That is what I do. However, it requires knowing a potential censoring time for each individual. For example, if the study ends on Dec 31 and I know they entered the study on June 3, I know exactly when they _would_ have been censored, even if I observe the event to occur in August (and thus do not directly observe the censoring “event” time). Alternately, you may know that everyone will be followed up for a maximum of 90 days. For some censoring mechanisms, this knowledge is not possible, unfortunately.

---

<div class="post-metadata">

**Author:** ![ermeel](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ermeel/32/4180_2.png) [@ermeel](https://discourse.mc-stan.org/u/ermeel)\
**Post date:** [November 15, 2018, 4:17am UTC](https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998/67 "2018-11-15T04:17:23Z")

</div>

> [@lcomm](#):
>
> The `beta` here matches the Stan manual parameterization and would be `exp(intercept+X_uncensored[i,]*betas)` in your notation.

Thanks, this useful!

Your sedond comment is applicable when there is no drop out, right? Otherwise what do you think about the approach to have another event model, that is conditionally independent of the first one, given covariates. One would then use this to describe the censoring / dropout process (inverting event status indicators). In the generative posterior predictive sampling step, one could then use the minimum of the two event times generated, and cap this by the max study length…

---

<div class="post-metadata">

**Author:** ![lcomm](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/lcomm/32/1238_2.png) [@lcomm](https://discourse.mc-stan.org/u/lcomm)\
**Post date:** [November 15, 2018, 11:36am UTC](https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998/68 "2018-11-15T11:36:33Z")

</div>

Yes, I am thinking of administrative censoring, which is the form that my data set has and the easiest to justify as being a noninformative censoring mechanism. Everyone is censored at 90 days if they are still alive to be censored. Maybe a better way of thinking about this distinction is not as _non-informative_ vs. _informative_ as I’ve been doing, but _fixed_ vs. _random_.

All of the posterior predictive options I’ve seen have required some amount of outside information about censoring. I have the easy case where the censoring time C\_i = 90 \ \forall \ i, but colleagues have also worked with random study entry time, which results in known – if individual-specific – values of C\_i regardless of whether censoring is observed to occur. (I will ask my colleagues if there is a draft they could put on arXiv.)

If I’m understanding you correctly, your goal is to develop a posterior predictive approach for _noninformative but random_ censoring? That would certainly enhance the usefulness of posterior predictive checks for survival data. Building a model for the censoring process would be a good start, but I’m not sure how you could deal with the fact that you know your event time model shouldn’t be trusted beyond T^\*. It would be annoying to discover your model fit is “bad” because of predictions you never intended to take seriously anyway. Perhaps that is a justification for building a dropout/“early censoring” model and enforcing an administrative censoring event at T^\* in your replicate data sets. You’d need to make sure any summary measures calculated from the posterior predictions account for censoring.

---

<div class="post-metadata">

**Author:** ![sambrilleman](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sambrilleman/32/1703_2.png) [@sambrilleman](https://discourse.mc-stan.org/u/sambrilleman)\
**Post date:** [November 16, 2018, 5:54am UTC](https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998/69 "2018-11-16T05:54:21Z")

</div>

Cheers @ermeel @lcomm.

> [@ermeel](#):
>
> This one here uses a uniform prior over the simplex, which we could of course change to Dirichlet priors, if required.

Excuse my ignorance, but my understanding was that B-spline basis terms sum to 1. But is it also true that the coefficients on the M-splines must sum to 1? I didn’t realise that. If so, the simplex declaration definitely makes sense. But does this constraint also hold true if the user specifies B-splines on the log baseline hazard?

But if we have a simplex declaration for the spline coefficients, with a Dirichlet prior, it might cause some confusion with the priors specified in the `stan_surv` call. At the moment the `prior_aux` argument refers to the spline coefficients when `basehaz` is set to M-splines or B-splines, and refers to the shape parameter when `basehaz` is Weibull or Gompertz. Of course, for the latter, the Dirichlet prior makes no sense, so it is not clear what we would make the default value for `prior_aux` since the available choices of prior depend on the chosen type of parametric model.

Any thoughts?

> [@ermeel](#):
>
> I think log(E/PT)log(E/PT) is just the log of the MLE of the baseline hazard (assuming a constant baseline hazard, i.e. survival times have exponential law). It’s a good heuristic, I also sometimes use.

Yeah this prior sounds like a good idea for the intercept. In `stan_jm` I basically use something equivalent in which I subtract the log of the crude event rate, i.e. log(E / PT), from the linear predictor, rather than centring the prior on it. It is mathematically equivalent of course, but just means that the user can still define their `prior_intercept` as being centred on zero which is a bit more intuitive I think if, say, they just want to change the prior distribution to `cauchy` or something (and avoids us messing with their choice of prior internally, aside from the autoscaling).

I think centring the intercept (or equivalently, the prior on the intercept) around the crude event rate is probably reasonable and justified for the same reasons **rstanarm** uses the autoscaling. But going as far as centering the shape parameter would definitely start feeling too much like a double use of the data I think. Although it’s not exactly clear where to draw the line! But I think the prior on the shape parameter is far less likely to be problematic compared with the prior on the intercept (famous last words).

> [@ermeel](#):
>
> I call this _conditioned_ posterior predictive checks, because posterior samples a generated conditioned to be less equal the length of study (or maximal event time observed). This is because I’ve seen cases where there was rather large support of the survival distribution on times larger then study length. Here _unconditioned_ or standard posterior predictive checks would indicate a bad model, see [this](https://ermeel86.github.io/case_studies/surv_stan_example_20181028.html) intro-level case study on Bayesian Survival Models in Stan, that I wrote, where I show this for a public dataset.

This seems really cool. In the examples you linked to, the predicted survival times appear more positively skewed than the observed event times, even with the conditioning – is this just because of the model fit? Or is it something related to the approach?

> [@ermeel](#):
>
> Lastly, what would we do in the case we would like to generate samples from the posterior predictive distribution without the conditioning. I am a bit concerned about the fact that the i-splines all are zero (or undefined?) beyond the right most knot (in time). One option would be to pick the right most knot large enough in advance (meaning larger than the observed samples), but this workaround seems suboptimal to me.

Should we even be allowing “unconditional” predictions (i.e. beyond the final event time)? Wouldn’t those simulated event times be drawn from a part of the survival time distribution for which we have no observed data (i.e. we don’t know how the cumulative hazard behaves in that region, yet we are taking it’s estimated value and inverting it), and hence that type of extrapolation seems kinda dangerous? Not to mention the quality of those predictions couldn’t be assessed directly against any observed data anyway?

---

<div class="post-metadata">

**Author:** ![lcomm](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/lcomm/32/1238_2.png) [@lcomm](https://discourse.mc-stan.org/u/lcomm)\
**Post date:** [November 16, 2018, 1:49pm UTC](https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998/71 "2018-11-16T13:49:27Z")

</div>

> [@sambrilleman](#):
>
> I think centring the intercept (or equivalently, the prior on the intercept) around the crude event rate is probably reasonable and justified for the same reasons **rstanarm** uses the autoscaling. But going as far as centering the shape parameter would definitely start feeling too much like a double use of the data I think. Although it’s not exactly clear where to draw the line! But I think the prior on the shape parameter is far less likely to be problematic compared with the prior on the intercept (famous last words).

I once felt that way. I became less principled when I wanted to run a simulation study and everything kept blowing up due to outrageous values for the shape parameter. It is particularly bothersome when specifying multiple hazard models like I have in my application. I agree it feels too much like double dipping for general use, but do not discount the ability of a bad Weibull shape value to totally wreck your model. Shape parameters can be finicky! For me I know I have crossed my personal line because I have to use an iterative procedure to obtain the MLE. If it is a simple rescaling by a formula I can write out like the \log(E/PT), I am okay with it. But if Step 1 of my Bayesian procedure is to verify that my _frequentist_ estimation procedure converged, I feel I am firmly in “double dip” territory.

---

<div class="post-metadata">

**Author:** ![ermeel](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ermeel/32/4180_2.png) [@ermeel](https://discourse.mc-stan.org/u/ermeel)\
**Post date:** [November 18, 2018, 2:31pm UTC](https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998/72 "2018-11-18T14:31:26Z")

</div>

> [@sambrilleman](#):
>
> Excuse my ignorance, but my understanding was that B-spline basis terms sum to 1. But is it also true that the coefficients on the M-splines must sum to 1? I didn’t realise that. If so, the simplex declaration definitely makes sense. But does this constraint also hold true if the user specifies B-splines on the log baseline hazard?

W.l.o.g. assume the covariate vector (per individual) has dimension 1, then we have for the cumulative hazard at time t for an individual with covariate x the following expression

H(t;x;\gamma) = e^{\beta\_0 + \beta\_1 x}\sum\_{k=1}^{m} \gamma\_k I\_k(t)

Now, suppose the spline coefficients do not sum to 1 and denote by 1\neq s=\sum\_{k=1}^m \gamma\_m the sum of spline coefficients. Then we can define normalised spline coefficients \tilde{\gamma}\_k \equiv \gamma\_k / s and obtain

H(t;x;\tilde{\gamma}) = e^{\beta\_0+\log(s) + \beta\_1 x}\sum\_{k=1}^{m} \tilde{\gamma}\_k I\_k(t)

Thus we have an identifiably problem for \beta\_0 and \log(s), if we keep the intercept \beta\_0 in our model. So either we leave the intercept out and use un-normalised spline coefficients or we keep the intercept and use normalised spline coefficients. This is true if we use B- or M- splines on the baseline hazard.

When we use B- or M-splines on the log-hazard scale (think Royston & Parmar), the above argument does not apply (so there is no apparent identifiability problem for the intercept).

> [@sambrilleman](#):
>
> But if we have a simplex declaration for the spline coefficients, with a Dirichlet prior, it might cause some confusion with the priors specified in the `stan_surv` call. At the moment the `prior_aux` argument refers to the spline coefficients when `basehaz` is set to M-splines or B-splines, and refers to the shape parameter when `basehaz` is Weibull or Gompertz. Of course, for the latter, the Dirichlet prior makes no sense, so it is not clear what we would make the default value for `prior_aux` since the available choices of prior depend on the chosen type of parametric model.

What about we use a logit or probit transformation to get from raw coefficients to simplex coefficients. The (unconstrained) normal priors currently used as default in `stan_surv`, would then be for the raw coefficients (see Ch 22.2., Unconstrained Priors: Probit and Logit in the Stan Manual 2.17.)

> [@sambrilleman](#):
>
> This seems really cool. In the examples you linked to, the predicted survival times appear more positively skewed than the observed event times, even with the conditioning – is this just because of the model fit? Or is it something related to the approach?

To be honest, I do not know. I haven’t made up my mind about how to improve this / what the cause is.

> [@sambrilleman](#):
>
> Should we even be allowing “unconditional” predictions (i.e. beyond the final event time)? Wouldn’t those simulated event times be drawn from a part of the survival time distribution for which we have no observed data (i.e. we don’t know how the cumulative hazard behaves in that region, yet we are taking it’s estimated value and inverting it), and hence that type of extrapolation seems kinda dangerous? Not to mention the quality of those predictions couldn’t be assessed directly against any observed data anyway?

Good point, but what would you do in cases where you need to “extrapolate” in an ongoing study using the survival model (potentially in a decision making framework; e.g. adjustments of treatments)? That is you are right into a study (where you haven’t observed many events yet, and have thus mostly censoring) and would like to predict survival times of patients.

Regarding default prior choices for the intercept, I think we should not spend too much time thinking about good default priors, but we should rather encourage (and support) users to do prior predictive checks to find problem specific appropriate priors.

---

<div class="post-metadata">

**Author:** ![sambrilleman](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sambrilleman/32/1703_2.png) [@sambrilleman](https://discourse.mc-stan.org/u/sambrilleman)\
**Post date:** [November 20, 2018, 6:44am UTC](https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998/73 "2018-11-20T06:44:30Z")

</div>

@ermeel @lcomm Hey guys, just to update… today I had a go at making the following changes to `stan_surv`:

1. I centered the predictor matrix.
2. I centered the baseline hazard around the crude event rate (i.e. \log(E/PT)).
3. I removed the intercept (`intercept = TRUE`) from the spline basis when `basehaz = "ms"` or `basehaz = "bs"` and included an intercept parameter in the linear predictor instead.
4. I extended upon 3. to deal with the identifiability issue when `basehaz = "ms"` by taking the raw M-spline coefficients defined on the positive real line, and transforming them to a simplex using Stan’s `softmax` function.

Steps 1. and 2. seemed to work fine. The result is that the prior for `prior_intercept` (for exponential, Weibull and Gompertz) or `prior_aux` (for M-spline and B-spline) corresponds to the (log) baseline hazard when centered around the (log) crude event rate and at the mean values of each of the covariates. I then decided to back-transform the intercept (for exponential, Weibull and Gompertz) or spline coefficients (for M-spline and B-spline) in the generated quantities block before they are returned to the user, so that the centering adjustments aren’t required when making posterior predictions. You can see these changes on [this](https://github.com/stan-dev/rstanarm/compare/feature/survival...center-stansurv-around-crude-log-event-rate) comparison view (in particular, the generated quantities block in the surv.stan file).

However, Steps 3. and 4. seemed to cause problems with the M-spline baseline hazard. Removing the intercept from the spline basis and including an explicit intercept in the linear predictor seems to have a detrimental effect on the sampling. Fitting a simple model to a simulated dataset with one binary covariate and a Weibull baseline hazard, I found that the constrained (simplex) spline coefficients led to sampling times that were about 20 times longer than when I used the unconstrained spline coefficients (with the intercept included in the spline basis). It may be something to do with priors, but I have tried a variety of scale parameters for `prior_intercept` and/or `prior_aux` and couldn’t seem to improve things. The model does return the correct parameter estimates (i.e. the estimated hazard ratio is similar to the value used in the simulation model) and baseline hazard (checked informally using the `ps_check` function), but sampling is very slow. If you want to see how I implemented the constrained spline coefficients, then I have pushed those changes on [this](https://github.com/stan-dev/rstanarm/tree/revise-priors-for-stansurv) branch – but specifically you might want to check out the comparison view that shows how I now generate the [constrained spline coefficients](https://github.com/stan-dev/rstanarm/compare/feature/survival...revise-priors-for-stansurv#diff-20e52bd83eb477b5e070fb2a84bee5bdR605).

Because of the difficulties sampling I’m tempted to just keep the intercept absorbed into the spline basis. The model just seems to perform better, but I can’t pin down why. I think this must have been the behaviour that led me to adopt the `intercept = TRUE` approach in the first place. But, I’m perfectly happy to adopt the other approach if we can get it to work. Since, I do agree, that it is a bit more intuitive for the user to have `prior_intercept` be relevant for all types of baseline hazard.

---

<div class="post-metadata">

**Author:** ![sambrilleman](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sambrilleman/32/1703_2.png) [@sambrilleman](https://discourse.mc-stan.org/u/sambrilleman)\
**Post date:** [November 21, 2018, 3:31am UTC](https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998/74 "2018-11-21T03:31:46Z")

</div>

Just another follow on from this. I did a bit more digging today and it turns out the M-spline model samples fine when using the constrained spline coefficients with a uniform prior over the simplex, as originally proposed by @ermeel.

The sampling difficulties I mentioned yesterday are only encountered when trying to define a prior (e.g. normal, Cauchy, t) on the raw unconstrained spline coefficients, and then transforming them to the constrained space (i.e. a unit simplex) using the `softmax` function.

So… if we opt for an explicit intercept parameter for all types of baseline hazard that will be ok. But the question that remains is what to make the default prior for `prior_aux`. Perhaps the best option is to have the default value for `prior_aux` so-called “missing” in the function signature, e.g.

```no-highlight
stan_surv(formula, data, basehaz = "ms", basehaz_ops, qnodes = 15,
  prior = normal(), prior_intercept = normal(), prior_aux,
  prior_smooth = exponential(autoscale = FALSE), prior_PD = FALSE,
  algorithm = c("sampling", "meanfield", "fullrank"),
  adapt_delta = 0.95, ...)

```

and just describe the default value in the documentation, something like…

> `prior_aux`: The prior distribution for “auxiliary” parameters related to the baseline hazard. The relevant parameters differ depending on the type of baseline hazard specified in the basehaz argument. The following applies (however, for further technical details, refer to the stan\_surv: Survival (Time-to-Event) Models vignette):
> 
> - `basehaz = "ms"`: the auxiliary parameters are the coefficients for the M-spline basis terms on the baseline hazard. These coefficients are defined using a simplex; that is, they are all between 0 and 1, and constrained to sum to 1. This constraint is necessary for identifiability of the intercept in the linear predictor. The default prior is a Dirichlet distribution with a scalar concentration parameter set equal to 1. That is, a uniform prior over all points defined within the support of the simplex. Specifying a scalar concentration parameter \> 1 supports a more even distribution (i.e. a smoother spline function), while specifying a scalar concentration parameter \< 1 supports a more sparse distribution (i.e. a less smooth spline function).
> - `basehaz = "bs"`: the auxiliary parameters are the coefficients for the B-spline basis terms on the log baseline hazard. These parameters are unbounded. The default is a normal distribution with mean 0 and scale 20.
> - `basehaz = "exp"`: there is no auxiliary parameter, since the log scale parameter for the exponential distribution is incorporated as an intercept in the linear predictor.
> - `basehaz = "weibull"`: the auxiliary parameter is the Weibull shape parameter, while the log scale parameter for the Weibull distribution is incorporated as an intercept in the linear predictor. The auxiliary parameter has a lower bound at zero. The default is a half-normal distribution with mean 0 and scale 2.
> - `basehaz = "gompertz"`: the auxiliary parameter is the Gompertz scale parameter, while the log shape parameter for the Gompertz distribution is incorporated as an intercept in the linear predictor. The auxiliary parameter has a lower bound at zero. The default is a half-normal distribution with mean 0 and scale 2.
> 
> Currently, `prior_aux` can be a call to dirichlet, normal, student\_t, cauchy, or exponential. See `priors` for details on these functions. Note that not all prior distributions are allowed with all types of baseline hazard. To omit a prior —i.e., to use a flat (improper) uniform prior— set `prior_aux` to `NULL`.

Maybe @bgoodri also has some advice on this.

---

<div class="post-metadata">

**Author:** ![ermeel](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ermeel/32/4180_2.png) [@ermeel](https://discourse.mc-stan.org/u/ermeel)\
**Post date:** [November 25, 2018, 6:39am UTC](https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998/75 "2018-11-25T06:39:13Z")

</div>

@sambrilleman, I think the documentation for `prior_aux` is fine. I propose we stick with it and see how it is received by users and in case it causes too much confusion, we change it.

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [November 25, 2018, 4:16pm UTC](https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998/76 "2018-11-25T16:16:15Z")

</div>

That seems fine, but I would do `prior_aux = dirichlet()` in the signature. We already have that in **rstanarm** for variance components.

---

<div class="post-metadata">

**Author:** ![sambrilleman](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sambrilleman/32/1703_2.png) [@sambrilleman](https://discourse.mc-stan.org/u/sambrilleman)\
**Post date:** [November 29, 2018, 7:06am UTC](https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998/77 "2018-11-29T07:06:26Z")

</div>

Fyi, I’ve just merged AFT models into the [survival](https://github.com/stan-dev/rstanarm/tree/feature/survival) branch. Covariate effects can be incorporated using time-fixed or time-varying acceleration factors. Although, at the moment only an exponential or Weibull can be specified for the baseline survival distribution. I’ve added a bunch of details on the model formulations to the vignette.

> [@bgoodri](#):
>
> That seems fine, but I would do `prior_aux = dirichlet()` in the signature. We already have that in **rstanarm** for variance components.

The difficulty with doing this is that if the user changes the baseline hazard from the default setting (which they are quite likely to do), for example, they specify `basehaz = "weibull"`, then the Dirichlet prior isn’t valid and they will get an error. So it is a bit of pain. Every time they change the baseline hazard they would have to change `prior_aux` too.

I’ve actually been playing around with the M-spline model more, based on the simplex and Dirichlet prior for the M-spline coefficients. But it does something that surprises (worries?) me… Increasing the number of knots in the M-spline function for the baseline hazard doesn’t necessarily increase the flexibility of the curve. For example, say I increase the number of knots from 5 to 50, I would expect the curve to be much more “wiggly”. Instead, it isn’t. In the examples I ran, the Dirichlet prior forces most of the coefficients to zero. In other words, the Dirichlet prior seems to enforce smoothing. Is this expected? If so, that seems ok. But I didn’t realise that the Dirichlet would act as a smoothing prior.

---

<div class="post-metadata">

**Author:** ![ermeel](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ermeel/32/4180_2.png) [@ermeel](https://discourse.mc-stan.org/u/ermeel)\
**Post date:** [December 1, 2018, 9:01am UTC](https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998/78 "2018-12-01T09:01:09Z")

</div>

> [@sambrilleman](#):
>
> In the examples I ran, the Dirichlet prior forces most of the coefficients to zero. In other words, the Dirichlet prior seems to enforce smoothing. Is this expected? If so, that seems ok. But I didn’t realise that the Dirichlet would act as a smoothing prior.

That is an interesting observation! I do not have an explanation for this, but maybe one can understand this by considering how “probability-mass” is distributed in high-dimensional simplexes through the Dirichlet, maybe @betanalpha has an intuition here?

> [@sambrilleman](#):
>
> For example, say I increase the number of knots from 5 to 50, I would expect the curve to be much more “wiggly”.

Regarding the expected “wiggliness”: This might also be less severe for this particular case, since we have M-splines here that are non-negative and are also combined using non-negative weights. I would expect (this is a very vague intuition) that the “wiggliness” is therefore suppressed, since we cannot have cancellation effects of splines and coefficients…

Out of curiosity, did you find different choices of priors for the m-/i-spline coefficients that can lead to a very “wiggly” curve, as the number of knots is increased?

---

<div class="post-metadata">

**Author:** ![betanalpha](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/betanalpha/32/4_2.png) [@betanalpha](https://discourse.mc-stan.org/u/betanalpha)\
**Post date:** [December 2, 2018, 2:31am UTC](https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998/79 "2018-12-02T02:31:25Z")

</div>

> [@ermeel](#):
>
> That is an interesting observation! I do not have an explanation for this, but maybe one can understand this by considering how “probability-mass” is distributed in high-dimensional simplexes through the Dirichlet, maybe @betanalpha has an intuition here?

The behavior of the Dirichlet distribution as the dimensionality increases depends on the configuration of the hyperparameters. Let’s assume that they’re all set to the same value, \alpha. For \alpha \le 1 then as the dimensionality increases the Dirichlet will concentrate around the corners where most of the values hover around zero. For \alpha \> 1 there will be a soft concentration in the center of the simplex.

---

<div class="post-metadata">

**Author:** ![ermeel](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ermeel/32/4180_2.png) [@ermeel](https://discourse.mc-stan.org/u/ermeel)\
**Post date:** [December 2, 2018, 6:44am UTC](https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998/80 "2018-12-02T06:44:56Z")

</div>

Thanks!

Here are some examples with 13 knots and different values of \alpha. The figures show prior samples of the baseline hazard curve and induced baseline survival curves, that is

h(t) = e^{\beta\_0} \sum\_{k=1}^{m} \gamma\_k M\_k(t)

and

S(t) = \exp\left(-e^{\beta\_0} \sum\_{k=1}^{m} \gamma\_k I\_k(t) \right)

where the gammas have a joint Dirichlet prior with common \alpha parameter and \beta\_0 has a Normal(-5,2) prior. The M\_k,I\_k are, respectively, the corresponding m-/i-spline basis functions at time point t. Each of the figures shows 8 samples.

\alpha=1:

 ![alpha_1](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/b/b730921734f2abcc4cd1570273ec0e603aae2c0c.png)

\alpha=2:

 ![alpha_2](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/e/ebb09f8da125abd35aab470729ae243429b994c9.png)

\alpha=4:

 ![alpha_4](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/7/791f053e55f24103eed4ea41fba0566e00053b79.png)

\alpha=8:

 ![alpha_8](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/1/1a00a7b52f627bf30c1fe006713a36d941779732.png)

\alpha=0.5

 ![5](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/2/2a32c6212ad56380773120a1b949ad7401bbb593.png)

\alpha=0.25:

 ![25](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/8/8f312f675ba8cc2112f0cfd3cc31d525b41d3f0e.png)

\alpha=0.125

 ![125](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/0/0f03e23c0d3de5c16e17a02dc6c8e8adde4f7a6b.png)

I guess what it is happening for low values of \alpha, we effectively see a randomly chosen basis function here, where as for large values of \alpha we see, roughly, an average basis function…

To quote [Wikipedia](https://en.m.wikipedia.org/wiki/Concentration_parameter):

> A concentration parameter of 1 (or _k_ , the dimension of the Dirichlet distribution, by the definition used in the topic modelling literature) results in all sets of probabilities being equally likely, i.e. in this case the Dirichlet distribution of dimension _k_ is equivalent to a uniform distribution over a [_k-1_ -dimensional simplex](https://en.m.wikipedia.org/wiki/Standard_simplex). Note that this is _not_ the same as what happens when the concentration parameter tends towards infinity. In the former case, all resulting distributions are equally likely (the distribution over distributions is uniform). In the latter case, only near-uniform distributions are likely (the distribution over distributions is highly peaked around the uniform distribution). Meanwhile, in the limit as the concentration parameter tends towards zero, only distributions with nearly all mass concentrated on one of their components are likely (the distribution over distributions is highly peaked around the _k_ possible [Dirac delta distributions](https://en.m.wikipedia.org/wiki/Dirac_delta_distribution) centered on one of the components, or in terms of the _k_ -dimensional simplex, is highly peaked at corners of the simplex).

---

<div class="post-metadata">

**Author:** ![lcomm](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/lcomm/32/1238_2.png) [@lcomm](https://discourse.mc-stan.org/u/lcomm)\
**Post date:** [December 5, 2018, 5:29pm UTC](https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998/81 "2018-12-05T17:29:25Z")

</div>

Let’s say you have a \mathrm{Dirichlet}(\alpha\_1, \dots, \alpha\_k) prior. According to [this](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3081791/) paper, the effective sample size of that prior information is \sum\_i^{k} \alpha\_i. So if you keep \alpha\_1 = \alpha\_2 = \dots = \alpha\_k constant as you increase k you are effectively asserting a stronger a priori belief that the basis weights are equal. Is that desirable?

[Previous page](https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998.md?page=3)

[Next page](https://discourse.mc-stan.org/t/survival-models-in-rstanarm/3998.md?page=5)