using Pkg, NTRsurv
pkg_dir = pkgdir(NTRsurv)
Pkg.activate(joinpath(pkg_dir, "examples"))
Pkg.instantiate()
# Load packages needed for tutorial
using Random, CSV, DataFrames, Distributions, Survival, NTRsurv, Plots, LaTeXStrings
# Set seed for reproducibility
Random.seed!(1234)This tutorial shows how to set informative and non-informative priors in NTR models and their impact for inference with small samples using a real data example with the Julia package NTRsurv.jl. For theoretical details of NTR modelling see the Weibull simulation study tutorial.
Tongue tumor data, without the type of DNA profile covariate, from Klein and Moeschberger (1997) Survival Analysis Techniques for Censored and truncated data. Obtained from R package survival.
# Load data
tongue_path = joinpath( dirname(dirname(pathof(NTRsurv))), "examples", "data", "tongue.csv")
tongue = CSV.read(tongue_path, DataFrame)
first(tongue,7)| Row | time | status |
|---|---|---|
| Float64 | Int64 | |
| 1 | 1.0 | 1 |
| 2 | 3.0 | 1 |
| 3 | 3.0 | 1 |
| 4 | 4.0 | 1 |
| 5 | 10.0 | 1 |
| 6 | 13.0 | 1 |
| 7 | 13.0 | 1 |
# Observed times
tongue_T = float.(tongue[!,:time])
# Observed status
tongue_δ = tongue[!,:status];# Number of observations
length(tongue_T)80# Exact observations
tongue_Tₑ = tongue_T[ tongue_δ .== 1 ]
# Censored observations
tongue_Tₒ = tongue_T[ tongue_δ .== 0 ];# Number of exact observations
sum(tongue_δ)53# Kaplan-Meier fit via Survival.jl
km = fit(KaplanMeier, tongue_T, tongue_δ)KaplanMeier{Float64, Float64}(EventTable{Float64}([1.0, 3.0, 4.0, 5.0, 8.0, 10.0, 12.0, 13.0, 16.0, 18.0 … 129.0, 131.0, 150.0, 157.0, 167.0, 176.0, 181.0, 231.0, 240.0, 400.0], [2, 3, 2, 2, 1, 1, 1, 3, 2, 1 … 1, 0, 0, 1, 1, 0, 1, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 0, 0, 0 … 0, 1, 1, 0, 0, 1, 0, 2, 1, 1], [80, 78, 75, 73, 71, 69, 68, 67, 64, 62 … 11, 10, 9, 8, 7, 6, 5, 4, 2, 1]), [0.975, 0.9375, 0.9125000000000001, 0.8875000000000001, 0.8750000000000001, 0.8623188405797103, 0.8496376811594205, 0.8115942028985509, 0.7862318840579712, 0.7735507246376814 … 0.2883605848313903, 0.2883605848313903, 0.2883605848313903, 0.2523155117274665, 0.21627043862354273, 0.21627043862354273, 0.1730163508988342, 0.1730163508988342, 0.1730163508988342, 0.1730163508988342], [0.01790287185098582, 0.028867513459481287, 0.034621238235890715, 0.03980586693257064, 0.042257712736425826, 0.04470842219836276, 0.04709919104246099, 0.05401815123202459, 0.058497802698822765, 0.06071575813144969 … 0.21297949710464092, 0.21297949710464092, 0.21297949710464092, 0.25143072414501894, 0.2950032760048817, 0.2950032760048817, 0.37017149114108233, 0.37017149114108233, 0.37017149114108233, 0.37017149114108233])# Kaplan-Meier fit illustration
plot(km.events.time,km.survival, seriestype = :steppost, ylim=(-0.1,1), xlabel = "Time", ylabel = "S(t)", label = "", title="Kaplan-Meier estimator")
scatter!( tongue_Tₑ, zeros(sum(tongue_δ)) .- 0.025, markerstrokewidth=0.0, color = 1, alpha=0.5, label="Exact. obs")
scatter!( tongue_Tₒ, zeros(length(tongue_δ)-sum(tongue_δ)) .- 0.05, marker=:cross, color = :black, alpha=0.75, label="Cens. obs.")
We set an informative prior using an empirical Bayesian approach were the baseline is the parametric exponential model survival with maximum-likelihood fit for the survival data.
# Create NTRsurv data structure
tongue_data = SurvivalData(tongue_T,tongue_δ)Right-censored survival data
Observations: 80
Exact events: 53 (66.2%)
Right-censored: 27 (33.8%)
Distinct times: 59 (73.8%)# Create NTRsurv baseline structure
info_baseline = EmpiricalBayesBaseline(tongue_data)Neutral-to-the-right baseline
Specification: Empirical ExponentialBaseline(rate=0.0089739)# Informative baseline mean
println("The baseline mean is ", 1/info_baseline.κ(1.0))The baseline mean is 111.43396226415095# Informative prior beliefs illustration
# Plotting time grid
t = collect(LinRange(0.0,400,200))
# Prior credible band
prior_band_info = prior_credible_band(t,1.0,info_baseline)Credible band
Credible level: 95.0%
Estimate: Prior mean
Time grid: 0.0 to 400.0 (200 points)
Monte-Carlo paths: 3000
Median band width: 0.89882
Maximum band width: 0.97474# Plot informative prior
info_prior_plt = plot( prior_band_info, fillalpha=0.3, xlabel="\$t\$", ylabel="\$S_0(t)\$", label="Prior mean baseline \nwith 95% band", title="Informative pior", size=(600,400))
# Illustration of NTR fit with informative model
# Informative NTR model structure
model_info = NeutralToTheRightModel( 1.0, info_baseline, tongue_data)Neutral-to-the-right model
Gamma process α: 1.0
Center: Empirical ExponentialBaseline(rate=0.0089739)
Data: SurvivalData(obs.=80, exact=53, censored=27, distinct=59)# Posterior credible band computation
posterior_band_info = posterior_credible_band(t,model_info)Credible band
Credible level: 95.0%
Estimate: Posterior median
Time grid: 0.0 to 400.0 (200 points)
Monte-Carlo paths: 3000
Median band width: 0.28831
Maximum band width: 0.29636print(posterior_band_info)CredibleBand( level=0.050000000000000044, estimate=posterior median, draws=3000)# Plot informative posterior
post_info_plot = plot( posterior_band_info, ylim=(-0.1,1), fillalpha=0.3, label="Posterior NTR mean \n with 95% credible band",size=(600,400))
plot!( post_info_plot, km.events.time,km.survival, seriestype = :steppost, c=3, xlabel = "\$t\$", ylabel = "\$S(t)\$", label = "Kaplan-Meier estimator", title="Posterior under informative prior")
scatter!( post_info_plot, tongue_Tₑ, zeros(sum(tongue_δ)) .- 0.025, markerstrokewidth=0.0, color = 1, alpha=0.5, label="Exact. obs")
scatter!( post_info_plot, tongue_Tₒ, zeros(length(tongue_δ)-sum(tongue_δ)) .- 0.05, marker=:cross, color = :black, alpha=0.75, label="Cens. obs.")
For a flatter tail behaviour we can select a more non-informative baseline.
# Create non-informative NTRsurv exponential baseline
noninfo_baseline = ExponentialBaseline(0.001)Neutral-to-the-right baseline
Specification: ExponentialBaseline(rate=0.001)# Non-nformative baseline mean
println("The baseline mean is ", 1/noninfo_baseline.κ(1.0))The baseline mean is 1000.0# Non-nformative prior beliefs illustration
# Prior credible band
prior_band_noninfo = prior_credible_band(t,1.0,noninfo_baseline)Credible band
Credible level: 95.0%
Estimate: Prior mean
Time grid: 0.0 to 400.0 (200 points)
Monte-Carlo paths: 3000
Median band width: 0.87594
Maximum band width: 0.87867# Plot non-informative prior baseline credible band
plot( prior_band_noninfo, fillalpha=0.3, xlabel="\$t\$", ylabel="\$S_0(t)\$", label="Prior mean baseline \n with 95% band", title="Non-informative pior", size=(600,400))
# Illustration of NTR fit with informative model
# Informative NTR model structure
model_noninfo = NeutralToTheRightModel( 1.0, noninfo_baseline, tongue_data)Neutral-to-the-right model
Gamma process α: 1.0
Center: ExponentialBaseline(rate=0.001)
Data: SurvivalData(obs.=80, exact=53, censored=27, distinct=59)# Posterior credible band computation
posterior_band_noninfo = posterior_credible_band(t,model_noninfo)Credible band
Credible level: 95.0%
Estimate: Posterior median
Time grid: 0.0 to 400.0 (200 points)
Monte-Carlo paths: 3000
Median band width: 0.3
Maximum band width: 0.30494# Plot non-informative posterior
post_noninfo_plot = plot( posterior_band_noninfo, ylim=(-0.1,1), fillalpha=0.3, label="Posterior mean (MC) \n with 95% band",size=(600,400))
plot!( post_noninfo_plot, km.events.time,km.survival, seriestype = :steppost, xlabel = "\$t\$", ylabel = "\$S(t)\$", c=3, label = "Kaplan-Meier estimator", title="Posterior under non-informative prior")
scatter!( post_noninfo_plot, tongue_Tₑ, zeros(sum(tongue_δ)) .- 0.025, markerstrokewidth=0.0, color = 1, alpha=0.5, label="Exact. obs")
scatter!( post_noninfo_plot, tongue_Tₒ, zeros(length(tongue_δ)-sum(tongue_δ)) .- 0.05, marker=:cross, color = :black, alpha=0.75, label="Cens. obs.")