Review data block

Let’s review the data block first.

This will tell us the structure and format of data input to the model.

print_stan_file(stan_file, section = 'data')

## data { 
##   int<lower=0> Nobs; 
##   int<lower=0> Ncen; 
##   vector[Nobs] yobs; 
##   vector[Ncen] ycen; 
## }

The censored & observed data points are provided as separate input vectors.

observed data points

• Nobs: number of observed data points
• yobs: times to observed events

censored data points

• Ncen: number of censored data points
• ycen: times to censored events

Recall that this is a NULL model (with no covariate values), so the number & values of observed covariates are not needed.

Review model block

The stan code contains an implicit constant term, in the linear predictor mu.

print_stan_file(stan_file, section = 'model')

## model { 
##   yobs ~ weibull(alpha, exp(-(mu)/alpha)); 
##   increment_log_prob(weibull_ccdf_log(ycen, alpha, exp(-(mu)/alpha))); 
##  
##   alpha_raw ~ normal(0.0, 1.0); 
##   mu ~ normal(0.0, tau_mu); 
## }

Observe how the ccdf (complementary cumulative distribution function) is used to compute the log probability of the censored observations.

What does the ccdf represent in this scenario?

How does the model address the censoring process?

Review parameters block

Our stan code also contains a reparameterization of the alpha term, in the transformed parameters block.

Observe:

print_stan_file(stan_file, section = 'transformed parameters')

## transformed parameters { 
##   real alpha; 
##   alpha <- exp(tau_al * alpha_raw); 
## }

(recall that tau_al is a constant scaling term set to 10, and alpha_raw is a parameter with a normal(0, 1) prior distribution).

This reparameterization achieves two things :

1. The use of tau_al * alpha_raw is an example of a non-centered parameterization.
   • It would have been mathematically equivalent to define a (non-transformed) parameter alpha with a prior normal(0, 10).
   • However, this parameterization yields a parameter (alpha_raw) which is on a similar scale as other parameters in our model. The exp() transformation makes the difference between these two scales even more dramatic.
   • In general, having all parameters on a similar scale makes the sampling more efficient.
2. The exp() transformation of this parameter additionally allows us to put a prior on log alpha.
   • we want to put a prior on log alpha since alpha enters into our likelihood in the exponent.

This seems like a lot of gymnastics to be doing.

However, it has practical implications for our modeling efficiency.

Observe that, for a single value of alpha_raw (e.g. 0.2), the transformation yields:

alpha_raw <- 0.2
tau_al <- 10
log_alpha <- alpha_raw * tau_al
alpha <- exp(log_alpha)
print(alpha)

## [1] 7.389056

which may seem silly.

However

Consider the resulting distribution of alpha over a range of values for alpha_raw sampled from our normal(0, 1) prior:

alpha_raw <- rnorm(1000, 0, 1)
tau_al <- 10
log_alpha <- alpha_raw * tau_al
alpha <- exp(log_alpha)
ggplot(data.frame(alpha = alpha, alpha_raw = alpha_raw), 
       aes(x = alpha)) + 
    geom_density() + 
    scale_x_log10(labels = scientific)

Notice how alpha ranges from 1e-10 to 1e+10 on a log scale. We have to truncate this dramatically to even consider plotting it on its original scale.

Sampling this parameter space may require different step sizes & different tuning parameter values throughout this distribution.

The alpha_raw scale, by comparison, is a lot friendlier.

ggplot(data.frame(alpha = alpha, alpha_raw = alpha_raw), 
       aes(x = alpha, y = alpha_raw)) + 
    geom_density2d() + 
    scale_x_log10(labels = scientific)

This distribution is centered at 0 and has more consistent behavior throughout its range of values.

What’s important to note here is that while the non-centered parameterization is mathematically equivalent to the standard parameterization, it is (in some ways) a different model. Consider that the reparameterization will impact the posterior estimates of parameter values.

Packages like rstanarm which provide easy wrappers to a variety of standard models implemented in Stan use a non-centered parameterization by default.

More information on non-centered parameterization:

1. discussion on stan-dev list
2. Gelman, 2004. Parameterization and Bayesian Modeling

Simulate data for arbitrary input values

We can use this function to simulate a dataset for hypothetical parameter values of alpha & mu.

test_alpha <- 0.8
test_mu <- -3

## sample size from TCGA blca data
test_n <- nrow(clin_data)

## test these inputs for arbitrary values of alpha & mu
simulated_data <- 
    weibull_sim_data(alpha = test_alpha,
                 mu = test_mu,
                 n = test_n
                 ) 
head(simulated_data)

