updated surv from cox ph with alternative default

based on aalen-johansen risk

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

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"""

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)

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