Tutorial environment, package and extra code loading, and random seed setups
usingPkg, NTRsurvpkg_dir =pkgdir(NTRsurv)Pkg.activate(joinpath(pkg_dir, "examples"))Pkg.instantiate()#Load packages and extra code needed for tutorialusingRandom, Distributions, Survival, StatsModels, DataFrames, HCubature, NTRsurv, Plots, Measures, LaTeXStrings# Code for frequentist confidence bandsfrequentist_confidence_bands_path =joinpath( dirname(dirname(pathof(NTRsurv))), "examples", "frequentist_confidence_bands.jl")include(frequentist_confidence_bands_path)# Set seed for reproducibilityRandom.seed!(1234)
Weibull Simulation Study
This notebook provides an in-depth presentation of a Neutral-to-the-Right (NTR) survival modeling workflow. We review key theoretical results from the literature, including posterior consistency at the optimal rate relative to the Kaplan–Meier estimator and a Bernstein–von Mises result for the Cox proportional hazards model, which justify the modeling approach. We also describe efficient posterior simulation schemes and demonstrate a full implementation using synthetic Weibull data with the Julia package NTRsurv.jl.
For a more hands-on introduction focused on applied data analysis, see the companion tutorials that emphasize practical modeling without extensive theoretical development.
NTR distributions for Kaplan-Meier alternative
We say that a random probability distribution supported on \(\mathbb{R}^+\), with cumulative distribution function \(F\), is neutral to the right if for positive real numbers \(0\leq t_1<\ldots <t_n\) the random vector
Such distributions can be characterized by having survival function
\[
S(t)=1-F(t)=\exp\left(-Y(t) \right)
\]
where \(Y(t)\) is an almost surely non-decreasing function with independent increments, known in the stochastic processes literature as a subordinator. We say that a random variable with survival function \(S\) as above has a \(\text{NTR}(Y)\) distribution. The so called subordinator \(Y\) is in turn characterized by a Poisson process \(N\) with intensity \(\nu(\mathrm{d}w,\mathrm{d}u)\) in \(\mathbb{R}^+\times\mathbb{R}^⁺\), the “w” axis corresponding to weights in the jumps of \(Y\) while the “u” axis corresponds to locations of those jumps, such that we can express
where the so called Lévy intensity satisfies the integrability condition ${ [0,t] +} (1 - e^{-w} )(w,u) < $ for any \(t>0\). Such subordinators have an almost sure representation
for some r.v.’s \(U_i,W_i\stackrel{\text{a.s.}}{>}0\) determined by the intensity \(\nu\) such that \(Y(t)\) is a.s. finite for any \(t>0\).
Particular choices of subordinators can have know distribution of their increments, usually providing easier agorithms for stochastic simulation. One such example is the following Gamma process.
NTR with Gamma processes workflow for survival analysis
The Laplace exponent of a Gamma process is given by \[\begin{align*}
\psi_t(\lambda) = -\log\left( \mathbb{E}\left[ e^{-\lambda Y(t)} \right] \right) = \beta \log\left( 1 +\frac{\lambda}{\alpha} \right) \kappa(t).
\end{align*}\]
We choose \[\beta = \frac{1}{\log\left( 1+\frac{1}{\alpha} \right)}\] with which
so \(\kappa\) serves as a cumulative hazard baseline for centering the model a priori.
On the other hand \[
\mathbb{E}\left[ (S(t))^2 \right]=e^{-\psi_t(2)} = exp\left( -\frac{\log\left(1+\frac{2}{\alpha} \right)}{\log\left(1+\frac{1}{\alpha} \right)} \kappa(t) \right)
\] and
For a Bayesian analysis we can work with the reparaetrization \(\tilde{\alpha}=\frac{1}{\alpha}\) as it is simpler to set a prior distribution on \(\tilde{\alpha}\) which concentrates around \(0\). Alternatively, in the proposed framework \(\alpha\) is fixed and should reflect the a priori believes on the variance of the centering baseline survival.
We recommend drawing prior simulations of survival curves to assess the a priori choices of the NTR model. For example, we visualize prior behavior under an unit rate exponential baseline, corresponding to cumulative hazard \(\kappa(t) = t\), and variance modulating hyperparameter \(\alpha=5\). We sample draws from the induced prior survival curves on a fixed time grid.
with \(X_1,\dots,X_n\) i.i.d. nonnegative random variables representing the event times of interest and \(C_1,\dots,C_n\) i.i.d. nonnegative censoring times, independent of the \(X_i\)’s. We denote the distinct order statistics for observed survival data \(\left\{ T_i=t_i\right\}_{i=1}^n\) as \(0 < t_{(1)} < t_{(2)} < \cdots < t_{(k_n)}\), with $k_n $ the observed number of distinct values. At each event time \(t_{(k)}\) we define the number of exact observations at time \(t_{(k)}\) as \[
d_k = \sum_{i=1}^n \mathbf{1}\left[ \delta_i=1,\;T_i = t_{(k)}\right\},
\] the number of individuals at risk just prior to \(t_{(k)}\) as \[
r_k = \sum_{i=1}^n \mathbf{1}\{ T_i \ge t_{(k)}\}.
\] and the number of individuals at risk just prior to \(t_{(k)}\) excluding exact observations at \(t_{(k)}\) as \[
r^\star_k = r_k - d_k.
\] We recall that the Kaplan–Meier estimator of the survival function \(S(t)\) is defined as \[
\widehat{S}_{KM}(t)
= \prod_{k \, : \, t_{(k)} \le t}
\left(1 - \frac{d_k}{r_k}\right),
\qquad t \ge 0,
\] with the convention that an empty product equals 1.
Let \(t_f>t_{(k_n)}\) be a final time up to which we conduct our survival analysis, modeling with a \(\text{NTR}(Y)\) prior where the subordinator \(Y\) has intensity \(\nu(\mathrm{d}w,\mathrm{d}u)\), we denote \(\eta_t(\mathrm{d}w)=\nu(\mathrm{d}w,(0,t])\) If \[
\eta_{t_0}'( \mathrm{d}w)=\frac{\partial \eta_t(\mathrm{d}w)}{\partial t}\Bigg|_{t=t_0}\] exists for all \(t_0\in(0,t_f]\) then, restricted to \([0,t_f]\), the posterior distribution given survival data \(\pmb{\mathcal{D}}\) is \(\text{NTR}(\tilde{Y}_{\pmb{\mathcal{D}}})\) with
and \(\left\{ W_k \right\}_{k=1}^{k_n}\perp Y^\circ\). This result is known as the posterior chracterization for NTR models and is akin to conjugacy in Bayesian parametric models not only conceptually but also as a property that provides computational tractability.
Under the simplifyng assumption of \(\nu\) being homogeneous and such that
Bernstein-von Mises theorem and its asymtotic consequences
Under technical conditions on the intensity \(\nu(\mathrm{d}w,\mathrm{d}u)\) underlying a NTR model a Bernstein-von Mises type result holds. This means that in the asymptotic regime where the sample size of survival data tends to infinity then Bayesian credible sets derived from the NTR posterior coincide with frequentist confidence statements based on the Kaplan–Meier estimator. More precisely, let \(S_0\) denote the true survival function of the event times \(\{X_i\}_{i=1}^n\), which we assume to be continuous, and let \(Q_0\) denote the survival function of the censoring times \(\{C_i\}_{i=1}^n\). That is \[S_0(t)=\mathbb{P}\left[ X_1 > t\right], \quad Q_0(t)=\mathbb{P}\left[ C_1 > t\right].\] Let \(\kappa_0\) be the cumulative hazard function associated to \(S_0\), define \[
U_0(t) = \int_0^t \frac{ 1 }{ S_0(w)Q_0(w) } \mathrm{d}\kappa_0(w),
\] and let \(B\) be a standard brownian motion in \(\mathbb{R}\), then
\[
\mathcal{L}\left( \sqrt{n} \left( S(\cdot) - \widehat{S}_{KM}^{(n)}(\cdot) \right) \, \big| \, \pmb{\mathcal{D}}_n \right) \stackrel{d}{\to } S_0(\cdot) B\left( U_0(\cdot) \right)
\] in the space of right continuous with left limit functions in \([0,t_f)\), endowed with the Skorokhod topology, almost surely with respect to the infinite product probability measure \(P^\infty_{1-S_0,1-Q_0}\), which generates the infinite sequence of survival data \(\pmb{\mathcal{D}}_\infty\).
So frequentist and Bayesian nonparametric inferences are equivalent in the big sample size limit as it is known that
and observed approximations based on the survival data are used for inference as
\[
S_0(\cdot) \approx \widehat{S}_{KM}^{(n)}(\cdot) \exp\left( -n^{-\frac{1}{2}} B\left( \hat{U}_n(\cdot) \right) \right),
\] where \[
\hat{U}_n(t) = n\sum_{i\,:\,T_i \leq t}^n \frac{d_i}{r_i^2},
,
\] more detail on the frequentist approach derivations is provided in the following section for the Cox model.
The NTR model with a Gamma process subordinator can be shown to satisfy such a Bernstein–von Mises result. Recalling that the Kaplan–Meier estimator is consistent for continuous \(S_0\), it follows that the posterior mean and posterior median survival estimators are also consistent. This consistency occurs at the optimal rate, since the \(\sqrt{n}\) scaling coincides with the convergence rate for the asymptotic Gaussian distribution of the Kaplan–Meier estimator, which results from a martingale central limit theorem,
For survival datasets of large size, frequentist uncertainty quantification based on the asymptotic Gaussian process distribution yields confidence bands for the true survival function that are asymptotically equivalent to credible bands drawn from the NTR posterior, under technical conditions. The main frequentist confidence bands in the literature are Hall–Wellner, although alternatives such as the equal-precision, also known as Nair, bands exist. For a full treatment of frequentist confidence bands see Klein & Moeschberger 2008. Alternatively, confidence bands can be constructed by bootstrapping the Kaplan–Meier estimator,(Akritas, 1986), yielding bands that are asymptotically equivalent to those derived from the martingale central limit theorem. Later on, we will visualize the Bernstein-von Mises result for our Gamma process workflow.
Likelihood computation for NTR models
Computation of the likelihood for NTR models with survival data is is provided by
In the particular case that there are no repetitions in the exact observations, i.e. \(d_k=1\) for \(1\leq k \leq k_n\), then no use of a binomial expansion is necessary and we obtain the simplification
From the posterior characterization we can obtain the posteror mean survival function, which can be used as an estimator for the survival function of the event times of interest,
where \(t_{(k_n+1)}=\infty\) and \(\psi^\circ_t(\lambda)= -\log\left( \mathbb{E}\left[ e^{-\lambda Y^\circ(t)} \right] \right)\) is the Laplace exponent of \(Y^\circ(t)\).
Again, under the case of \(\nu\) being homogeneous and such that $( w,(0,t]) = (t) (w) w $ we have the simplification
We simulate 333, independent, Weibull survival times without censoring. Using the prior specification described above, we fit a NTR model and compare the Kaplan–Meier estimator with the posterior mean survival curve. This comparison illustrates the asymptotic agreement between the two procedures.
Show Julia code
# Implementation of Weibull simulation study and posterior mean illustration# Synthetic data generationn =333λ=0.7# Weibull rate parameterk =2.5# Weibull shape parameterT =rand( Weibull(k,λ), n)δ =zeros(Int64,n) .+1# no censoring# Compute Kaplan-Meier estimator via Survival packagekm =fit(KaplanMeier, T, δ)# NTR survival data structdata =SurvivalData(T,δ)# NTR model structmodel =NeutralToTheRightModel( α, baseline, data)# Plotting time grid updated for better visualizationt =collect(LinRange(0.0,1.75,100))# Posterior mean survival computationSm =mean_posterior_survival(t,model)# Posterior computations plotpost_mean_plot =plot( t, Sm, xlabel="\$t\$", ylabel="\$S(t)\$", linewidth=2, label="Posterior mean (analytical)",left_margin =5.5mm, bottom_margin =5mm)plot!( t, ccdf.( Weibull(k,λ), t), color=:black, linestyle=:dot, linewidth=2, label="True survival \$S_0\$")plot!(km.events.time,km.survival,seriestype =:steppost, linestyle=:dash, c=3, linewidth=2,label="Kaplan-Meier estimator")
Posterior simulation of NTR survival functions
From the posterior characterization of the NTR model, we have that simulation of posterior survival draws is feasible by drawing
It is seen from the posterior characterization that in order to draw simulations of \(\tilde{Y}_{\pmb{\mathcal{D}}}\) in \([0,t_f]\) we only need to:
Draw Gamma processes in \((t_{(k-1)},t_{(k)}]\) with parameters \(\kappa\), \(\beta_\alpha\) and \(\alpha^\circ=\alpha+r_k\).
Draw random weights, corresponding to fixed jump location \(t_{(k)}\), \(W_k\sim\mathcal{L}(f_k)\).
For the random weights in our Gamma process workflow we use an efficient acceptance-rejection sampler in the sense that the probability of acceptance goes to one as \(r_k^\star\) goes to infinity, provided that \(\frac{1}{n}\sum_{i=1}^n \delta_i\) converges to a positive limit. Complexity of the draw for random weights depends on the associated jump location \(t_{(k)}\) corresponds to repeated exact observations on the survival data or not. Details on the acceptance-rejection scheme follow next.
No repetitions case
Let the kernel of \(f_k\) be denoted by \[
f_k^{\text{ker}}= e^{-\left( \alpha + r^\star_k \right)s}\left( 1-e^{-s} \right)s^{-1}
\] Use proposal distribution \(\text{Exp}(\alpha+r_k^\star)\) and accept a draw \(V\) from it if for an independent \(U\sim\text{Uniform}(0,1)\)\[
\frac{1-e^{-V}}{V}\leq U.
\]
Repetitions case
In this case \[
f_k^{\text{ker}}= e^{-\left( \alpha + r^\star_k \right)s}\left( 1-e^{-s} \right)^{d_k}s^{-1}
\] and we use proposal distribution \(\text{Gamma}(d_k,\alpha + r_k^\star)\) and accept a draw \(V\) from it if for an independent \(U\sim\text{Uniform}(0,1)\)\[
\left( \frac{1-e^{-V}}{V} \right)^{d_k}\leq U.
\]
We draw posterior samples and construct credible bands which illustrate the Bernsetin-von Mises theorem for the NTR model Gamma process workflow.
We can consider the inclussion of covariates in NTR based models. Denote such covariates as \(\pmb{Z}\) for \(T\) a possibly censored to the right random variable. The Cox regression model for NTR Cox-baseline is determined by covariate conditional survival function \[
S(t;\pmb{Z})= \exp\left( -e^{g(\pmb{c},\pmb{Z})}Y(t) \right)
\]
We will draw covariate dependent survival data from a Weibull Cox regression model.
We observe how the Cox covariate dependence structure drives survival data away from the baseline distribution.
Show Julia code
# Implementation of Weibull simulation study with covariates and histogram illustration# Covariate distributionDz =MixtureModel(Normal[Normal(-2.5, 1.2),Normal(2.5, 1.2),Normal(0.0, 1.0)], [1/3, 1/3, 1/3])# Quantiles of covariate distributionq_v =quantile(Dz,[1/3,2/3,])# Synthetic data generationm =999Z = [ [rand(Dz)] for _ in1:m ]randXgivenZ(z::Vector{Float64};b::Float64,k::Float64,λ::Float64) =rand(Weibull(k, λ /exp(b * z[1] / k)))T = [ randXgivenZ(z;b,k,λ) for z in Z ]δ =zeros(Int64,m) .+1histogram(T,normalize=true,label="Data with covariate dependence", xlabel="\$t\$", ylabel="\$f(t)\$")plot!( collect(minimum(T):0.01:maximum(T)), pdf.(Weibull(k,λ),collect(minimum(T):0.01:maximum(T))), label="Cox-baseline density \n (no covarate dependence)" )
Show Julia code
# Covariate histogramhistogram([ z[1] for z in Z], c=3, normalize=true)plot!( collect(-6:0.1:6), pdf.(Dz,collect(-6:0.1:6)),c=4, linewidth=4 )
Posterior characterization (conjugacy structure) for Cox NTR priors
For the Cox regression setting we define \[\begin{align*}
h^{(e)}_k(\pmb{c}) &= \sum_{\{i\,:\, t_i = t_{(k)},\, \delta_i=1 \}} g(\pmb{c},\pmb{Z}_i), \quad \bar{h}^{(c)}_k(\pmb{c}) = \sum_{ \{i\,:\, t_i = t_{(k)},\, \delta_i=0 \} }g(\pmb{c},\pmb{Z}_i)
\\
\bar{h}^{(e)}_k(\pmb{c}) &= \sum_{i=k}^{k_n} h^{(e)}_i(\pmb{c},\widehat{\pmb{Z}}), \quad \hphantom{} \bar{h}^{(c)}_k(\pmb{c}) = \sum_{i=k}^{k_n} h^{(c)}_i(\pmb{c},\widehat{\pmb{Z}}).
\end{align*}\] So the quantities related to the numbers of at risk patients become \[\begin{align*}
h_k(\pmb{c})=\bar{h}^{(c)}_k(\pmb{c})+\bar{h}^{(e)}_k(\pmb{c}), &\quad h^\star_k(\pmb{c})=\bar{h}^{(c)}_k(\pmb{c})+\bar{h}^{(e)}_{k+1}(\pmb{c})
\end{align*}\] where \(\bar{h}^{(e)}_{k_n+1}(\pmb{c})=0\).
The posterior distribution of the Cox-baseline \(Y\) is determined by
Consider a posteriori the survival time \(X^{(\text{new})}\) corresponding to a new patient with associated covariate vector \(\pmb{Z}^{(\text{new})}\). The posterior distribution of \(g(\pmb{c},\pmb{Z}^{(\text{new})})Y\) in \((0,t_f]\) is characterized by
the random jump weights \(W_k\), corresponding to a jump at fixed location \(t_{(k)}\), have probability density function, with respect to dominating measure \(\eta'_{t_{(k)}}\circ m ^{-1}\),
Under suitable regularity conditions on the prior and the baseline hazard, it is satisfied that \[
\mathcal{L}\!\left( \sqrt{n}\left( S(\cdot) - \widehat{S}_{BR}^{(n)}(\cdot) \right),\;
\sqrt{n}\left( \mathbf{c} - \hat{\pmb{c}}_{ML}^{(n)} \right)
\;\middle|\; \mathbf{\mathcal{D}}_n
\right)
\;\xrightarrow{d}\; \left( S_0(\cdot) \left( B\left( U_0(\,\cdot\,;\pmb{c}_0) \right) \right) - \pmb{V}_0'(\, \cdot \, ; \pmb{c}_0)\pmb{I}(\mathbf{c}_0)^{-\frac{1}{2}}\pmb{W},\,
\pmb{I}(\mathbf{c}_0)^{-\frac{1}{2}}\pmb{W} \right)
\] almost surely under the true data-generating law. So frequentist and Bayesian nonparametric inferences are again equivalent for the Cox model in the big sample size limit as it is also known that
In the frequentist approach observed aproximations are used as \[
\sqrt n\left( \hat{\kappa}_{BR}^{(n)}(\cdot)-\kappa_0(\cdot)\right)
\;\xrightarrow{d}\;
\hat{\pmb{V}}_n'( \cdot )\hat{\pmb{I}}_n^{-\frac{1}{2}}\pmb{W}
- B\left( \hat{U}_n(\cdot)
\right)
\] and, obtained from application of the delta method in the aproximation above, \[
\sqrt n\left( \hat{S}_{BR}^{(n)}(\cdot)-S_0(\cdot)\right)
\;\xrightarrow{d}\;
\hat{S}_{BR}^{(n)}(\cdot) \left( \pmb{V}_n'( \cdot)\hat{\pmb{I}}_n^{-\frac{1}{2}}\pmb{W}
- B\left( \hat{U}_n(\cdot) \right)
\right).
\] So for frequentist uncertainty quantification that naturally compares with the Bayesian NTR approach we can sample the asymptotic survival curves to construct via Monte-Carlo confidence bands by drawing from the Gaussian approximation on the log-survival scale and then exponentiating, which keeps the simulated curves in the survival scale. \[
\widehat{S}_{BR}^{(n)}(\cdot) \exp\left( n^{-\frac{1}{2}}\left( \hat{\pmb{V}}_n(\cdot) \hat{\pmb{I}}_n^{-\frac{1}{2}}\pmb{W} - B\left( \hat{U}_n(\cdot) \right) \right) \right)
\] For a new subject with corresponding covariate vector \(\pmb{Z}^{(\text{new})}\) we estimate its survival curve \(S(t;\pmb{Z}^{(\text{new})} )=S_0(t)^{\exp\left( \pmb{c}'\pmb{Z}^{(\text{new})} \right) }\) using \(\hat{S}_{BR}^{(n)}(t)=\exp\left( -\hat{\kappa}_{BR}(t) \right)\) so a nataural estimator for the subject-specific survival is \[
\hat{S}_{BR}^{(n)}(t; \pmb{Z}^{(\text{new})}) = \exp\left( -e^{ {\hat{\pmb{c}}_{ML}^{(n)}}' \pmb{Z}^{(\text{new})} } \hat{\kappa}_{BR}(t) \right)
\] Aplying the delta method to the log-survival function in the Cox model given by \[
\gamma_t(\kappa,\pmb{c}) = \log\left( S(t;\pmb{z},\kappa,\pmb{c}) \right) = -\exp\left( \pmb{c}'\pmb{z} \right)\kappa(t)
\] so \[
\frac{\partial \gamma_t}{\partial \kappa}
\bigg|_{(\hat \kappa,\hat{\pmb c})}
=
-\exp\left( \hat{\pmb{c}}'\pmb{z} \right),
\]
and
\[
\nabla_{\pmb c}\gamma_t(\kappa,\pmb c)
\bigg|_{(\hat \kappa,\hat{\pmb c})}
=
-\hat \kappa(t) \exp\left( \hat{\pmb{c}}'\pmb{z} \right) \pmb z'.
\] We apply the delta method; informally let \[
\sqrt{n}\Delta \kappa(\cdot) = \sqrt{n}\left( \kappa_0(\cdot) - \hat{\kappa}_{BR}^{(n)}(\cdot) \right) \approx
B\left( \hat{U}_n(\cdot) \right) - \hat{\pmb{V}}_n'( \cdot )\hat{\pmb{I}}_n^{-\frac{1}{2}}\pmb{W}
\] and
\[
\sqrt{n}\Delta \pmb{c} = \sqrt{n}\left( \pmb{c}_0 - \hat{\pmb{c}}_{ML}^{(n)} \right) \approx \hat{\pmb{I}}_n^{-\frac{1}{2}}\pmb{W}
\] so from the Taylor expansion
Likelihood computation with Gamma processes workflow for survival analysis
Denoting \(l_k^\star = \underset{\{ l\,:\, t_l=t_{(k)},\, \delta_l = 1 \}}{\arg\min} g\left( \pmb{c},\pmb{Z}_l \right)\), doing a product expasion and using the Laplace exponent formula for Gamma processes we obtain
In the particular case that there are no repetitions in the exact observations, i.e. \(d_k=1\) for \(1\leq k \leq k_n\), then no use of a product expansion is necessary and we obtain the simplification
We use the above likelihood to compute the posterior density over the regression coefficient and compare with both the frequentist partial likelihood scaled to be a density and the asymptotic Gaussian distribution whch appears also in the Bernstein-von Mises theorem.
Show Julia code
# Cox-NTR model setup and illustration of posterior distribution in regression coefficent# DataFrame for Survival.jl fitdf =DataFrame(time = T, status = δ, z = [ z[1] for z in Z ])# Add EvenTime column to DataFramedf.event =EventTime.( df.time, df.status)# Fit fequentist Cox modelfreq_cox_model = Survival.coxph(@formula(event ~ z), df)# Frequentist fit of Cox regression coefficient, c_freq =coef(freq_cox_model)# Regression data struct for Cox-NTR modelingdataregre =RegressionSurvivalData(T,δ,Z)# Asymptotic vraince components computationfisherI_obs, V_obs, U_obs =cox_var_obs(t,c_freq[1],dataregre)σ_freq_obs =1/sqrt( m*fisherI_obs )# Cox-NTR model variance modualting parameterα =5.0# Compile logLlikelihood evaluationNTRsurv.loglikelihood([0.3],α,baseline,dataregre) # In this notebbok loglikelihood funcion is also loaded by Survival so it has to be sepcified.# Illustrative plot of proportional density on bpost_k(v) =exp(NTRsurv.loglikelihood( v, α, baseline, dataregre) - NTRsurv.loglikelihood([b],α,baseline,dataregre))post_C =hcubature( post_k, [-1.5], [2.0] )[1]post(z) =post_k(z)/post_Ccoxpll_k(v) =exp( cox_partial_loglik(v,dataregre) -cox_partial_loglik(b,dataregre) )freq_C =hcubature( coxpll_k, [-1.5], [2.0] )[1]coxpll(z) =coxpll_k(z)/freq_C# Plot of posterior density on regression coefficient compared with density scaled partial loglikelihood and true coefficient `b=0.3`plot(collect(-0.5:0.001:0.75), [ post([x]) for x incollect(-0.5:0.001:0.75) ], xlabel="\$b\$", ylabel="\$\\,f\\,(b)\$", label="NTR posterior", xlim=(0.0,0.75) )plot!(collect(-0.5:0.001:0.75), [ coxpll([x]) for x incollect(-0.5:0.001:0.75) ], label="Cox partial likelihood (density scaled)", xlim=(0.0,0.75) )plot!(collect(-0.5:0.001:0.75), pdf.(Normal(c_freq[1],σ_freq_obs), collect(-0.5:0.001:0.75)), linestyle=:dash,label="Frequentist asymptotic")vline!([b],label="True b")
Posterior mean survival for Cox-NTR models
Let \(X^{(\text{new})}\) be a new observation with corresponding covariate vector \(\pmb{Z}^{(\text{new})}\). From the posterior characterization of the Cox-NTR model we can obtain a posteror mean estimator for the survival function of \(X^{(\text{new})} \,|\,\pmb{Z}^{(\text{new})}\)
where \(t_{(k_n+1)}=\infty\) and \(\psi^{\circ \pmb{Z}^{(\text{new})}}_t(\lambda)= -\log\left( \mathbb{E}\left[ \exp\left( -\lambda Y_{\pmb{c},\pmb{Z}^{(\text{new})}}^\circ(t) \right)\right] \right)\) is the Laplace exponent of \(Y_{\pmb{c},\pmb{Z}}^\circ(t)\).
Again, under the case of \(\nu\) being homogeneous and such that $(w,(0,t]) = (t)(w) $ we have the simplification
and for the Gamma process case in particular we obtain, denoting again \(l_k^\star = \underset{\{ l\,:\, t_l=t_{(k)},\, \delta_l = 1 \}}{\arg\min} g\left( \pmb{b},\pmb{Z}_l \right)\), doing a product expasion and using the corresponding Laplace exponent formula,
We use Metropolis-Hastings, with Robbins-Monro stochastic approximation algorithm for tuning of proposal variance, for marginal simulation of the regression coefficient \(\pmb{c}\). We can evaluate the posterior mean estimator above using a point estimate for \(\pmb{c}\). In this example we use the posterior mean obtained from the Markov chain.
Show Julia code
# MCMC run# Compilation runrobbins_monro_mh_tune(3,3,x->NTRsurv.loglikelihood([x],α,baseline,dataregre),1.0,1.0,0.4,0.7)# Robbins-Monro algorithm for variance of proposal distribution tuningsd_prop_tuned, s₀_tune, lliks₀_tune =robbins_monro_mh_tune(500,250,x->NTRsurv.loglikelihood([x],α,baseline,dataregre),1.0,1.0,0.6,0.7)# Metropolis-Hastings chain runchain_s, _ =random_walk_mh( 3000, x->NTRsurv.loglikelihood([x],α,baseline,dataregre), s₀_tune, lliks₀_tune, sd_prop_tuned[end])# Acceptance rate in the MH chain, target value was 0.6acceptance_rate(chain_s)# Plot MH chainplot(chain_s,label="Tune M-H chain")
# Illustration of posterior mean comparisson with frequentist Breslow estimator survival curves# Cox-NTR model setup with MCMC posterior mean regression coefficient and illustration of posterior mean survival# Posterior mean estimate of regression coefficientc_post =mean(chain_s)# NTR model for Cox regressionplugin_cox_model =PluginCoxNeutralToTheRightModel( [c_post], α, baseline, dataregre)# Breslow estimator for cumulative hazardHfreq =BreslowEstimCoxModel(freq_cox_model)# Frequentists fits of survival curveszv = [ 0.0, 2.0, -2.0 ]gv =exp.( c_freq[1].*zv )Sfreq0 = [1.0; exp.( - gv[1]*Hfreq )]Sfreq1 = [1.0; exp.( - gv[2]*Hfreq )]Sfreq2 = [1.0; exp.( - gv[3]*Hfreq )]# Plot posterior mean survival curves in Cox-NTR modelSm0 =mean_posterior_survival(t, [0.0], plugin_cox_model);post_mean_regre_plot =plot(t,Sm0,c=1,xlabel="\$t\$", ylabel="\$S(t)\$",label="Plugin NTR post. mean surv. Z=0")plot!(post_mean_regre_plot, t, ccdf.( Weibull( k, λ), t), c=:blue, linestyle=:dash, label="True surv. Z=0" ) plot!( post_mean_regre_plot,[0.0;dataregre.T], Sfreq0, seriestype =:steppost, c=:cyan, linestyle=:dot,alpha=0.75,label="Frequentist Cox estimator Z=0")Sm1 =mean_posterior_survival(t,[2.0],plugin_cox_model)plot!(post_mean_regre_plot ,t, Sm1,c=2,label="Plugin NTR post. mean surv. Z=2")plot!(post_mean_regre_plot, t, ccdf.( Weibull( k, λ/exp(b*2.0/k)), t), c=:red, linestyle=:dash, label="True surv. Z=2" ) plot!( post_mean_regre_plot,[0.0;dataregre.T], Sfreq1, seriestype =:steppost, c=:magenta, linestyle=:dot, alpha=0.75,label="Frequentist Cox estimator Z=2")Sm2 =mean_posterior_survival(t,[-2.0],plugin_cox_model)plot!( post_mean_regre_plot,t,Sm2,c=3,label="Plugin NTR post. mean surv. Z=-2")plot!( post_mean_regre_plot, t, ccdf.( Weibull( k, λ/exp(-b*2.0/k)), t), c=:green, linestyle=:dash, label="True surv. Z=-2" ) plot!( post_mean_regre_plot, [0.0;dataregre.T], Sfreq2, seriestype =:steppost, c=:lime, linestyle=:dot, alpha=0.75,label="Frequentist Cox estimator Z=-2")
For comparisson we fitted the frequentist Cox model using Survival.jl and plotted the survival curves based on the Breslow estimator for the baseline cumulative hazard function of the frequentist Cox model.
Posterior simulation of NTR survival functions
From the posterior characterization of the NTR model, we have that simulation of posterior survival draws is feasible by drawing
It is seen from the posterior characterization that in order to draw simulations of \(\tilde{Y}_{\pmb{c},\pmb{Z}^{(\text{new})},\pmb{\mathcal{D}}}\) in \([0,t_f]\) we only need to:
Draw Gamma processes in \((t_{(k-1)},t_{(k)}]\) with parameters \(\kappa\), \(\beta_\alpha\) and \(\alpha^\circ= (\alpha+h_k)/g\left(\pmb{c},\pmb{Z}^{(\text{new})} \right)\).
Draw random weights, corresponding to fixed jump location \(t_{(k)}\), \(W_k\sim\mathcal{L}(f_{k,\pmb{Z}^{(new)}})\).
For the sampling the random weights in the posterior characterization we use an efficient acceptance-rejection sampler in the sense that the probability of acceptance goes to one as \(h_k^\star(\pmb{c})\) goes to infinity, provided that \(\frac{1}{n}\sum_{i=1}^n \delta_i\) converges to a strictly positive limit. Complexity of the draw for random weights depends on the associated jump location \(t_{(k)}\) corresponds to repeated exact observations on the survival data or not. Details on the acceptance-rejection scheme follow next.
No repetitions case
Let the kernel of \(f_{k,\pmb{Z}^{(\text{new})}}\) be denoted by \[
f_{k,\pmb{Z}^{(\text{new})}}^{\text{ker}}= e^{-\left( \alpha + r^\star_k \right)s}\left( 1-e^{-s} \right)^{d_k}s^{-1}
\] Use proposal distribution \(\text{Exp}(\alpha+r_k^\star)\) and accept a draw \(V\) from it if for an independent \(U\sim\text{Uniform}(0,1)\)\[
\frac{1-e^{-V}}{V}\leq U.
\]
Alternatively, a fully Bayesian model can be considered where instead of plugging an estimator of the regression coefficient, we treat it as random with samples coming from the Markov chain that was computed above. We ccompare the posterior mean survivals and credible bands between the non-fully and fully Bayesian models to see that there is not much difference.
Show Julia code
# Computation and posterior mean illustration of fully Bayesian Cox-NTR modelcox_model =CoxNeutralToTheRightModel([ [v] for v in chain_s], α, baseline, dataregre)Sfbm0 =mean_posterior_survival(t, [0.0], cox_model)Sfbm1 =mean_posterior_survival(t,[2.0], cox_model)Sfbm2 =mean_posterior_survival(t,[-2.0], cox_model)compare_post_regre_plot =plot(t,Sm0,c=1, linewidth=1.5, xlabel="\$t\$", ylabel="\$S(t)\$",label="NTR post. mean surv. Z=0")plot!( compare_post_regre_plot, t , Sfbm0, c=:black, linestyle=:dash, label="NTR post. mean surv. Z=0")plot!( compare_post_regre_plot, t, Sm1, c=2, linewidth=1.5, label="NTR post. mean surv. Z=2")plot!( compare_post_regre_plot, t , Sfbm1, c=:yellow, linestyle=:dash, label="NTR post. mean surv. Z=2")plot!( compare_post_regre_plot, t,Sm2,c=3,linewidth=1.5, label="NTR post. mean surv. Z=-2")plot!( compare_post_regre_plot, t, Sfbm2, c=:purple, linestyle=:dash, label="NTR post. mean surv. Z=-2")
Show Julia code
# Computation and illustration of fully Bayesian Cox-NTR model credible bandscred0 =posterior_credible_band( t, [0.0], cox_model)cred1 =posterior_credible_band( t, [2.0], cox_model)cred2 =posterior_credible_band( t, [-2.0], cox_model)post_fully_bayes_regre_draws_plot =plot( plugin_cred0, xlabel="\$t\$", ylabel="\$S(t)\$", fillalpha=0.3, label="Plugin NTR post. mean surv. \nwith 95% credible band Z=0",size=(600,400))plot!(post_fully_bayes_regre_draws_plot, plugin_cred1, c=2,label="Plugin NTR post. mean surv. \nwith 95% credible band Z=2")plot!(post_fully_bayes_regre_draws_plot, plugin_cred2, c=3,label="Plugin NTR post. mean surv. \nwith 95% credible band Z=-2")plot!(post_fully_bayes_regre_draws_plot, t, cred0.d, c=:black, linestyle=:dash, label="NTR post. mean surv. \n95% credible band Z=0")plot!(post_fully_bayes_regre_draws_plot, t, cred0.u, c=:black, linestyle=:dash, label="")plot!(post_fully_bayes_regre_draws_plot, t, cred1.d, c=:yellow, linestyle=:dash, label="NTR post. mean surv. \n95% credible band Z=2")plot!(post_fully_bayes_regre_draws_plot, t, cred1.u, c=:yellow, linestyle=:dash, label="")plot!(post_fully_bayes_regre_draws_plot, t, cred2.d, c=:purple, linestyle=:dash, label="NTR post. mean surv. \n95% credible band Z=-2")plot!(post_fully_bayes_regre_draws_plot, t, cred2.u, c=:purple, linestyle=:dash, label="",title="",size=(600,400))
The similarity between the fully Bayesian Cox-NTR survival estimates and the plug-in Cox-NTR approach should not be surprising as the Bernstein-von Mises theorem implies that for relatively big sample the posterior density of the regression coefficients accumulates near the truth. However also in small sampe size regimes, since computational complexity of the fully Bayes approach requires computing Cox risk-set sufficient statistics across posterior draws of the regression coefficient, the plug-in approach is to be preferred as it provides a substantially cheaper default workflow. We therefore recommend the plug-in Cox-NTR analysis for routine use, while retaining the fully Bayesian version as an optional sensitivity analysis for small samples, weakly identified regression effects, or predictions at extreme covariate values.