##   surv_months censor_months os_status   os_months
## 1   0.5995253      15.65978  DECEASED   0.5995253
## 2  24.4544603      62.24349  DECEASED  24.4544603
## 3   5.9572821      87.10862  DECEASED   5.9572821
## 4  29.3960125      32.69540  DECEASED  29.3960125
## 5  19.1953454      94.13512  DECEASED  19.1953454
## 6 207.2250962     343.31931  DECEASED 207.2250962

The simulated survival curve looks like:

## plot KM curve from simulated data
simulated_data <- 
    simulated_data %>%
    dplyr::mutate(os_deceased = os_status == 'DECEASED')

autoplot(survival::survfit(Surv(os_months, os_deceased) ~ 1,
                      data = simulated_data
                      ), conf.int = F) + 
    ggtitle('Simulated KM curve')

preparing data for stan

Stan takes data input as a list. The contents of the list should match those of the data block in the stan code.

E.g. looking at the data block -

print_stan_file(stan_file, section = 'data')

## data { 
##   int<lower=0> Nobs; 
##   int<lower=0> Ncen; 
##   vector[Nobs] yobs; 
##   vector[Ncen] ycen; 
## }

our input list to Stan should contain dimensions & values for observed & censored data, separately.

observed_data <- simulated_data %>%
    dplyr::filter(os_status == 'DECEASED')

censored_data <- simulated_data %>%
    dplyr::filter(os_status != 'DECEASED')

stan_data <- list(
    Nobs = nrow(observed_data),
    Ncen = nrow(censored_data),
    yobs = observed_data$os_months,
    ycen = censored_data$os_months
)
rm(censored_data)
rm(observed_data)
str(stan_data)

## List of 4
##  $ Nobs: int 292
##  $ Ncen: int 117
##  $ yobs: num [1:292] 0.6 24.45 5.96 29.4 19.2 ...
##  $ ycen: num [1:117] 33.17 9.95 17.57 6.31 150.28 ...

(wrap this prep-data process in a function gen_stan_data for later)

gen_stan_data <- function(data) {
    observed_data <- data %>%
        dplyr::filter(os_status == 'DECEASED')
    
    censored_data <- data %>%
        dplyr::filter(os_status != 'DECEASED')
    
    stan_data <- list(
        Nobs = nrow(observed_data),
        Ncen = nrow(censored_data),
        yobs = observed_data$os_months,
        ycen = censored_data$os_months
    )
}

test simulated values with stan

Let’s call stan:

recover_simulated <- 
    rstan::stan(stan_file,
                data = gen_stan_data(simulated_data),
                chains = 4,
                iter = 1000,
                seed = 1328025050
                )
print(recover_simulated)

## Inference for Stan model: weibull_survival_null_model.
## 4 chains, each with iter=1000; warmup=500; thin=1; 
## post-warmup draws per chain=500, total post-warmup draws=2000.
## 
##                    mean      se_mean           sd          2.5%
## alpha_raw  6.000000e-02 9.000000e-02 1.300000e-01 -3.000000e-02
## mu        -2.630000e+00 7.100000e-01 1.010000e+00 -3.500000e+00
## alpha      5.210000e+00 5.350000e+00 7.570000e+00  7.600000e-01
## lp__      -1.220654e+45 1.494616e+45 2.114767e+45 -4.892679e+45
##                     25%      50%      75%    97.5% n_eff    Rhat
## alpha_raw -2.000000e-02    -0.02     0.07     0.29     2   33.82
## mu        -3.270000e+00    -3.14    -2.17    -0.89     2    7.38
## alpha      8.200000e-01     0.86     5.31    18.32     2  228.85
## lp__      -1.217518e+45 -1400.52 -1399.71 -1399.34     2 1240.90
## 
## Samples were drawn using NUTS(diag_e) at Wed Jun 22 22:20:40 2016.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

What’s wrong with this picture?

(A: poor convergence) (A: in some chains, we see a lot of numerical problems.)

Setting initial values

This step is usually optional, but may be necessary for some models.

In this case, it may be useful to set initial values. Recall the projected range of our transformed parameter alpha?

By default, Stan chooses a random initial value for each parameter on the unconstrained scale between -2 and 2. This random initialization is on the unconstrained support for each parameter. This guarantees that initial values are consistent with the constrained range.

When we pass the initial values in, however, these are on the constrained scale. See the Stan manual for more details about transformations applied to constrained variables.

print_stan_file(stan_file, section = 'parameters')

## parameters { 
##   real alpha_raw; 
##   real mu; 
## }

gen_inits <- function() {
      list(
        alpha_raw = 0.01*rnorm(1),
        mu = rnorm(1)
      )
}

stan code with initial values

Let’s try fitting our stan model again with our initial values function.

recover_simulated2 <- 
    rstan::stan(stan_file,
                data = gen_stan_data(simulated_data),
                chains = 4,
                iter = 1000,
                init = gen_inits
                )
print(recover_simulated2)

