I came across this tutorial on accelerated time failure survival regression modeling: http://austinrochford.com/posts/2017-10-02-bayes-param-survival.html. I’m interested in the Weibull part.
It’s not clear to me how to model the case with an intercept. See http://personal.psu.edu/drh20/525/weekly/weibull.pdf. The regression model seems to be this one:
\log (T) = \mu + \sigma W + \gamma^\top X
I’m ultimately interested in the Weibull parameters \beta = 1/\sigma (the inverse of the scale) and \eta = \exp(\mu) (the exponential of the intercept). In addition, if I have no predictors, i.e., X = 0, how can I capture that in the Gumbel distribution (not forgetting the intercept \mu)?
Here’s some code, but I don’t think it’s correct:
import pymc3 as pm from theano import shared, tensor as tt y = np.log(time) y_std = (y - y.mean()) / y.std() # Not sure if I need this because of the intercept mu cens = event == 0 # cens_ = shared(cens) # There are no predictors def gumbel_sf(y, mu, sigma): return 1.0 - tt.exp(-tt.exp(-(y - mu) / sigma)) # Build Bayesian model with pm.Model() as model: # http://personal.psu.edu/drh20/525/weekly/weibull.pdf # log(T) = mu + sigmaW + gamma^\top z # where gamma = 0 (no predictors), mu ~ N(0, 10), and sigmaW ~ Gumbel(0, s) with s ~ HalfNormal(5) # Hyperprior s = pm.HalfNormal("s", tau=5.0) # Priors mu = pm.Normal("mu", mu=0, sd=10) # Likelihood for uncensored and censored survival times y_obs = pm.Gumbel("y_obs", mu=mu, beta=s, observed=y_std[~cens]) y_cens = pm.Bernoulli("y_cens", p=gumbel_sf(y_std[cens], mu=mu, sigma=s), observed=np.ones(cens.sum())) # Initialization start = pm.find_MAP() # Perform MC sampling with model: trace = pm.sample(draws=2000, start=start, tune=1000)
How do I obtain \beta and \eta after the sampling is done?