From 9d41c1a82866b2d460ba3f60c4018c6923d6d426 Mon Sep 17 00:00:00 2001 From: Alex Keil Date: Fri, 29 Sep 2023 09:50:23 -0400 Subject: [PATCH] updated surv from cox ph with alternative default based on aalen-johansen risk --- src/coxmodel.jl | 20 ++++++++++++++------ src/docstr.jl | 6 ++++++ src/examples.jl | 3 ++- test/runtests.jl | 2 +- 4 files changed, 23 insertions(+), 8 deletions(-) diff --git a/src/coxmodel.jl b/src/coxmodel.jl index 0b3569a..93ef974 100644 --- a/src/coxmodel.jl +++ b/src/coxmodel.jl @@ -778,10 +778,10 @@ end # fitting functions for PHSurv objects ##################################################################################################################### -function _fit!(m::M; coef_vectors = nothing, pred_profile = nothing) where {M<:PHSurv} +function _fit!(m::M; coef_vectors = nothing, pred_profile = nothing, method="aalen-johansen") where {M<:PHSurv} hr = ones(Float64, length(m.eventtypes)) ch::Float64 = 0.0 - lsurv::Float64 = 1.0 + surv_previous::Float64 = 1.0 if (isnothing(coef_vectors)) coef_vectors = [coef(fit) for fit in m.fitlist] end @@ -795,14 +795,22 @@ function _fit!(m::M; coef_vectors = nothing, pred_profile = nothing) where {M<:P @inbounds for (j, d) in enumerate(m.eventtypes) if m.event[i] == d m.basehaz[i] *= hr[j] # baseline hazard times hazard ratio - m.risk[i, j] = lci[j] + m.basehaz[i] * lsurv + m.risk[i, j] = lci[j] + m.basehaz[i] * surv_previous else m.risk[i, j] = lci[j] end end - ch += m.basehaz[i] - m.surv[i] = exp(-ch) - lsurv = m.surv[i] + if lowercase(method[1:3]) == "aal" + ## aalen-johansen + m.surv[i] = surv_previous - surv_previous * m.basehaz[i] + elseif lowercase(method[1:3]) == "che" + #Cheng, Fine and Wei + ch += m.basehaz[i] + m.surv[i] = exp(-ch) + else + throw("method $method not recognized (use 'aalen-johansen' or 'cheng-fine-wei')") + end + surv_previous = m.surv[i] lci = m.risk[i, :] end m.fit = true diff --git a/src/docstr.jl b/src/docstr.jl index 4c40c06..205591c 100644 --- a/src/docstr.jl +++ b/src/docstr.jl @@ -481,7 +481,13 @@ using Random # Risk at x=1, z=0.5 res3 = risk_from_coxphmodels([ft1, ft2], pred_profile=[1.0, 0.5]) + + # method using Aalen-Johansen analog (default) + res3 = risk_from_coxphmodels([ft1, ft2], pred_profile=[1.0, 0.5]) + # method using Cheng, Fine and Wei's cumulative hazards based approach + res3b = risk_from_coxphmodels([ft1, ft2], pred_profile=[1.0, 0.5], method="cheng-fine-wei") ``` +Cheng SC, Fine JP, Wei LJ. Prediction of Cumulative Incidence Function under the Proportional Hazards Model. Biometrics. 1998;54:219–228. """ DOC_E_YEARSOFLIFELOST = raw""" diff --git a/src/examples.jl b/src/examples.jl index c3245ef..f8c5d8f 100644 --- a/src/examples.jl +++ b/src/examples.jl @@ -122,7 +122,7 @@ hcat(d[1:15], event[1:15]) # cause-specific Cox models ft1 = fit(PHModel, X, int, outt, d .* (event .== 1)) -ft2 = fit(PHModel, X, int, outt, d .* (event .== 2)) +ft2 = fit(PHModel, X[findall(event .!== 1),:], int[findall(event .!== 1)], outt[findall(event .!== 1)], d .* (event .== 2)[findall(event .!== 1)]) ft2_equivalent = fit(PHModel, X, int, outt, (event .== 2)) @@ -140,6 +140,7 @@ risk2 = aalen_johansen(int, outt, event) # risk at the baseline values of each covariate (note high risk inferred by both covariates being protective of both event types) fit(PHSurv, [ft1, ft2]) risk1_ref = risk_from_coxphmodels([ft1, ft2]) +sum(risk1_ref.risk, dims=2) # risk at mean predictor values Xpred = mean(X, dims = 1) # (a "profile" for the average individual) diff --git a/test/runtests.jl b/test/runtests.jl index 02c0eda..ee3045d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -534,7 +534,7 @@ using Random, Tables rfromc = risk_from_coxphmodels( [ft1, ft2], coef_vectors = [coef(res), coef(res2)], - pred_profile = [0 0 -1], + pred_profile = [0 0 -1], method="che" ) println(rfromc) refrisk = rfromc.risk