# Tutorial environment setup
using Pkg, NTRsurv
pkg_dir = pkgdir(NTRsurv)
Pkg.activate(joinpath(pkg_dir, "examples"))
Pkg.instantiate()
# Load neutral to the righ survival modeling package package
using Random, Plots, LaTeXStrings, CSV, DataFrames, Distributions, NTRsurv, Survival, Measures
# Set random number generator seed for reproducibility
Random.seed!(1234)This tutorial shows how to perform a meta-analysis with cohort frailties NTR-Cox regression models using a cohort depending Cox regression function with the Julia package NTRsurv.jl. For theoretical details of NTR modelling see the Weibull simulation study tutorial.
In this example we use data from the R package mcsurvdata corresponding to patients with breast cancer in different cohorts. Observations with missing cohort value were removed from the data.
We consider a NTR model with LogNormal frailties in the different cohorts as
\[\begin{align*} \epsilon_j & \stackrel{\text{i.i.d.}}{\sim}\text{Normal}(m,s^2),\quad j\in\{1,\ldots,c\}, \\ X_{i} | Y, Z_i & \stackrel{\text{i.i.d.}}{\sim} \text{NTR}\left( \exp(\epsilon_{Z_i})Y \right),\quad i\in\{1,\ldots,n\}, \end{align*}\]
Denoting \(\Gamma_j=\exp(\epsilon_j)\) and \(\pmb{\Gamma}=\left( \Gamma_1,\ldots, \Gamma_c \right)\) we have that a priori \[ S_i(t) = \mathbb{P}\left[ X_i>t \,\Big| Y,Z_i,\pmb{\Gamma}\right] = \exp\left( -\Gamma_{Z_{ i}} Y(t) \right). \]
As \(Y\) already controls the baseline specification we set \(m=0\) to avoid indetifiability issues and assign \(s\sim \text{Normal(0,1)}\) for the variance specification on \(\pmb{\Gamma}\).
# Load data
data_path = joinpath( dirname(dirname(pathof(NTRsurv))), "examples", "data", "brcamcsurvdata.csv")
data = CSV.read(data_path, DataFrame)
first(data,7)| Row | cohort | status | event |
|---|---|---|---|
| Float64 | Float64 | Float64 | |
| 1 | 1.0 | 0.0 | 84.6333 |
| 2 | 1.0 | 1.0 | 163.1 |
| 3 | 1.0 | 0.0 | 164.933 |
| 4 | 1.0 | 1.0 | 7.8 |
| 5 | 1.0 | 0.0 | 164.333 |
| 6 | 1.0 | 0.0 | 99.5333 |
| 7 | 1.0 | 1.0 | 128.0 |
# Array of data in each cohort
dataArr = [ data[ data[:,1] .== i, :] for i in collect(1.0:1.0:5.0) ]5-element Vector{DataFrame}:
398×3 DataFrame
Row │ cohort status event
│ Float64 Float64 Float64
─────┼────────────────────────────
1 │ 1.0 0.0 84.6333
2 │ 1.0 1.0 163.1
3 │ 1.0 0.0 164.933
4 │ 1.0 1.0 7.8
5 │ 1.0 0.0 164.333
6 │ 1.0 0.0 99.5333
7 │ 1.0 1.0 128.0
8 │ 1.0 0.0 36.5667
9 │ 1.0 1.0 36.2667
10 │ 1.0 1.0 132.033
11 │ 1.0 0.0 163.533
⋮ │ ⋮ ⋮ ⋮
389 │ 1.0 0.0 140.5
390 │ 1.0 1.0 41.3667
391 │ 1.0 0.0 108.933
392 │ 1.0 1.0 42.7
393 │ 1.0 0.0 118.2
394 │ 1.0 1.0 28.5
395 │ 1.0 1.0 22.0333
396 │ 1.0 1.0 39.2
397 │ 1.0 1.0 55.4667
398 │ 1.0 0.0 110.967
377 rows omitted
217×3 DataFrame
Row │ cohort status event
│ Float64 Float64 Float64
─────┼────────────────────────────
1 │ 2.0 1.0 10.8667
2 │ 2.0 0.0 241.6
3 │ 2.0 0.0 179.433
4 │ 2.0 0.0 87.5333
5 │ 2.0 0.0 113.833
6 │ 2.0 0.0 236.633
7 │ 2.0 0.0 228.6
8 │ 2.0 1.0 50.0333
9 │ 2.0 0.0 138.333
10 │ 2.0 0.0 228.8
11 │ 2.0 0.0 103.467
⋮ │ ⋮ ⋮ ⋮
208 │ 2.0 0.0 258.333
209 │ 2.0 0.0 269.333
210 │ 2.0 0.0 189.1
211 │ 2.0 1.0 81.1
212 │ 2.0 0.0 256.5
213 │ 2.0 0.0 242.567
214 │ 2.0 0.0 108.767
215 │ 2.0 0.0 236.7
216 │ 2.0 0.0 33.0
217 │ 2.0 1.0 24.8667
196 rows omitted
585×3 DataFrame
Row │ cohort status event
│ Float64 Float64 Float64
─────┼────────────────────────────
1 │ 3.0 0.0 28.0667
2 │ 3.0 1.0 49.6667
3 │ 3.0 1.0 106.567
4 │ 3.0 0.0 188.367
5 │ 3.0 1.0 44.7333
6 │ 3.0 0.0 11.7
7 │ 3.0 0.0 193.167
8 │ 3.0 0.0 139.3
9 │ 3.0 1.0 105.667
10 │ 3.0 1.0 15.2
11 │ 3.0 0.0 298.033
⋮ │ ⋮ ⋮ ⋮
576 │ 3.0 0.0 159.733
577 │ 3.0 1.0 104.767
578 │ 3.0 0.0 153.3
579 │ 3.0 0.0 195.867
580 │ 3.0 0.0 202.767
581 │ 3.0 0.0 220.233
582 │ 3.0 1.0 164.567
583 │ 3.0 0.0 229.833
584 │ 3.0 1.0 36.6333
585 │ 3.0 1.0 48.4333
564 rows omitted
175×3 DataFrame
Row │ cohort status event
│ Float64 Float64 Float64
─────┼────────────────────────────
1 │ 4.0 0.0 110.1
2 │ 4.0 0.0 123.533
3 │ 4.0 0.0 123.533
4 │ 4.0 0.0 122.4
5 │ 4.0 0.0 125.833
6 │ 4.0 0.0 113.367
7 │ 4.0 0.0 124.0
8 │ 4.0 1.0 65.1667
9 │ 4.0 0.0 163.167
10 │ 4.0 0.0 156.333
11 │ 4.0 0.0 124.267
⋮ │ ⋮ ⋮ ⋮
166 │ 4.0 0.0 168.6
167 │ 4.0 1.0 46.6667
168 │ 4.0 1.0 30.1667
169 │ 4.0 0.0 208.967
170 │ 4.0 1.0 109.0
171 │ 4.0 0.0 86.4
172 │ 4.0 0.0 158.633
173 │ 4.0 0.0 207.467
174 │ 4.0 1.0 175.967
175 │ 4.0 0.0 201.9
154 rows omitted
131×3 DataFrame
Row │ cohort status event
│ Float64 Float64 Float64
─────┼────────────────────────────
1 │ 5.0 0.0 135.3
2 │ 5.0 0.0 63.5667
3 │ 5.0 0.0 70.5333
4 │ 5.0 1.0 69.1
5 │ 5.0 0.0 218.233
6 │ 5.0 1.0 211.133
7 │ 5.0 0.0 144.433
8 │ 5.0 0.0 191.933
9 │ 5.0 0.0 143.533
10 │ 5.0 0.0 70.9
11 │ 5.0 1.0 143.167
⋮ │ ⋮ ⋮ ⋮
122 │ 5.0 0.0 200.333
123 │ 5.0 0.0 203.533
124 │ 5.0 0.0 122.8
125 │ 5.0 0.0 202.233
126 │ 5.0 1.0 58.4667
127 │ 5.0 0.0 213.367
128 │ 5.0 0.0 2.3
129 │ 5.0 1.0 30.0667
130 │ 5.0 0.0 264.6
131 │ 5.0 0.0 29.2333
110 rows omitted
# Data sample sizes by cohort
c_sizes = size.(dataArr)5-element Vector{Tuple{Int64, Int64}}:
(398, 3)
(217, 3)
(585, 3)
(175, 3)
(131, 3)c_exfracs = [ sum( dataArr[i].status )/c_sizes[i][1] for i in 1:5]5-element Vector{Float64}:
0.2562814070351759
0.2903225806451613
0.37435897435897436
0.21142857142857144
0.24427480916030533# Kaplan-Meier fits in each cohort
kmArr = [ fit(KaplanMeier, dataArr[i][:,3], Int.(dataArr[i][:,2])) for i in 1:5];kmplt = plot(kmArr[1].events.time,kmArr[1].survival, seriestype = :steppost, xlabel = "Time", ylabel = "S(t)", label = "", ylim=(0,1), title="Kaplan-Meier estimator")
for i in 2:5
plot!(kmArr[i].events.time,kmArr[i].survival, seriestype = :steppost, label="")
end
kmplt
# Survival data elements
T = data[:,3]
δ = Int64.(data[:,2])
Z = [ [v] for v in data[:,1] ];# Survival data Struct for regression
datameta = RegressionSurvivalData( T, δ, Z)Right-censored survival data with covariates
Observations: 1506
Exact events: 453 (30.1%)
Right-censored: 1053 (69.9%)
Covariates: 1
Distinct times: 1359 (90.2%)
Distinct covariates: 5 (0.3%)# NTR prior specification
λ = 1/500
α = 10.010.0# Cox-NTR model baseline
baseline = ExponentialBaseline(λ)Neutral-to-the-right baseline
Specification: ExponentialBaseline(rate=0.002)# Plotting grid
t = collect(LinRange(0.0,356.0,200))
# Prior survival band
prior_cred = prior_credible_band( t, α, baseline )Credible band
Credible level: 95.0%
Estimate: Prior mean
Time grid: 0.0 to 356.0 (200 points)
Monte-Carlo paths: 3000
Median band width: 0.55965
Maximum band width: 0.57368# Prior baseline illustration
plt_prior = plot( prior_cred, fillalpha=0.3, label="Posterior NTR mean \n with 95% credible band",size=(600,400))
# Meta-analysis Cox regression fuction
f(c::Vector{Float64},z::Vector{Float64}) = exp(c[Int(z[1])]) f (generic function with 1 method)# Compilation runs
robbins_monro_mh_within_gibbs_tune(3, 3, x->NTRsurv.loglikelihood(x,α,baseline,f,datameta) + logpdf(Normal(0.0,exp(x[6])),x[1]) + logpdf(Normal(0.0,exp(x[6])),x[2]) + logpdf(Normal(0.0,exp(x[6])),x[3]) + logpdf(Normal(0.0,exp(x[6])),x[4]) + logpdf(Normal(0.0,exp(x[6])),x[5]) + logpdf(Normal(),x[6]),ones(6), ones(6), 0.4.*ones(6), 0.7)([[1.0, 1.0, 1.0, 1.0, 1.0, 1.0], [1.0634913484723583, 1.4467842341727166, 0.7817432769595243, 0.7817432769595243, 0.7817432769595243, 1.4467842341727166], [0.8835317934793192, 1.2019654612446036, 0.8188268965505427, 0.6494599514369258, 0.6494599514369258, 1.2019654612446036], [0.9176537961413386, 1.0329165932722524, 0.7036640533654728, 0.6745420982881861, 0.5581175018210232, 1.2483853738849104]], [1.227774143907035, 0.6933166388222056, 1.209116268891461, 1.0293289801565728, 1.0, 0.3046607807211055], -4751.9795043782415)# Robbins-Monro algorithm run for tuning of Metropolis-Hastings variance proposal
sd_prop_tuned, s₀_tune, lliks₀_tune = robbins_monro_mh_within_gibbs_tune(100, 50, x->NTRsurv.loglikelihood(x,α,baseline,f,datameta) + logpdf(Normal(0.0,exp(x[6])),x[1]) + logpdf(Normal(0.0,exp(x[6])),x[2]) + logpdf(Normal(0.0,exp(x[6])),x[3]) + logpdf(Normal(0.0,exp(x[6])),x[4]) + logpdf(Normal(0.0,exp(x[6])),x[5]) + logpdf(Normal(),x[6]),ones(6), ones(6), 0.4.*ones(6), 0.7)Progress: 100%[==================================================] Time: 0:01:02
([[1.0, 1.0, 1.0, 1.0, 1.0, 1.0], [0.8753210025791485, 0.8429457352093522, 0.8536022048457613, 0.9438496342853365, 0.8429457352093522, 1.005037726123811], [0.7918178072762002, 0.777093042648433, 0.757701166438233, 0.8538089976319071, 0.777093042648433, 1.0088473596933005], [0.7351629279038937, 0.7104186291945208, 0.6715915900632083, 0.8113244836195033, 0.7441583897262535, 1.0600342691364435], [0.6899291191011414, 0.6492977773975395, 0.6514637462396358, 0.7464435717184968, 0.7123685453416724, 1.0558353790120527], [0.6448705812457092, 0.6175865723869846, 0.6124730098725343, 0.7058660563687154, 0.6815336077185276, 1.0832335616254845], [0.613311176280374, 0.5904404819873905, 0.5886202565129157, 0.6819311260438822, 0.6653427470285829, 1.079841780357553], [0.5859139059192108, 0.5587200362436088, 0.5677052850171113, 0.6483746656754465, 0.6356211703996777, 1.061493830487597], [0.559311382724815, 0.5404125931184823, 0.5539387251624381, 0.6354313161612338, 0.6094274241344172, 1.0312227291861467], [0.5291794737691576, 0.5196949764016201, 0.5327025616483385, 0.6085877985902344, 0.6005589454576115, 1.0245260051538774] … [0.3192008649443277, 0.34925355208548836, 0.27958823592058657, 0.4046032012819604, 0.43419467012236046, 1.0044616726733036], [0.31931002510246215, 0.3496717622794503, 0.2782530746365226, 0.40508768936311335, 0.4332307497159202, 1.0030888253110197], [0.3188768612389196, 0.35008734834074234, 0.2781114642496415, 0.40316760551120984, 0.4319095358420953, 1.005986564381399], [0.3192530488484297, 0.34961593569344834, 0.27773697124889146, 0.40228578339216464, 0.4309648483797949, 1.003786235803886], [0.3188262938034092, 0.3485652945515271, 0.2769023356993194, 0.4054596171380096, 0.4311089897061228, 1.005802290449828], [0.31840318750600416, 0.3469489758398213, 0.27584717312248314, 0.4052577552034352, 0.4305368758875397, 1.0036341811600715], [0.3171985453541898, 0.34649184029453867, 0.27639325148139077, 0.40505743677973205, 0.4296154419217674, 1.0072797167940297], [0.3186031173353495, 0.34660525865281056, 0.27580589981249687, 0.4038660879558661, 0.4294045869626435, 1.0043171208830721], [0.3189657183322564, 0.34531229713146183, 0.2758955472819499, 0.40235952307126455, 0.4291953126137319, 1.0046435623713679], [0.31752781316217477, 0.3476605388386746, 0.2757620200689175, 0.4008689769954856, 0.4300272022566353, 1.0057790026656148]], [0.10643054723367396, 0.003515543190531073, 0.35271045119375, 0.031484393236567515, -0.1942295199548059, -1.1539916283790934], -4744.923512457405)# Metropolis-Hasting within Gibbs algorithm run
chain_s, _ = random_walk_mh_within_gibbs( 3000, x->NTRsurv.loglikelihood(x,α,baseline,f,datameta) + logpdf(Normal(0.0,exp(x[6])),x[1]) + logpdf(Normal(0.0,exp(x[6])),x[2]) + logpdf(Normal(0.0,exp(x[6])),x[3]) + logpdf(Normal(0.0,exp(x[6])),x[4]) + logpdf(Normal(0.0,exp(x[6])),x[5]) + logpdf(Normal(),x[6]), s₀_tune, lliks₀_tune, sd_prop_tuned[end])([[0.10643054723367396, 0.003515543190531073, 0.35271045119375, 0.031484393236567515, -0.1942295199548059, -1.1539916283790934], [0.3066969709056687, 0.003515543190531073, 0.43786270453451376, 0.3966936314153411, -0.1942295199548059, -1.1539916283790934], [0.3066969709056687, 0.027585044970382458, 0.43786270453451376, 0.10816824513692275, 0.20829011424683155, -1.1539916283790934], [0.3066969709056687, 0.027585044970382458, 0.43786270453451376, 0.10816824513692275, 0.20829011424683155, -1.1539916283790934], [0.3066969709056687, 0.027585044970382458, 0.43786270453451376, -0.1919781574286542, 0.20829011424683155, -1.227849915071102], [0.3066969709056687, 0.027585044970382458, 0.30216819342093115, -0.1919781574286542, -0.2091922553964424, -1.227849915071102], [0.07001719742212964, 0.027585044970382458, 0.27329195730624767, -0.1284565843914947, -0.2091922553964424, -1.00073343658083], [0.06979186736553307, -0.2312050052281685, 0.27329195730624767, -0.1284565843914947, 0.013310189818043217, -1.00073343658083], [0.06979186736553307, -0.2312050052281685, 0.27329195730624767, -0.1284565843914947, -0.14410214693805234, -1.5536507400274986], [0.06979186736553307, -0.2312050052281685, 0.27329195730624767, -0.1284565843914947, -0.14410214693805234, -1.8710460298115295] … [0.5622154408797904, -0.036211292070641665, 0.5818073867731198, 0.1482258938824208, 0.02998037582672844, -1.1979819643962857], [0.5622154408797904, -0.036211292070641665, 0.5818073867731198, 0.1482258938824208, 0.035700199555728424, -1.3198703388466195], [0.5253379088132468, 0.040607674179362285, 0.5818073867731198, -0.04296230207432103, 0.035700199555728424, -1.4918033323438304], [0.5253379088132468, 0.040607674179362285, 0.5818073867731198, -0.04296230207432103, 0.035700199555728424, -1.4918033323438304], [0.5253379088132468, 0.040607674179362285, 0.5059628127967446, -0.04296230207432103, 0.2935294927929757, -1.4918033323438304], [0.5253379088132468, -0.05810255853977141, 0.5059628127967446, -0.14174808848676762, 0.2935294927929757, -0.8421093452373483], [0.5253379088132468, -0.05810255853977141, 0.6017508349138226, 0.5188038103281862, 0.2935294927929757, -0.4732792269509778], [0.5253379088132468, 0.05783504318857073, 0.6017508349138226, 0.5188038103281862, 0.14174189304790955, -0.4732792269509778], [0.5253379088132468, -0.018269039890121533, 0.5198009146681734, 0.5188038103281862, 0.15426976777301812, -0.48802152307329033], [0.5253379088132468, 0.09261179485032216, 0.6498817004000894, 0.0580350753231999, 0.15426976777301812, -0.5682655219869355]], -4745.773738220244)# Acceptance rate of MCMC chains
acceptance_rate(chain_s)6-element Vector{Float64}:
0.3627875958652884
0.38779593197732576
0.3614538179393131
0.39613204401467156
0.37279093031010335
0.3844614871623875# Plot of MCMC chains for assessing converegence of posterior frailties
plot([ v[1] for v in chain_s],label="Tune M-H chain 1",legend=:outerright)
plot!([ v[2] for v in chain_s],label="Tune M-H chain 2")
plot!([ v[3] for v in chain_s],label="Tune M-H chain 3")
plot!([ v[4] for v in chain_s],label="Tune M-H chain 4")
plot!([ v[5] for v in chain_s],label="Tune M-H chain 5")
# Plot of MCMC chain for frailty variance parameters
plot([ exp(v[6]) for v in chain_s],label="Tune M-H chain 5")
# Cohort covariate vector
zv = [ [1.0], [2.0], [3.0], [4.0], [5.0] ];# Fully Bayesian model fit
cox_model = CoxNeutralToTheRightModel( chain_s, α, baseline, f, datameta)Cox Neutral-to-the-right model
Coefficient draws: 3000
Risk score: Custom g
Gamma process α: 10.0
Center at Z=zeros: ExponentialBaseline(rate=0.002)
Data: SurvivalData(obs.=1506, exact=453, censored=1053, distinct T = 1359, distinct Z = 5)# Computation of 95% posterior credible bands
NTR_bands = posterior_credible_band( t, zv, cox_model)5-element Vector{NTRsurv.CredibleBand}:
CredibleBand( level=0.050000000000000044, estimate=posterior median, covariate Z=[1.0], model=fully Bayesian, draws=3000)
CredibleBand( level=0.050000000000000044, estimate=posterior median, covariate Z=[2.0], model=fully Bayesian, draws=3000)
CredibleBand( level=0.050000000000000044, estimate=posterior median, covariate Z=[3.0], model=fully Bayesian, draws=3000)
CredibleBand( level=0.050000000000000044, estimate=posterior median, covariate Z=[4.0], model=fully Bayesian, draws=3000)
CredibleBand( level=0.050000000000000044, estimate=posterior median, covariate Z=[5.0], model=fully Bayesian, draws=3000)# Illustration of posterior cohort survival estimations
p_meta_1 = plot( NTR_bands[1], fillalpha=0.3, label="",size=(600,400))
plot!(p_meta_1,kmArr[1].events.time,kmArr[1].survival, seriestype = :steppost, xlabel = "t", ylabel = "S(t)", label = "", c=:black, ylim=(0,1), title="Cohort 1", leftmargin=5mm)
p_meta_2 = plot( NTR_bands[2], fillalpha=0.3, label="",size=(600,400))
plot!(p_meta_2,kmArr[2].events.time,kmArr[2].survival, seriestype = :steppost, xlabel = "t", ylabel = "", label = "", c=:black, ylim=(0,1), title="Cohort 2")
p_meta_3 = plot( NTR_bands[3], fillalpha=0.3, label="",size=(600,400))
plot!(p_meta_3,kmArr[3].events.time,kmArr[3].survival, seriestype = :steppost, xlabel = "t", ylabel = "", label = "", c=:black, ylim=(0,1), title="Cohort 3")
p_meta_4 = plot( NTR_bands[4], fillalpha=0.3, label="",size=(600,400))
plot!(p_meta_4,kmArr[4].events.time,kmArr[4].survival, seriestype = :steppost, xlabel = "t", ylabel = "S(t)", label = "", c=:black, ylim=(0,1), title="Cohort 4", leftmargin=5mm)
p_meta_5 = plot( NTR_bands[5], fillalpha=0.3, label="",size=(600,400))
plot!(p_meta_5,kmArr[5].events.time,kmArr[5].survival, seriestype = :steppost, xlabel = "t", ylabel = "", label = "", c=:black, ylim=(0,1), title="Cohort 5")
p_legend = plot(
[0], [0],
ribbon = (NTR_bands[5].m .- NTR_bands[5].d, NTR_bands[5].u .- NTR_bands[5].m),
fillalpha = 0.3,
label = "Posterior NTR mean \nwith 95% credible band",
framestyle = :none,
axis = false,
legend=(0.15,0.5),
foreground_color_legend=:transparent
)
plot!([0],[0],c=:black, label="Kaplan-Meier estimator ")
p_meta_cohorts = plot(p_meta_1, p_meta_2, p_meta_3, p_meta_4, p_meta_5, p_legend, layout = (2, 3), size = (900, 600),sup_title="")