## Inference for Stan model: weibull_survival_null_model.
## 4 chains, each with iter=1000; warmup=500; thin=1; 
## post-warmup draws per chain=500, total post-warmup draws=2000.
## 
##               mean se_mean   sd     2.5%      25%      50%      75%
## alpha_raw    -0.02    0.00 0.00    -0.03    -0.02    -0.02    -0.01
## mu           -3.20    0.01 0.17    -3.55    -3.31    -3.20    -3.09
## alpha         0.84    0.00 0.04     0.76     0.81     0.84     0.87
## lp__      -1400.41    0.07 1.17 -1403.77 -1400.84 -1400.03 -1399.58
##              97.5% n_eff Rhat
## alpha_raw    -0.01   229 1.00
## mu           -2.89   224 1.00
## alpha         0.92   233 1.00
## lp__      -1399.33   311 1.01
## 
## Samples were drawn using NUTS(diag_e) at Wed Jun 22 22:20:45 2016.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

Now we see fewer numerical problems, and better R-hat values.

checking convergence

Assessing convergence can be a tricky business, since every model & every scenario is different.

In general, I look at three things when assessing convergence :

1. Rhat values (are they close to 1)?
2. Review traceplots for lp__ & key parameters
3. launch shinystan for further checking.

rstan::traceplot(recover_simulated2, 'lp__')

rstan::traceplot(recover_simulated2, c('alpha','mu'), ncol = 1)

if (interactive())
    shinystan::launch_shinystan(recover_simulated2)

extracting parameters from the Stanfit object

We use the rstan::extract() function to extract parameters from the stanfit object.

E.g. to extract alpha & mu:

pp_alpha <- rstan::extract(recover_simulated2,'alpha')$alpha
pp_mu <- rstan::extract(recover_simulated2,'mu')$mu

Each of these is a 1xD vector of values, where D = the number of posterior (post-warmup) draws.

In this case, we have 2000: 4 chains * 1000 iterations / 2

For example, we can plot the posterior distribution of alpha, in context of test_alpha:

ggplot(data.frame(alpha = pp_alpha, mu = pp_mu)) + 
    geom_density(aes(x = alpha)) + 
    geom_vline(aes(xintercept = test_alpha), colour = 'red') +
    ggtitle('Posterior distribution of alpha\nshowing true value in red')

and, repeat the same for mu:

ggplot(data.frame(alpha = pp_alpha, mu = pp_mu)) + 
    geom_density(aes(x = mu)) + 
    geom_vline(aes(xintercept = test_mu), colour = 'red') +
    ggtitle('Posterior distribution of mu\nshowing true value in red')

The posterior estimates of the parameters do contain those used to simulate our data, but they are not necessarily at the mode of each distribution.

Also, we have a high degree of correlation between mu and alpha:

ggplot(data.frame(alpha = pp_alpha, mu = pp_mu)) + 
    geom_density2d(aes(x = alpha, y = mu)) +
    geom_point(aes(x = test_alpha, y = test_mu), colour = 'red', size = 2) +
    ggtitle('Posterior distributions of mu and alpha\nshowing true parameter values in red')

It does seem like these values are on the edge of the posterior density for the parameter values.

How likely are the parameter values used to simulate our data, according to this model?

One nice thing about sampling from the full posterior is that we can compute posterior probabilities directly, by summarizing over the values.

To compute the probability of seeing a value of alpha >= 0.8:

mean(pp_alpha >= test_alpha)

## [1] 0.841

