Probabilistic sensitivity analysis and survival models
Implementing probabilistic sensitivity analysis in survival-analysis contexts — relevant for health-economic modelling.
Today we’re going to take a little look into probabilistic sensitivity analysis (PSA), and how it can be implemented within the context of survival analysis. Now PSA is used extensively in health-economic modelling, where uncertain parameters are assigned probability distributions, and values are repeatedly drawn from them, to represent the uncertainty in those parameters. We can then assess the impact of such changes on various predictions from our model, such as survival proportions or (restricted) mean survival time. Or, from a micro-simulation perspective, we are also interested in generating individual survival times from our specified/adapted models – this becomes more relevant within a multi-state context, which is a natural extension.
The mechanics of PSA
To conduct PSA, we need to have an estimated survival model, and have the ability to alter its parameters.
We’ll simulate our own data, recurrence-free survival times on 686 patients with breast cancer, so that we know the truth. About 35% of patients, at random, receive hormonal therapy (hormon), and the times to recurrence, in years, follow a Weibull proportional hazards model, , with scale , shape and log hazard ratio . Each patient is followed for between 1 and 7 years, drawn uniformly, and censored at the end of it. avalon does the simulating, with avalon standard:
. clear
. set seed 686
. set obs 686
Number of observations (_N) was 0, now 686.
. gen hormon = runiform() < 0.35 // hormonal therapy
. gen cens = runiform(1,7) // follow-up, 1 to 7 years
. avalon standard rectime censrec, /// time to recurrence, event
> distribution(weibull) /// Weibull baseline
> lambda(0.1) gamma(1.3) /// scale and shape
> covariates(hormon -0.4) /// log hazard ratio
> maxtime(cens) // censored at end of follow-upand we stset the data:
. stset rectime, failure(censrec)
Survival-time data settings
Failure event: censrec!=0 & censrec<.
Observed time interval: (0, rectime]
Exit on or before: failure
--------------------------------------------------------------------------
686 total observations
0 exclusions
--------------------------------------------------------------------------
686 observations remaining, representing
278 failures in single-record/single-failure data
2,132.133 total analysis time at risk and under observation
At risk from t = 0
Earliest observed entry t = 0
Last observed exit t = 6.986002278 of the 686 patients have an observed recurrence. We use merlin to fit a proportional hazards Weibull model, adjusting for hormonal therapy. The response is stset’s _t, and family(weibull, failure(_d)) names the distribution and the event indicator.
. merlin (_t hormon, family(weibull, failure(_d))), nolog
Fitting full model:
Fixed effects regression model Number of obs = 686
Log likelihood = -825.08363
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
_t: |
hormon | -.521471 .137144 -3.80 0.000 -.7902683 -.2526736
_cons | -2.29398 .1133219 -20.24 0.000 -2.516087 -2.071873
log(gamma) | .2732089 .0512124 5.33 0.000 .1728344 .3735834
------------------------------------------------------------------------------and we get a log hazard ratio of -0.521, against the -0.4 we simulated from, less than one standard error (0.137) away. The baseline parameters are close too: _cons, the log of the scale, is -2.294 (log 0.1 = -2.303), and log(gamma) is 0.273 (log 1.3 = 0.262). Now I’m going to stick to merlin’s own predict in this post: it predicts at whatever times, and for whatever covariate values, we give it, and the same fitted merlin models are what pendragon takes in a more general multi-state setting, so hopefully you can see how to extend it. We first store our model object:
. estimates store m1so that we can get it back after we’ve altered it. I create a time variable at which to calculate predictions,
. gen tvar = 5 in 1
(685 missing values generated)I’m going to predict the survival probability for a patient in the treated group, at t = 5.
We call predict with the survival option, using at(hormon 1) and timevar(tvar).
. predict s5_fit, survival at(hormon 1) timevar(tvar)
(685 missing values generated)Nice and simple. We get one new variable, holding S(5) in the row where tvar is 5:
. list tvar s5_fit in 1
+------------------+
| tvar s5_fit |
|------------------|
1. | 5 .60871659 |
+------------------+which is the probability of still being recurrence free at 5 years, 0.609. As we simulated the data, we know the true value too:
. // the true S(5) for a treated patient
. display exp(-0.1 * 5^1.3 * exp(-0.4))
.580898460.581, so our estimate sits comfortably within sampling error of it: its 95% confidence interval, which we’ll see below, runs from 0.535 to 0.675.
Altering a fitted coefficient
Now let’s say I want to alter my log hazard ratio to assess its impact. Well, I want to set it to be -1, and recalculate my predicted survival at 5 years.
First I’ll make a copy of my estimated coefficient vector, which is stored in e(b).
. matrix b1 = e(b)My log hazard ratio is in the first element, so I’ll replace it with a -1, and check it’s worked:
. matrix b1[1,1] = -1
. matrix list b1
b1[1,3]
_cmp_1_1_1: cons1: dap1_1:
_cons _cons _cons
y1 -1 -2.2939799 .27320891All good. Now the tricky part. There are lots of options for Stata’s maximum likelihood (ml) engine, that most of the time feel like only I know about…but not many spend their time delving into Stata’s source code. But everything here is fully documented! Essentially I want a Weibull model object, with coefficients set to the values in b1. So, we can pass those as initial values using the from() option, but tell it not to proceed with estimation beyond the initial setup step, using iter(0):
. merlin (_t hormon, family(weibull, failure(_d))), from(b1) iter(0)
Fitting full model:
Iteration 0: Log likelihood = -832.15571
convergence not achieved
Fixed effects regression model Number of obs = 686
Log likelihood = -832.15571
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
_t: |
hormon | -1 .1653616 -6.05 0.000 -1.324103 -.6758973
_cons | -2.29398 .1198196 -19.15 0.000 -2.528822 -2.059138
log(gamma) | .2732089 .0558581 4.89 0.000 .1637291 .3826887
------------------------------------------------------------------------------which gives us exactly what we want. We of course get the convergence not achieved message, but our ‘estimated’ parameter vector is set to the values we specified! Check in e(b) as well:
. matrix list e(b)
e(b)[1,3]
_cmp_1_1_1: cons1: dap1_1:
_cons _cons _cons
y1 -1 -2.2939799 .27320891Note, we could have used Ben Jann’s erepost command to make this more elegant, but I prefer to show the steps. So now that we have our merlin model object, with our new hormon coefficient, we can simply predict again:
. predict s5_alt, survival at(hormon 1) timevar(tvar)
(685 missing values generated)
. list tvar s5_alt in 1
+------------------+
| tvar s5_alt |
|------------------|
1. | 5 .73519627 |
+------------------+S(5) goes up from 0.609 to 0.735. predict doesn’t care if your model converged or not, it simply uses the values it finds in e(b). Even if I hadn’t fitted a first model, I could’ve just passed a matrix of my own values to merlin, laid out like e(b), as from() fills them in by position. The only things to make sure you have in memory are the variables the model names: a survival time and an event indicator (here stset’s _t and _d), and a variable called hormon. Those things can be easily created with some dummy data.
Repeating it over a distribution
We can now generalise this to assess the assumption that our treatment effect comes from , by making a large number of draws from the normal distribution, and calculating survival at 5 years for each run, as follows:
. local Nsim = 200
. if _N < `Nsim' set obs `Nsim' // a row for each draw's prediction
. //to store our predictions
. gen s5 = .
(686 missing values generated)
. //loop over Nsim draws
. set seed 13490 // <- don't forget!
. forvalues i=1/`Nsim' {
2. local draw = rnormal(-1,0.1)
3. mat b1[1,1] = `draw'
4. qui merlin (_t hormon, family(weibull, failure(_d))), ///
> from(b1) iter(0)
5. //double check our hormon coefficient matches the random draw
. assert [_cmp_1_1_1][_cons]==`draw'
6. //predict survival at 5 years
. qui predict s5_draw, survival at(hormon 1) timevar(tvar)
7. //store my new estimate
. qui replace s5 = s5_draw[1] in `i'
8. drop s5_draw
9. }which gives us our distribution of estimates of S(5):
. summarize s5
Variable | Obs Mean Std. dev. Min Max
-------------+---------------------------------------------------------
s5 | 200 .7316271 .0231414 .6745485 .7964941Our 200 draws give a mean S(5) of 0.732, with a standard deviation of 0.023, ranging from 0.675 to 0.796. Hopefully that shows the basis for how to manipulate a fitted model in Stata. Any parameter can be altered, and any prediction can be obtained: I simply tricked merlin into turning it back into an estimates object, which means we can use a post-estimation command, such as predict. Here only the treatment effect is varied, from an assumed distribution, with the baseline parameters held at their estimates.
Drawing every parameter with merlin psa
To propagate the estimation uncertainty in all the parameters of a fitted model, draw the whole coefficient vector from a multivariate normal distribution with mean e(b) and variance-covariance matrix e(V), which preserves the correlation between the parameters. merlin psa does exactly that. First, though, a word on what’s in memory: after the loop, the model in memory holds the final draw, not our fitted model. m1 is untouched, as the loop never stored anything, so we restore it:
. display "hormon coefficient in memory: " [_cmp_1_1_1][_cons]
hormon coefficient in memory: -1.097083
. estimates restore m1
(results m1 are active now)
. display "hormon coefficient in memory: " [_cmp_1_1_1][_cons]
hormon coefficient in memory: -.52147096merlin psa then draws the parameter vector 1,000 times from , on the scale merlin estimates on (so it is , not , that is drawn from a normal distribution), and replays each draw through the fitted model. It returns the prediction at e(b), with the 2.5th and 97.5th percentiles of the draws as _lci and _uci; genmean() and gensd() keep the mean and standard deviation of the draws as well:
. merlin psa s5_psa, survival at(hormon 1) timevar(tvar) ///
> reps(1000) seed(4217) ///
> genmean(s5_psa_mean) gensd(s5_psa_sd)
. list s5_psa s5_psa_lci s5_psa_uci s5_psa_mean s5_psa_sd in 1, abbrev(12)
+---------------------------------------------------------------+
| s5_psa s5_psa_lci s5_psa_uci s5_psa_mean s5_psa_sd |
|---------------------------------------------------------------|
1. | .60871659 .53618273 .67115205 .6061554 .03466903 |
+---------------------------------------------------------------+The point estimate is the 0.609 we had from predict, and the draws give a band from 0.536 to 0.671, with a mean of 0.606 and a standard deviation of 0.035. The true 0.581 sits well inside it. This is the estimation uncertainty in all three parameters, rather than an assumed distribution for the treatment effect alone, which is why it is centred on our estimate rather than on a treatment effect of -1. For a model this simple it agrees closely with the delta-method confidence interval from predict’s ci option,
. predict s5_dm, survival at(hormon 1) timevar(tvar) ci
. list s5_dm s5_dm_lci s5_dm_uci in 1, abbrev(10)
+-----------------------------------+
| s5_dm s5_dm_lci s5_dm_uci |
|-----------------------------------|
1. | .60871659 .53472609 .6746006 |
+-----------------------------------+which runs from 0.535 to 0.675. merlin psa can also replay draws you supply, from whatever distributions your analysis assigns: thetadraws() takes a matrix with one row per draw, in the column order of e(b), so the loop above could be a single call. And savedraws() keeps every draw’s prediction, ready to feed into a larger model.
Micro-simulation
So far every prediction has been analytic. But we can use simulation if we so wish! If we’re interested in conducting a micro-simulation and utilising how long each observation spends recurrence free (clearly relevant in a health-economic context), we can use avalon model, which simulates a survival time for every row of the data in memory, from a stored merlin model, using each row’s covariate values. It needs the variables the model was fitted on to be in memory, so rather than start a new dataset I keep one row of ours, put it in the treated group, and copy it 100,000 times. maxtime(5) censors anyone still recurrence free at 5 years, and seed() makes the simulation reproducible:
. keep in 1 // one treated patient...
(685 observations deleted)
. replace hormon = 1
(1 real change made)
. expand 100000 // ...copied 100,000 times
(99,999 observations created)
. avalon model stime event, model(m1) maxtime(5) seed(3418736)Inspecting the simulated times
The simulated times go straight into the data, as new variables, so there is nothing to save and read back in. Let’s take a look:
. list hormon stime event in 1/5
+----------------------------+
| hormon stime event |
|----------------------------|
1. | 1 5 0 |
2. | 1 3.9464051 1 |
3. | 1 3.4894514 1 |
4. | 1 1.5425275 1 |
5. | 1 1.7809361 1 |
+----------------------------+We get the simulated time to recurrence, or to censoring at 5 years, and the event indicator. To calculate our survival probability of interest in this case, we can simply count how many observations are still recurrence free, and divide by the total.
. count if event==0
60,899
. di "S(5) = " r(N)/_N
S(5) = .60899That’s 0.6090, against 0.6087 from the analytic prediction for the same model, m1: pretty close, as they should be. With 100,000 simulated patients the Monte Carlo standard error is about 0.0015, and we should always increase the number simulated until predictions stabilise.
To carry the parameter uncertainty into the micro-simulation as well, add uncertainty(mvn) to avalon model: it draws the parameter vector from and simulates at that draw. By default one draw serves the whole dataset, so repeating the call, with a new seed each time, traces the parameter uncertainty draw by draw; ndraws(#) instead spreads # draws across the rows of a single simulated dataset.
All of the above extends to a more general multi-state setting – fit a merlin model for each transition, and pass them, with your transition matrix, to pendragon for predictions on a fine time grid, where ci cimethod(psa) draws each transition’s parameters from its own , or to avalon msm, which simulates a transition from a fitted merlin model with hazard#(model()).