And for mu (testing mu >= -3:

mean(pp_mu >= test_mu)

## [1] 0.12

The joint probability of seeing this combination of parameter values, however, isn’t as encouraging:

mean(pp_mu >= test_mu & pp_alpha >= test_alpha)

## [1] 0.018

We would need more iterations in order to estimate this with greater precision.

simulating data for each posterior draw

We can use hadley’s purrr::map2 to simulate data for each pair of mu*alpha values.

pp_newdata <- 
    purrr::map2(.x = pp_alpha,
                .y = pp_mu,
                .f = ~ weibull_sim_data(alpha = .x, 
                                mu = .y,
                                n = test_n
                                )
                )

If you’re not familiar with purrr, I recommend you inspect the resulting object.

We now have a list of D datasets, each containing a simulation according to that draw’s parameter values for mu & alpha.

Let’s plot the time to event in the posterior draws, and compare this to the test dataset we used to fit our model.

ggplot(pp_newdata %>%
           dplyr::bind_rows() %>%
           dplyr::mutate(type = 'posterior predicted values') %>%
           bind_rows(simulated_data %>% dplyr::mutate(type = 'actual data'))
       , aes(x = os_months, group = os_status, colour = os_status, fill = os_status)) +
    geom_density(alpha = 0.5) +
    facet_wrap(~type, ncol = 1)

These look pretty similar.

summarizing posterior predictive draws

Next we might ask about the posterior estimates of the survival curve. How would we estimate this?

One way (there may be several) is to:

1. compute the cumulative survival at each observed timepoint for each draw from the posterior
2. aggregate the cumulative survival estimates to discrete units of time
3. summarize the cumulative survival for each interval, over the posterior draws.

This is the method we will use here.

## cumulative survival rate for each posterior draw
pp_survdata <-
    pp_newdata %>%
    purrr::map(~ dplyr::mutate(., os_deceased = os_status == 'DECEASED')) %>%
    purrr::map(~ survival::survfit(Surv(os_months, os_deceased) ~ 1, data = .)) %>%
    purrr::map(fortify)

## summarize cum survival for each unit time (month), summarized at 95% confidence interval
pp_survdata_agg <- 
    pp_survdata %>%
    purrr::map(~ dplyr::mutate(., time_group = floor(time))) %>%
    dplyr::bind_rows() %>%
    dplyr::group_by(time_group) %>%
    dplyr::summarize(surv_mean = mean(surv)
                     , surv_p50 = median(surv)
                     , surv_lower = quantile(surv, probs = 0.025)
                     , surv_upper = quantile(surv, probs = 0.975)
                     ) %>%
    dplyr::ungroup()

Finally, we overlay the posterior predictive simulations of the survival curve with that from our original test dataset.

## km-curve for test data 
test_data_kmcurve <- 
    fortify(
        survival::survfit(
            Surv(os_months, os_deceased) ~ 1, 
            data = simulated_data %>% 
                dplyr::mutate(os_deceased = os_status == 'DECEASED')
            )) %>%
    dplyr::mutate(lower = surv, upper = surv)

ggplot(pp_survdata_agg %>%
           dplyr::mutate(type = 'posterior predicted values') %>%
           dplyr::rename(surv = surv_p50, lower = surv_lower, upper = surv_upper, time = time_group) %>%
           bind_rows(test_data_kmcurve %>% dplyr::mutate(type = 'actual data')),
       aes(x = time, group = type, linetype = type)) + 
    geom_line(aes(y = surv, colour = type)) +
    geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2) +
    xlim(c(0, 200))

Here we see that the survival curve from our original (simulated) data matches the posterior predictive density pretty closely. This is a good thing since we are working with simulated data!

Saving as a function

As before, we will want to wrap this in a function so that it can be reused in future steps, e.g. when we work with our TCGA data.

pp_predict_surv <- function(pp_alpha, pp_mu, n,
                            level = 0.9,
                            plot = F, data = NULL,
                            sim_data_fun = weibull_sim_data
                            ) {
    pp_newdata <- 
        purrr::map2(.x = pp_alpha,
                    .y = pp_mu,
                    .f = ~ sim_data_fun(alpha = .x, mu = .y, n = n)
                    )
    
    pp_survdata <-
        pp_newdata %>%
        purrr::map(~ dplyr::mutate(., os_deceased = os_status == 'DECEASED')) %>%
        purrr::map(~ survival::survfit(Surv(os_months, os_deceased) ~ 1, data = .)) %>%
        purrr::map(fortify)
    
    ## compute quantiles given level 
    lower_p <- 0 + ((1 - level)/2)
    upper_p <- 1 - ((1 - level)/2)
    
    pp_survdata_agg <- 
        pp_survdata %>%
        purrr::map(~ dplyr::mutate(.,
                                   time_group = floor(time))) %>%
        dplyr::bind_rows() %>%
        dplyr::group_by(time_group) %>%
        dplyr::summarize(surv_mean = mean(surv)
                         , surv_p50 = median(surv)
                         , surv_lower = quantile(surv,
                                                 probs = lower_p)
                         , surv_upper = quantile(surv,
                                                 probs = upper_p)
                         ) %>%
        dplyr::ungroup()
    
    if (plot == FALSE) {
        return(pp_survdata_agg)
    } 
    
    ggplot_data <- pp_survdata_agg %>%
           dplyr::mutate(type = 'posterior predicted values') %>%
           dplyr::rename(surv = surv_p50,
                         lower = surv_lower,
                         upper = surv_upper, time = time_group)
    
    if (!is.null(data))
        ggplot_data <- 
            ggplot_data %>% 
            bind_rows(
                fortify(
                    survival::survfit(
                        Surv(os_months, os_deceased) ~ 1, 
                        data = data %>% 
                            dplyr::mutate(
                                os_deceased = os_status == 'DECEASED')
                        )) %>%
                dplyr::mutate(lower = surv,
                              upper = surv, type = 'actual data')
                )
    
    pl <- ggplot(ggplot_data,
                 aes(x = time, group = type, linetype = type)) + 
        geom_line(aes(y = surv, colour = type)) +
        geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2)
        
    pl 
}