diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index f8597d5..9a23481 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.11.1","generation_timestamp":"2024-11-15T13:53:14","documenter_version":"1.8.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.11.1","generation_timestamp":"2024-11-15T14:59:03","documenter_version":"1.8.0"}} \ No newline at end of file diff --git a/dev/Likelihood/index.html b/dev/Likelihood/index.html index 0651fb3..1047dfb 100644 --- a/dev/Likelihood/index.html +++ b/dev/Likelihood/index.html @@ -71,4 +71,4 @@ \\ S(t|\lambda,k)=& 1-I(k,\lambda t)\\ \ln S(t|\alpha,\kappa)=&\ln(1-I(\exp(\kappa)), \exp(\ln(t)-\alpha))\\ -\end{aligned}$

Where we define $I(k,s)$ as the upper incomplete gamma function ratio given by

\[I(k,s) = \frac{\int_o^s t^{k-1}\exp(-t)dt}{\Gamma(k)}\]

(see gamma_inc function from SpecialFunctions.jl)

Log-logistic distribution

(not yet implemented)

Semi-parametric partial likelihoods

In the Cox model, the partial-likelihoods are used in place of the likelihood function. These models are are modeled directly in terms of hazard ratios, allowing that the baseline hazard can be an arbitrary distribution. The Cox models implemented here are semi-parametric because they include a combination of parametric (hazard ratios) and non-parametric (baseline hazard) components. Cox's original likelihood is used here, and, in place of tied survival times, two different options are implemented for addressing ties. See the survival[surv] package vignette for original citations and methods for baseline hazard and partial-likelihood calculations.

Efron's partial likelihood

This is the default in coxph (documentation in progress)

Breslow's partial likelihood

Documentation in progress

Time-varying covariates

Documentation in progress

Numerical algorithms

Fitting algorithms include direct calculation, hard-coded Newton-Raphson algorithms, and optimization algorithms from the Optim.jl module:

Non-parametric analysis

Semi-parametric analysis

Parametric analysis

+\end{aligned}$

Where we define $I(k,s)$ as the upper incomplete gamma function ratio given by

\[I(k,s) = \frac{\int_o^s t^{k-1}\exp(-t)dt}{\Gamma(k)}\]

(see gamma_inc function from SpecialFunctions.jl)

Log-logistic distribution

(not yet implemented)

Semi-parametric partial likelihoods

In the Cox model, the partial-likelihoods are used in place of the likelihood function. These models are are modeled directly in terms of hazard ratios, allowing that the baseline hazard can be an arbitrary distribution. The Cox models implemented here are semi-parametric because they include a combination of parametric (hazard ratios) and non-parametric (baseline hazard) components. Cox's original likelihood is used here, and, in place of tied survival times, two different options are implemented for addressing ties. See the survival[surv] package vignette for original citations and methods for baseline hazard and partial-likelihood calculations.

Efron's partial likelihood

This is the default in coxph (documentation in progress)

Breslow's partial likelihood

Documentation in progress

Time-varying covariates

Documentation in progress

Numerical algorithms

Fitting algorithms include direct calculation, hard-coded Newton-Raphson algorithms, and optimization algorithms from the Optim.jl module:

Non-parametric analysis

Semi-parametric analysis

Parametric analysis

diff --git a/dev/coxmodel/index.html b/dev/coxmodel/index.html index 6c6869e..833f99b 100644 --- a/dev/coxmodel/index.html +++ b/dev/coxmodel/index.html @@ -294,4 +294,4 @@ ─────┼──────────────────────────────────────────── 1 │ 0.872755 0.954257 0.895928 0.912742 2 │ -0.443074 -0.439597 -0.475232 -0.485941 - 3 │ 1.32749 1.32809 1.28391 1.32294
+ 3 │ 1.32749 1.32809 1.28391 1.32294
diff --git a/dev/index.html b/dev/index.html index 84ce6a9..6ea5717 100644 --- a/dev/index.html +++ b/dev/index.html @@ -41,7 +41,7 @@ # can use the ID type to refer to units with multiple observations id, int, outt, data = dgm(MersenneTwister(), 1000, 10; regimefun = int_0) -LSurvivalResp(int, outt, data[:,4], ID.(id))

Index of functions

Function help

LSurvival.AbstractLSurvivalParms — Type

AbstractLsurvParms

Abstract type representing a model predictors and coefficient parameters

source
LSurvival.AbstractLSurvivalResp — Type

AbstractLsurvResp

Abstract type representing a model response vector

source
LSurvival.AbstractNPSurv — Type

Abstract type for non-parametric survival models, including Kaplan-Meier, Aalen Johansen, and Cox-model based estimates of survival using an Aalen-Johansen-like estimator

source
LSurvival.AbstractPH — Type

Abstract type for proportional hazards models

source
LSurvival.AbstractPSModel — Type

AbstractPS

Abstract type for parametric survival models

source
LSurvival.AbstractSurvDist — Type

AbstractSurvDist

Abstract type for parametric survival distributions

source
LSurvival.ID — Type

Type for identifying individuals in survival outcomes.

Used for the id argument in

  • Outcome types: LSurvivalResp, LSurvivalCompResp
  • Model types: PHModel, KMRisk, AJRisk

Accepts any Number or String. There is no significance to having this particular struct, but it enables easier use of multiple dispatch.

[ID(i) for i in 1:10]
source
LSurvival.LSurvivalCompResp — Type

Outcome type for competing risk survival outcomes subject to left truncation and right censoring (not generally needed for users)

Parameters

  • enter Time at observation start

  • exit Time at observation end

  • y event occurrence in observation

  • wts observation weights

  • eventtimes unique event times

  • origin origin on the time scale

  • id person level identifier (must be wrapped in ID() function)

  • eventtypes vector of unique event types

  • eventmatrix matrix of indicators on the observation level

    Signatures:

 struct LSurvivalCompResp{
+LSurvivalResp(int, outt, data[:,4], ID.(id))

Index of functions

Function help

LSurvival.AbstractNPSurv — Type

Abstract type for non-parametric survival models, including Kaplan-Meier, Aalen Johansen, and Cox-model based estimates of survival using an Aalen-Johansen-like estimator

source
LSurvival.ID — Type

Type for identifying individuals in survival outcomes.

Used for the id argument in

  • Outcome types: LSurvivalResp, LSurvivalCompResp
  • Model types: PHModel, KMRisk, AJRisk

Accepts any Number or String. There is no significance to having this particular struct, but it enables easier use of multiple dispatch.

[ID(i) for i in 1:10]
source
LSurvival.LSurvivalCompResp — Type

Outcome type for competing risk survival outcomes subject to left truncation and right censoring (not generally needed for users)

Parameters

  • enter Time at observation start

  • exit Time at observation end

  • y event occurrence in observation

  • wts observation weights

  • eventtimes unique event times

  • origin origin on the time scale

  • id person level identifier (must be wrapped in ID() function)

  • eventtypes vector of unique event types

  • eventmatrix matrix of indicators on the observation level

    Signatures:

 struct LSurvivalCompResp{
  E<:AbstractVector,
  X<:AbstractVector,
  Y<:AbstractVector,
@@ -83,7 +83,7 @@
  )
 LSurvivalCompResp(
   exit::X,
   y::Y,
-  ) where {X<:Vector,Y<:Union{Vector{<:Real},BitVector}}
source
LSurvival.LSurvivalResp — Type

Outcome type for survival outcome subject to left truncation and right censoring.

Will not generally be needed by users

Parameters

  • enter: Time at observation start
  • exit: Time at observation end
  • y: event occurrence in observation
  • wts: observation weights
  • eventtimes: unique event times
  • origin: origin on the time scale
  • id: person level identifier (must be wrapped in ID() function)
 struct LSurvivalResp{
+  ) where {X<:Vector,Y<:Union{Vector{<:Real},BitVector}}
source
LSurvival.LSurvivalResp — Type

Outcome type for survival outcome subject to left truncation and right censoring.

Will not generally be needed by users

Parameters

  • enter: Time at observation start
  • exit: Time at observation end
  • y: event occurrence in observation
  • wts: observation weights
  • eventtimes: unique event times
  • origin: origin on the time scale
  • id: person level identifier (must be wrapped in ID() function)
 struct LSurvivalResp{
  E<:AbstractVector,
  X<:AbstractVector,
  Y<:AbstractVector,
@@ -127,7 +127,7 @@
   y::Y,
   ) where {E<:Vector,X<:Vector,Y<:Union{Vector{<:Real},BitVector}}
 LSurvivalResp(exit::X, y::Y) where {X<:Vector,Y<:Vector}

Examples

  # no late entry
   LSurvivalResp([.5, .6], [1,0])
-
source
LSurvival.PHModel — Type

PHModel: Mutable object type for proportional hazards regression (not generally needed for users)

Parameters

  • R Survival response

  • P # parameters

  • ties String: "efron" or "breslow"

  • fit Bool: logical for whether the model has been fitted

  • bh AbstractMatrix: baseline hazard estimates

    Signatures

 mutable struct PHModel{G<:LSurvivalResp,L<:AbstractLSurvivalParms} <: AbstractPH
+
source
LSurvival.PHModel — Type

PHModel: Mutable object type for proportional hazards regression (not generally needed for users)

Parameters

  • R Survival response

  • P # parameters

  • ties String: "efron" or "breslow"

  • fit Bool: logical for whether the model has been fitted

  • bh AbstractMatrix: baseline hazard estimates

    Signatures

 mutable struct PHModel{G<:LSurvivalResp,L<:AbstractLSurvivalParms} <: AbstractPH
  R::G        # Survival response
  P::L        # parameters
  ties::String #"efron" or"breslow"
@@ -152,7 +152,7 @@
  R = LSurvivalResp(enter, t, Int.(d), wt)
  P = PHParms(X)
  mf = PHModel(R,P)
-  LSurvival._fit!(mf)
source
LSurvival.PHSurv — Type

Mutable type for proportional hazards models (not generally needed by users)

PHSsurv: Object type for proportional hazards regression

surv::Vector{Float64} risk::Matrix{Float64} basehaz::Vector{Float64} event::Vector{Float64}

  • fitlist: vector of PHSurv objects (Cox model fits)
  • eventtypes: vector of unique event types
  • times: unique event times
  • surv: Overall survival at each time
  • risk: Cause-specific risk at each time (1 for each outcome type)
  • basehaz: baseline hazard for a specific event type
  • event: value of event type that occurred at each time

Methods: fit, show

mutable struct PHSurv{G<:Array{T} where {T<:PHModel}} <: AbstractNPSurv
+  LSurvival._fit!(mf)
source
LSurvival.PHSurv — Type

Mutable type for proportional hazards models (not generally needed by users)

PHSsurv: Object type for proportional hazards regression

surv::Vector{Float64} risk::Matrix{Float64} basehaz::Vector{Float64} event::Vector{Float64}

  • fitlist: vector of PHSurv objects (Cox model fits)
  • eventtypes: vector of unique event types
  • times: unique event times
  • surv: Overall survival at each time
  • risk: Cause-specific risk at each time (1 for each outcome type)
  • basehaz: baseline hazard for a specific event type
  • event: value of event type that occurred at each time

Methods: fit, show

mutable struct PHSurv{G<:Array{T} where {T<:PHModel}} <: AbstractNPSurv
 fitlist::G        
 eventtypes::AbstractVector
 times::AbstractVector
@@ -163,13 +163,13 @@
 end
 
 PHSurv(fitlist::Array{T}, eventtypes) where {T<:PHModel}
-PHSurv(fitlist::Array{T}) where {T<:PHModel}
source
LSurvival.Strata — Type

Type for identifying individuals in survival outcomes. Used for the strata argument in PHModel (not yet implemented)

Accepts any Number or String. There is no significance to having this particular struct, but it enables easier use of multiple dispatch.

[Strata(i) for i in 1:10]
source
LSurvival._update_PHParms! — Method

Update the partial likelihood, gradient and Hessian values from a Cox model fit (used during fitting, not generally useful for users).

Uses Breslow's or Efron's partial likelihood.

Updates over all observations

Signature

_update_PHParms!(
+PHSurv(fitlist::Array{T}) where {T<:PHModel}
source
LSurvival.Strata — Type

Type for identifying individuals in survival outcomes. Used for the strata argument in PHModel (not yet implemented)

Accepts any Number or String. There is no significance to having this particular struct, but it enables easier use of multiple dispatch.

[Strata(i) for i in 1:10]
source
LSurvival._update_PHParms! — Method

Update the partial likelihood, gradient and Hessian values from a Cox model fit (used during fitting, not generally useful for users).

Uses Breslow's or Efron's partial likelihood.

Updates over all observations

Signature

_update_PHParms!(
  m::M,
  # big indexes
  ne::I,
  caseidxs::Vector{Vector{T}},
  risksetidxs::Vector{Vector{T}},
- ) where {M<:AbstractPH,I<:Int,T<:Int}

updatePHParms!(m, risksetidxs, caseidxs, ne, den)

source
LSurvival.aalen_johansen — Method

Aalen-Johansen estimator for cumulative cause-specific risk (in the presence of competing events)

Signatures

 StatsBase.fit!(m::T; kwargs...) where {T<:AbstractNPSurv}
+ ) where {M<:AbstractPH,I<:Int,T<:Int}

updatePHParms!(m, risksetidxs, caseidxs, ne, den)

source
LSurvival.aalen_johansen — Method

Aalen-Johansen estimator for cumulative cause-specific risk (in the presence of competing events)

Signatures

 StatsBase.fit!(m::T; kwargs...) where {T<:AbstractNPSurv}
 
  aalen_johansen(enter::AbstractVector, exit::AbstractVector, y::AbstractVector,
    ; kwargs...)
@@ -179,7 +179,7 @@
 enter = zeros(length(t));
    # event variable is coded 0[referent],1,2
 m = fit(AJSurv, enter, t, event)
-mw = fit(AJSurv, enter, t, event, wts=wt)

or, equivalently:

aalen_johansen(enter, t, event, wts=wt)
source
LSurvival.bootstrap — Method

Bootstrapping coefficients of a proportional hazards model

Signatures

# single bootstrap draw, keeping the entire object
+mw = fit(AJSurv, enter, t, event, wts=wt)

or, equivalently:

aalen_johansen(enter, t, event, wts=wt)
source
LSurvival.bootstrap — Method

Bootstrapping coefficients of a proportional hazards model

Signatures

# single bootstrap draw, keeping the entire object
 bootstrap(rng::MersenneTwister, m::PHModel)
 bootstrap(m::PHModel)
 # muliple bootstrap draws, keeping only coefficient estimates
@@ -233,7 +233,7 @@
 [stddev_finite(mb[:,i]) for i in 1:2]
 ## asymptotic standard error
 stderror(mainfit)
-
source
LSurvival.bootstrap — Method

Bootstrap methods for Aalen-Johansen cumulative risk estimator

Signatures

 # single bootstrap draw, keeping the entire object
+
source
LSurvival.bootstrap — Method

Bootstrap methods for Aalen-Johansen cumulative risk estimator

Signatures

 # single bootstrap draw, keeping the entire object
  bootstrap(rng::MersenneTwister, m::AJSurv)
  bootstrap(m::AJSurv)
 
@@ -254,7 +254,7 @@
 
 aj1.R
 aj2.R
-
source
LSurvival.bootstrap — Method

Bootstrap methods for Kaplan-Meier survival curve estimator

Signatures

 # single bootstrap draw, keeping the entire object
+
source
LSurvival.bootstrap — Method

Bootstrap methods for Kaplan-Meier survival curve estimator

Signatures

 # single bootstrap draw, keeping the entire object
  bootstrap(rng::MersenneTwister, m::KMSurv)
  bootstrap(m::KMSurv)
 
@@ -276,7 +276,7 @@
 
 km1.R
 km2.R
-
source
LSurvival.bootstrap — Method

Bootstrapping coefficients of a proportional hazards model

Signatures

# single bootstrap draw, keeping the entire object
+
source
LSurvival.bootstrap — Method

Bootstrapping coefficients of a proportional hazards model

Signatures

# single bootstrap draw, keeping the entire object
 bootstrap(rng::MersenneTwister, m::PHModel)
 bootstrap(m::PHModel)
 # muliple bootstrap draws, keeping only coefficient estimates
@@ -330,7 +330,7 @@
 [stddev_finite(mb[:,i]) for i in 1:2]
 ## asymptotic standard error
 stderror(mainfit)
-
source
LSurvival.bootstrap — Method

Bootstrap sampling of a proportional hazards predictor object

using LSurvival, Random
 
 id, int, outt, data =
 LSurvival.dgm(MersenneTwister(1212), 20, 5; afun = LSurvival.int_0)
@@ -346,27 +346,27 @@
 
 Mod = PHModel(R2, P2)
 LSurvival._fit!(Mod, start=Mod.P._B)
-
source
LSurvival.bootstrap — Method

Bootstrapping sampling of a competing risk survival response

Signatures

bootstrap(rng::MersenneTwister, R::T) where {T<:LSurvivalCompResp}
+
source
LSurvival.bootstrap — Method

Bootstrapping sampling of a competing risk survival response

Signatures

bootstrap(rng::MersenneTwister, R::T) where {T<:LSurvivalCompResp}
 bootstrap(R::T) where {T<:LSurvivalCompResp}
z,x,t,d,event,weights =
 LSurvival.dgm_comprisk(MersenneTwister(1212), 300)
 enter = zeros(length(event))
 
 # survival outcome:
 R = LSurvivalCompResp(enter, t, event, weights, ID.(collect(1:length(t))))    # specification with ID only
-bootstrap(R) # note that entire observations/clusters identified by id are kept
source
LSurvival.bootstrap — Method

Bootstrapping sampling of a survival response

id, int, outt, data =
+bootstrap(R) # note that entire observations/clusters identified by id are kept
source
LSurvival.bootstrap — Method

Bootstrapping sampling of a survival response

id, int, outt, data =
 LSurvival.dgm(MersenneTwister(1212), 20, 5; afun = LSurvival.int_0)
 
 d, X = data[:, 4], data[:, 1:3]
 weights = rand(length(d))
 
 # survival outcome:
-R = LSurvivalResp(int, outt, d, ID.(id))    # specification with ID only
source
LSurvival.calcp — Method

Two-tailed p-value for a (null) standard normal distribution depends on SpecialFunctions https://en.wikipedia.org/wiki/Normal_distribution

    calcp(z) = 1.0 - SpecialFunctions.erf(abs(z) / sqrt(2))
+R = LSurvivalResp(int, outt, d, ID.(id))    # specification with ID only
source
LSurvival.calcp — Method

Two-tailed p-value for a (null) standard normal distribution depends on SpecialFunctions https://en.wikipedia.org/wiki/Normal_distribution

    calcp(z) = 1.0 - SpecialFunctions.erf(abs(z) / sqrt(2))
 
-    calcp(1.96)

source
LSurvival.cdfchisq — Method

quantile function for a chi-squared distribution depends on SpecialFunctions https://en.wikipedia.org/wiki/Chi-squared_distribution

Source code, example:

    cdfchisq(df, x) = SpecialFunctions.gamma_inc(df / 2, x / 2, 0)[1]
+    calcp(1.96)

source
LSurvival.cdfchisq — Method

quantile function for a chi-squared distribution depends on SpecialFunctions https://en.wikipedia.org/wiki/Chi-squared_distribution

Source code, example:

    cdfchisq(df, x) = SpecialFunctions.gamma_inc(df / 2, x / 2, 0)[1]
 
-    cdfchisq(3, 3.45)
source
LSurvival.cdfnorm — Method

quantile function for a standard normal distribution depends on SpecialFunctions https://en.wikipedia.org/wiki/Normal_distribution

Source code, example:

    cdfnorm(z) = 0.5 * (1 + SpecialFunctions.erf(z / sqrt(2)))
+    cdfchisq(3, 3.45)
source
LSurvival.cdfnorm — Method

quantile function for a standard normal distribution depends on SpecialFunctions https://en.wikipedia.org/wiki/Normal_distribution

Source code, example:

    cdfnorm(z) = 0.5 * (1 + SpecialFunctions.erf(z / sqrt(2)))
     
-    cdfnorm(1.96)
source
LSurvival.coxph — Method

Fit method for AbstractPH objects (Cox models)

Keyword arguments (used here, and passed on to internal structs)

  • ties "breslow" or "efron" (default)

  • wts observation weights

  • ties "breslow" or "efron" (default)

  • offset not currently used at all

  • fitargs arguments passed to other structs, which include

    • id cluster or individual level ID (defaults to a unique value for each row of data) see note below on ID
    • contrasts StatsModel style contrasts (dicts) that can be used for variable transformations/indicator variable creation (e.g. https://juliastats.org/StatsModels.jl/stable/contrasts/)

  • Arguments passed onto fitting routine:

    • eps (default: Float64 = 1e-9) tolerance for declaring convergence. Model is determined to be converged when relative change in log-partial likelihood is < eps .
    • getbasehaz (default: true): estimate baseline hazard
    • start (default: nothing) nothing, or vector of floats corresponding to initial values for parameters. Note that this defaults to a vector of zeros when set to nothing, and setting to other values invalidates some of the test statistics reported by default with coxph.
    • keepx (default: true) logical. Keep design matrix in AbstractPH object output (set to false for slight computational gain).
    • keepy (default: true)logical. Keep outcome in AbstractPH object output (set to false for slight computational gain).
    • bootstrap_sample (default: false) Fit the model to a bootstrap sample of the data (not generally used by end-users, but provides some convenience toward bootstrap variance estimation).
    • bootstrap_rng (default: Random.MersenneTwister()) Random number seed used when drawing a bootstrap sample of the data (not generally used by end-users, but provides some convenience toward bootstrap variance estimation).

    Signatures

  fit(::Type{M},
+    cdfnorm(1.96)
source
LSurvival.coxph — Method

Fit method for AbstractPH objects (Cox models)

Keyword arguments (used here, and passed on to internal structs)

  • ties "breslow" or "efron" (default)

  • wts observation weights

  • ties "breslow" or "efron" (default)

  • offset not currently used at all

  • fitargs arguments passed to other structs, which include

    • id cluster or individual level ID (defaults to a unique value for each row of data) see note below on ID
    • contrasts StatsModel style contrasts (dicts) that can be used for variable transformations/indicator variable creation (e.g. https://juliastats.org/StatsModels.jl/stable/contrasts/)

  • Arguments passed onto fitting routine:

    • eps (default: Float64 = 1e-9) tolerance for declaring convergence. Model is determined to be converged when relative change in log-partial likelihood is < eps .
    • getbasehaz (default: true): estimate baseline hazard
    • start (default: nothing) nothing, or vector of floats corresponding to initial values for parameters. Note that this defaults to a vector of zeros when set to nothing, and setting to other values invalidates some of the test statistics reported by default with coxph.
    • keepx (default: true) logical. Keep design matrix in AbstractPH object output (set to false for slight computational gain).
    • keepy (default: true)logical. Keep outcome in AbstractPH object output (set to false for slight computational gain).
    • bootstrap_sample (default: false) Fit the model to a bootstrap sample of the data (not generally used by end-users, but provides some convenience toward bootstrap variance estimation).
    • bootstrap_rng (default: Random.MersenneTwister()) Random number seed used when drawing a bootstrap sample of the data (not generally used by end-users, but provides some convenience toward bootstrap variance estimation).

    Signatures

  fit(::Type{M},
   X::AbstractMatrix,#{<:FP},
   enter::AbstractVector{<:Real},
   exit::AbstractVector{<:Real},
@@ -435,7 +435,7 @@
  stderror(ft2w, type="robust")      # robust variance (INCORRECT)
  
  ft2w # the coefficient table now shows incorrect confidence intervals and test statistics
-
source
LSurvival.ddloglik! — Method

Hessian contribution for an observation in a parametric survival model

    d = m.d
     i = 1
     enter = m.R.enter[i]
     exit = m.R.exit[i]
@@ -443,7 +443,7 @@
     wts = m.R.wts[i]
     x = m.P.X[i,:]
     θ=[1,0,.4]
-
source
LSurvival.dgm — Method

Generating discrete survival data without competing risks

Signatures

dgm(rng::MersenneTwister, n::Int, maxT:Int; afun = int_0, yfun = yprob, lfun = lprob)
+
source
LSurvival.dgm — Method

Generating discrete survival data without competing risks

Signatures

dgm(rng::MersenneTwister, n::Int, maxT:Int; afun = int_0, yfun = yprob, lfun = lprob)
 
 dgm(n::Int, maxT::Int; kwargs...)

Usage: dgm(rng, n, maxT;afun=int0, yfun=yprob, lfun=lprob) dgm(n, maxT;afun=int0, yfun=yprob, lfun=lprob)

Where afun, yfun, and lfun are all functions that take arguments v,l,a and output time-specific values of a, y, and l respectively Example:


 expit(mu) =  inv(1.0+exp(-mu))
@@ -460,19 +460,19 @@
 expit(-3 + 2*v + 0*l + 2*a)
 end
   # 10 individuals followed for up to 5 times
-LSurvival.dgm(10, 5;afun=aprob, yfun=yprob, lfun=lprob)
source
LSurvival.dgm_comprisk — Method

Generating continuous survival data with competing risks

Signatures

dgm_comprisk(rng::MersenneTwister, n::Int)
+LSurvival.dgm(10, 5;afun=aprob, yfun=yprob, lfun=lprob)
source
LSurvival.dgm_comprisk — Method

Generating continuous survival data with competing risks

Signatures

dgm_comprisk(rng::MersenneTwister, n::Int)
 
 dgm_comprisk(n::Int)
    - rng = random number generator    
     - n = sample size

Example:

using LSurvival
 # 100 individuals with two competing events
 z,x,t,d,event,weights = LSurvival.dgm_comprisk(100)
     
-
source
LSurvival.dgm_phmodel — Method

Proportional hazards model from a Weibull distribution with scale parameter λ

dgm_phmodel(rng::MersenneTwister, n::Int; 
+
source
LSurvival.dgm_phmodel — Method

Proportional hazards model from a Weibull distribution with scale parameter λ

dgm_phmodel(rng::MersenneTwister, n::Int; 
     λ=1.25,
     β=[0.0, 0.0]
     )

keyword parameters:

  • λ: Weibull scale parameter
  • β: vector of regression coefficients
rng = MersenneTwister()
 X, t, d, _ = dgm_phmodel(2000; λ=1.25,β=[1.0, -0.5])
-coxph(@formula(Surv(t0,t,d)~x+z), (t=t,t0=t.*0,d=d,x=X[:,1],z=X[:,2]))
source
LSurvival.dloglik! — Method

Gradient contribution for an observation in a parametric survival model

    d = m.d
+coxph(@formula(Surv(t0,t,d)~x+z), (t=t,t0=t.*0,d=d,x=X[:,1],z=X[:,2]))
source
LSurvival.dloglik! — Method

Gradient contribution for an observation in a parametric survival model

    d = m.d
     i = 1
     enter = m.R.enter[i]
     exit = m.R.exit[i]
@@ -480,7 +480,7 @@
     wts = m.R.wts[i]
     x = m.P.X[i,:]
     θ=[1,0,.4]
-
source
LSurvival.dlpdf_gengamma — Method

α=0.1 ρ =-1.2 κ=1.9 t = 2.0

exp((log(t) - α)exp(-ρ) - ρ) - exp(κ - ρ) (α - log(t))exp(κ - ρ) + (log(t) - α)exp((log(t) - α)exp(-ρ) - ρ) - 1 (log(t) - α)exp(κ - ρ) - exp(κ)SpecialFunctions.digamma(exp(κ))

source
LSurvival.jackknife — Method

Obtain jackknife (leave-one-out) estimates from a Aalen-Johansen risk curve (risk at end of follow-up) by refitting the model n times

Signatures

jackknife(m::M;kwargs...) where {M<:AJSurv}
source
LSurvival.jackknife — Method

Obtain jackknife (leave-one-out) estimates from a Kaplan-Meier survival curve (survival at end of follow-up) by refitting the model n times

Signatures

jackknife(m::M;kwargs...) where {M<:KMSurv}
using LSurvival, Random, StatsBase
+
source
LSurvival.dlpdf_gengamma — Method

α=0.1 ρ =-1.2 κ=1.9 t = 2.0

exp((log(t) - α)exp(-ρ) - ρ) - exp(κ - ρ) (α - log(t))exp(κ - ρ) + (log(t) - α)exp((log(t) - α)exp(-ρ) - ρ) - 1 (log(t) - α)exp(κ - ρ) - exp(κ)SpecialFunctions.digamma(exp(κ))

source
LSurvival.jackknife — Method

Obtain jackknife (leave-one-out) estimates from a Aalen-Johansen risk curve (risk at end of follow-up) by refitting the model n times

Signatures

jackknife(m::M;kwargs...) where {M<:AJSurv}
source
LSurvival.jackknife — Method

Obtain jackknife (leave-one-out) estimates from a Kaplan-Meier survival curve (survival at end of follow-up) by refitting the model n times

Signatures

jackknife(m::M;kwargs...) where {M<:KMSurv}
using LSurvival, Random, StatsBase
 
 dat1 = (time = [1, 1, 6, 6, 8, 9], status = [1, 0, 1, 1, 0, 1], x = [1, 1, 1, 0, 0, 0])
 
@@ -503,7 +503,7 @@
 std(bs[:,1])
 stderror(m, type="jackknife")
 stderror(mc, type="jackknife")
-@assert jk == jkc
source
LSurvival.jackknife — Method

Obtain jackknife (leave-one-out) estimates from a Cox model by refitting the model n times

using LSurvival, Random, StatsBase
+@assert jk == jkc
source
LSurvival.jackknife — Method

Obtain jackknife (leave-one-out) estimates from a Cox model by refitting the model n times

using LSurvival, Random, StatsBase
 id, int, outt, data =
 LSurvival.dgm(MersenneTwister(112), 100, 10; afun = LSurvival.int_0)
 data[:, 1] = round.(data[:, 1], digits = 3)
@@ -544,7 +544,7 @@
 std(bs[:,1])
 stderror(m, type="jackknife")
 stderror(mc, type="jackknife")
-@assert jk == jkc
source
LSurvival.kaplan_meier — Method

Kaplan-Meier estimator for cumulative conditional risk

Signatures

StatsBase.fit!(m::T; kwargs...) where {T<:AbstractNPSurv}
+@assert jk == jkc
source
LSurvival.kaplan_meier — Method

Kaplan-Meier estimator for cumulative conditional risk

Signatures

StatsBase.fit!(m::T; kwargs...) where {T<:AbstractNPSurv}
 
 kaplan_meier(enter::AbstractVector, exit::AbstractVector, y::AbstractVector,
    ; kwargs...)

Keyword arguments

  • wts::Vector{<:Real} = similar(enter, 0); vector of case weights (or zero length vector) for each observation
  • id::Vector{<:AbstractLSurvivalID} = [ID(i) for i in eachindex(y)]; Vector of AbstractSurvID objects denoting observations that form a single unit (used in bootstrap and jackknife methods)
  • atol = 0.00000001; absolute tolerance for defining tied event times
  • censval = 0; value of the outcome to be considered a censored event
  • keepy = true; keep the outcome vector after fitting (may save memory with large datasets)
  • eps = 0.00000001; deprecated (replaced by atol)
using LSurvival
@@ -552,7 +552,7 @@
 z,x,t,d, event,wt = LSurvival.dgm_comprisk(MersenneTwister(1212), 1000);
 enter = zeros(length(t));
 m = fit(KMSurv, enter, t, d)
-mw = fit(KMSurv, enter, t, d, wts=wt)

or, equivalently:

kaplan_meier(enter, t, d, wts=wt)
source
LSurvival.lgh! — Method

dat1 = (time = [1, 1, 6, 6, 8, 9], status = [1, 0, 1, 1, 0, 1], x = [1, 1, 1, 0, 0, 0]) enter = zeros(length(dat1.time)) t = dat1.time d = dat1.status X = hcat(ones(length(dat1.x)), dat1.x) wt = ones(length(t))

dist = Exponential() P = PSParms(X[:,1:1], extraparms=length(dist)-1) P = PSParms(X, extraparms=length(dist)-1) P.B P.grad R = LSurvivalResp(dat1.time, dat1.status) # specification with ID only m = PSModel(R,P,dist)

λ=1 θ = rand(2) lgh!(m, θ) θ .+= inv(-m.P.hess) * m.P.grad * λ

lgh!(m, θ .+ [-.00, 0, 0]) lgh!(m, θ .+ [-.01, 0.05, 0.05])

m.P._LL

source
LSurvival.lgh_breslow! — Method

Update the partial likelihood, gradient and Hessian values from a Cox model fit (used during fitting, not generally useful for users).

Uses Breslow's partial likelihood.

Updates over all observations

Signature

lgh_breslow!(m::M, j, caseidx, risksetidx) where {M<:AbstractPH}
source
LSurvival.lgh_efron! — Method

Update the partial likelihood, gradient and Hessian values from a Cox model fit (used during fitting, not generally useful for users).

Uses Efron's partial likelihood.

Updates over all observations

Signature

lgh_efron!(m::M, j, caseidx, risksetidx) where {M<:AbstractPH}
source
LSurvival.loglik — Method

Log likelihood contribution for an observation in a parametric survival model

    d = m.d
+mw = fit(KMSurv, enter, t, d, wts=wt)

or, equivalently:

kaplan_meier(enter, t, d, wts=wt)
source
LSurvival.lgh! — Method

dat1 = (time = [1, 1, 6, 6, 8, 9], status = [1, 0, 1, 1, 0, 1], x = [1, 1, 1, 0, 0, 0]) enter = zeros(length(dat1.time)) t = dat1.time d = dat1.status X = hcat(ones(length(dat1.x)), dat1.x) wt = ones(length(t))

dist = Exponential() P = PSParms(X[:,1:1], extraparms=length(dist)-1) P = PSParms(X, extraparms=length(dist)-1) P.B P.grad R = LSurvivalResp(dat1.time, dat1.status) # specification with ID only m = PSModel(R,P,dist)

λ=1 θ = rand(2) lgh!(m, θ) θ .+= inv(-m.P.hess) * m.P.grad * λ

lgh!(m, θ .+ [-.00, 0, 0]) lgh!(m, θ .+ [-.01, 0.05, 0.05])

m.P._LL

source
LSurvival.lgh_breslow! — Method

Update the partial likelihood, gradient and Hessian values from a Cox model fit (used during fitting, not generally useful for users).

Uses Breslow's partial likelihood.

Updates over all observations

Signature

lgh_breslow!(m::M, j, caseidx, risksetidx) where {M<:AbstractPH}
source
LSurvival.lgh_efron! — Method

Update the partial likelihood, gradient and Hessian values from a Cox model fit (used during fitting, not generally useful for users).

Uses Efron's partial likelihood.

Updates over all observations

Signature

lgh_efron!(m::M, j, caseidx, risksetidx) where {M<:AbstractPH}
source
LSurvival.loglik — Method

Log likelihood contribution for an observation in a parametric survival model

    d = m.d
     i = 1
     enter = m.R.enter[i]
     exit = m.R.exit[i]
@@ -562,165 +562,165 @@
     θ=[1,0,.4]
     
     m.P._B
-
source
LSurvival.lpdf — Method

Log-likelihood calculation for Exponential regression: PDF

β = [-2, 1.2]
+
source
LSurvival.lpdf — Method

Log-likelihood calculation for Exponential regression: PDF

β = [-2, 1.2]
 x = [2,.1]
 ρ = -0.5
 t = 3.0
 α = dot(β,x)
 d = Exponential()
-lpdf_gradient(d, θ, t, x)
source
LSurvival.lpdf — Method

Log probability distribution function: Exponential distribution

    β = [-2, 1.2]
+lpdf_gradient(d, θ, t, x)
source
LSurvival.lpdf — Method

Log probability distribution function: Exponential distribution

    β = [-2, 1.2]
     x = [2,.1]
     ρ = -0.5
     t = 3.0
     α = dot(β,x)
     d = Exponential()
-    lpdf(d, t)
source
LSurvival.lpdf — Method

log probability distribution function, generalized gamma distribution

Location scale representation (Klein Moeschberger ch 12)

source
LSurvival.lpdf — Method

log probability distribution function, Gamma distribution

Location scale representation (Klein Moeschberger ch 12)

source
LSurvival.lpdf — Method

Log likelihood calculation for Lognormal regression: PDF

β = [-2, 1.2]
+    lpdf(d, t)
source
LSurvival.lpdf — Method

log probability distribution function, generalized gamma distribution

Location scale representation (Klein Moeschberger ch 12)

source
LSurvival.lpdf — Method

log probability distribution function, Gamma distribution

Location scale representation (Klein Moeschberger ch 12)

source
LSurvival.lpdf — Method

Log likelihood calculation for Lognormal regression: PDF

β = [-2, 1.2]
 x = [2,.1]
 ρ = -0.5
 t = 3.0
 α = dot(β,x)
 d = Lognormal()
-lpdf(d, vcat(θ,ρ), t, x)
source
LSurvival.lpdf — Method

log probability distribution function: Weibull distribution

Location scale representation (Klein Moeschberger ch 12)

source
LSurvival.lpdf — Method

Log-likelihood calculation for weibull regression: PDF

β = [-2, 1.2]
+lpdf(d, vcat(θ,ρ), t, x)
source
LSurvival.lpdf — Method

log probability distribution function: Weibull distribution

Location scale representation (Klein Moeschberger ch 12)

source
LSurvival.lpdf — Method

Log-likelihood calculation for weibull regression: PDF

β = [-2, 1.2]
 x = [2,.1]
 ρ = -0.5
 t = 3.0
 α = dot(β,x)
 d = Weibull()
-lpdf(d, vcat(θ,ρ), t, x)
source
LSurvival.lpdf — Method

Log probability distribution function: Weibull distribution

location scale representation (Klein Moeschberger ch 12)

α=0.1   # location
+lpdf(d, vcat(θ,ρ), t, x)
source
LSurvival.lpdf — Method

Log probability distribution function: Weibull distribution

location scale representation (Klein Moeschberger ch 12)

α=0.1   # location
 ρ=-1.2  # log(scale)
 time=2
-lpdf(Weibull(α, ρ), time)
source
LSurvival.lpdf_gradient — Method

Gradient calculation for Exponential regression: PDF

β = [-2, 1.2]
 x = [2,.1]
 ρ = -0.5
 t = 3.0
 α = dot(β,x)
 d = Exponential()
-lpdf_gradient(d, θ, t, x)
source
LSurvival.lpdf_gradient — Method

Gradient calculation for Lognormal regression: PDF

β = [-2, 1.2]
 x = [2,.1]
 ρ = -0.5
 t = 3.0
 α = dot(β,x)
 d = Lognormal()
-lpdf_gradient(d, vcat(θ,ρ), t, x)
source
LSurvival.lpdf_gradient — Method

Gradient calculation for weibull regression: PDF

β = [-2, 1.2]
 x = [2,.1]
 ρ = -0.5
 t = 3.0
 #α = dot(β,x)
 d = Weibull()
-lpdf_gradient(d, vcat(θ,ρ), t, x)
source
LSurvival.lpdf_hessian — Method

Hessian calculation for Exponential regression: PDF

β = [-2, 1.2]
+lpdf_gradient(d, vcat(θ,ρ), t, x)
source
LSurvival.lpdf_hessian — Method

Hessian calculation for Exponential regression: PDF

β = [-2, 1.2]
 x = [2,.1]
 ρ = -0.5
 t = 3.0
 α = dot(β,x)
 d = Exponential()
-lpdf_hessian(d, θ, t, x)
source
LSurvival.lpdf_hessian — Method

Hessian calculation for Weibull distribution: PDF

β = [-2, 1.2]
 x = [2,.1]
 ρ = -0.5
 t = 3.0
 α = dot(β,x)
 d = Exponential(α)
-lpdf_hessian(d, t)
source
LSurvival.lpdf_hessian — Method

Hessian calculation for Log-normal regression: PDF

β = [-2, 1.2]
 x = [2,.1]
 ρ = -0.5
 t = 3.0
 α = dot(β,x)
 d = Lognormal()
-lpdf_hessian(d, vcat(θ,ρ), t, x)
source
LSurvival.lpdf_hessian — Method

Hessian calculation for Log-normal distribution: PDF

β = [-2, 1.2]
+lpdf_hessian(d, vcat(θ,ρ), t, x)
source
LSurvival.lpdf_hessian — Method

Hessian calculation for Log-normal distribution: PDF

β = [-2, 1.2]
 x = [2,.1]
 ρ = -0.5
 t = 3.0
 α = dot(β,x)
 d = Lognormal(α, ρ)
-lpdf_hessian(d, t)
source
LSurvival.lpdf_hessian — Method

Hessian calculation for weibull regression: PDF

β = [-2, 1.2]
 x = [2,.1]
 ρ = -0.5
 t = 3.0
 α = dot(β,x)
 d = Weibull()
-lpdf_hessian(d, vcat(θ,ρ), t, x)
source
LSurvival.lsurv — Method

Log-likelihood calculation for Exponential regression: Survival

β = [-2, 1.2]
+lpdf_hessian(d, vcat(θ,ρ), t, x)
source
LSurvival.lsurv — Method

Log-likelihood calculation for Exponential regression: Survival

β = [-2, 1.2]
 x = [2,.1]
 ρ = -0.5
 t = 3.0
 α = dot(β,x)
 d = Exponential()
-lsurv(d, θ, t, x)
source
LSurvival.lsurv — Method

Log survival function: Exponential distribution

    β = [-2, 1.2]
+lsurv(d, θ, t, x)
source
LSurvival.lsurv — Method

Log survival function: Exponential distribution

    β = [-2, 1.2]
     x = [2,.1]
     ρ = -0.5
     t = 3.0
     α = dot(β,x)
     d = Exponential()
-    lsurv(d, t)
source
LSurvival.lsurv — Method

log probability distribution function, generalized gamma distribution

Location scale representation (Klein Moeschberger ch 12)

source
LSurvival.lsurv — Method

log probability distribution function, Gamma distribution

Location scale representation (Klein Moeschberger ch 12)

source
LSurvival.lsurv — Method

Log likelihood calculation for Log-normal regression: Survival

β = [-2, 1.2]
+    lsurv(d, t)
source
LSurvival.lsurv — Method

log probability distribution function, generalized gamma distribution

Location scale representation (Klein Moeschberger ch 12)

source
LSurvival.lsurv — Method

log probability distribution function, Gamma distribution

Location scale representation (Klein Moeschberger ch 12)

source
LSurvival.lsurv — Method

Log likelihood calculation for Log-normal regression: Survival

β = [-2, 1.2]
 x = [2,.1]
 ρ = -0.5
 t = 3.0
 α = dot(β,x)
 d = Lognormal()
-lsurv(d, vcat(θ,ρ), t, x)
source
LSurvival.lsurv — Method

log probability distribution function: Weibull distribution

Location scale representation (Klein Moeschberger ch 12)

source
LSurvival.lsurv — Method

Log-likelihood calculation for weibull regression: Survival

β = [-2, 1.2]
+lsurv(d, vcat(θ,ρ), t, x)
source
LSurvival.lsurv — Method

log probability distribution function: Weibull distribution

Location scale representation (Klein Moeschberger ch 12)

source
LSurvival.lsurv — Method

Log-likelihood calculation for weibull regression: Survival

β = [-2, 1.2]
 x = [2,.1]
 ρ = -0.5
 t = 3.0
 α = dot(β,x)
 d = Weibull()
-lsurv(d, vcat(θ,ρ), t, x)
source
LSurvival.lsurv — Method

Log survival distribution function: Weibull distribution

location, log(scale) representation (Klein Moeschberger ch 12)

α=0.1   # location
+lsurv(d, vcat(θ,ρ), t, x)
source
LSurvival.lsurv — Method

Log survival distribution function: Weibull distribution

location, log(scale) representation (Klein Moeschberger ch 12)

α=0.1   # location
 ρ=-1.2  # log(scale)
 time=2
-lsurv(Weibull(α, ρ), time)
source
LSurvival.lsurv_gradient — Method

Gradient calculation for Exponential regression: Survival

β = [-2, 1.2]
 x = [2,.1]
 ρ = -0.5
 t = 3.0
 α = dot(β,x)
 d = Exponential()
-lsurv_gradient(d, θ, t, x)
source
LSurvival.lsurv_gradient — Method

Gradient calculation for Log-normal regression: Survival

β = [-2, 1.2]
 x = [2,.1]
 ρ = -0.5
 t = 3.0
 α = dot(β,x)
 d = Lognormal()
-lsurv_gradient(d, vcat(θ,ρ), t, x)
source
LSurvival.lsurv_gradient — Method

Gradient calculation for weibull regression: Survival

β = [-2, 1.2]
 x = [2,.1]
 ρ = -0.5
 t = 3.0
 α = dot(β,x)
 d = Weibull()
-lsurv_gradient(d, vcat(θ,ρ), t, x)
source
LSurvival.lsurv_hessian — Method

Hessian calculation for Exponential regression: Survival

β = [-2, 1.2]
+lsurv_gradient(d, vcat(θ,ρ), t, x)
source
LSurvival.lsurv_hessian — Method

Hessian calculation for Exponential regression: Survival

β = [-2, 1.2]
 x = [2,.1]
 ρ = -0.5
 t = 3.0
 α = dot(β,x)
 d = Exponential()
-lsurv_hessian(d, θ, t, x)
source
LSurvival.lsurv_hessian — Method

Hessian calculation for Exponential distribution: Survival

β = [-2, 1.2]
 x = [2,.1]
 ρ = -0.5
 t = 3.0
 α = dot(β,x)
 d = Exponential(α)
-lsurv_hessian(d, t)
source
LSurvival.lsurv_hessian — Method

Hessian calculation for Log-normal regression: Survival

β = [-2, 1.2]
 x = [2,.1]
 ρ = -0.5
 t = 3.0
 α = dot(β,x)
 d = Lognormal()
-lsurv_hessian(d, vcat(θ,ρ), t, x)
source
LSurvival.lsurv_hessian — Method

Hessian calculation for Log-normal distribution: Survival

β = [-2, 1.2]
+lsurv_hessian(d, vcat(θ,ρ), t, x)
source
LSurvival.lsurv_hessian — Method

Hessian calculation for Log-normal distribution: Survival

β = [-2, 1.2]
 x = [2,.1]
 ρ = -0.5
 t = 3.0
 α = dot(β,x)
 d = Lognormal(α, ρ)
-lsurv_hessian(d, t)
source
LSurvival.lsurv_hessian — Method

Hessian calculation for weibull regression: Survival

β = [-2, 1.2]
 x = [2,.1]
 ρ = -0.5
 t = 3.0
 α = dot(β,x)
 d = Weibull()
-lsurv_hessian(d, vcat(θ,ρ), t, x)
source
LSurvival.lsurv_hessian — Method

Hessian calculation for Weibull distribution: Survival

β = [-2, 1.2]
+lsurv_hessian(d, vcat(θ,ρ), t, x)
source
LSurvival.lsurv_hessian — Method

Hessian calculation for Weibull distribution: Survival

β = [-2, 1.2]
 x = [2,.1]
 ρ = -0.5
 t = 3.0
 α = dot(β,x)
 d = Weibull(α, ρ)
-lsurv_hessian(d, t)
source
LSurvival.qstdnorm — Method

quantile function for a standard normal distribution depends on SpecialFunctions https://en.wikipedia.org/wiki/Normal_distribution

Source code, example:
    qstdnorm(p) = sqrt(2) * SpecialFunctions.erfinv(2.0 * p - 1.0)
+lsurv_hessian(d, t)
source
LSurvival.qstdnorm — Method

quantile function for a standard normal distribution depends on SpecialFunctions https://en.wikipedia.org/wiki/Normal_distribution

Source code, example:
    qstdnorm(p) = sqrt(2) * SpecialFunctions.erfinv(2.0 * p - 1.0)
     
-    qstdnorm(.975)
source
LSurvival.qweibull — Method

Quantile function for the Weibull distribution

\[F(t) = egin{cases} 1 - e^{{-t/ρ}^{α}}&, t ≥ 0\ 0&, t < 0 nd{cases} @@ -728,13 +728,13 @@ # cross reference the approach in the Distributions package quantile(Distributions.Weibull(.75, 1.1), .3) -LSurvival.qweibull(0.3, .75, 1.1)

source
LSurvival.randweibull — Method

Random draw from Weibull distribution

lightweight function used for simulation

Note that there is no checking that parameters α,ρ are positively bound, and errors will be given if this is not the case

Signatures:

randweibull(rng::MersenneTwister,α::Real,ρ::Real)
+LSurvival.qweibull(0.3, .75, 1.1)
source
LSurvival.randweibull — Method

Random draw from Weibull distribution

lightweight function used for simulation

Note that there is no checking that parameters α,ρ are positively bound, and errors will be given if this is not the case

Signatures:

randweibull(rng::MersenneTwister,α::Real,ρ::Real)
 randweibull(α::Real,ρ::Real)

Source code, example:

randweibull(rng, α, ρ) = qweibull(rand(rng), α, ρ)
 randweibull(α, ρ) = randweibull(MersenneTwister(), α, ρ)
 
 # cross reference the approach in the Distributions package
 rand(Distributions.Weibull(.75, 1.1))
-randweibull(.75, 1.1)
source
LSurvival.risk_from_coxphmodels — Method

Survival curve estimation using multiple cox models

Signatures

  risk_from_coxphmodels(fitlist::Vector{T}, args...; kwargs...) where {T<:PHModel}
+randweibull(.75, 1.1)
source
LSurvival.risk_from_coxphmodels — Method

Survival curve estimation using multiple cox models

Signatures

  risk_from_coxphmodels(fitlist::Vector{T}, args...; kwargs...) where {T<:PHModel}
 
   fit(::Type{M}, fitlist::Vector{T}, ; fitargs...) where {M<:PHSurv,T<:PHModel}

Optional keywords

  • coef_vectors = nothing(default) or vector of coefficient vectors from the cox models [will default to the coefficients from fitlist models]
  • pred_profile = nothing(default) or vector of specific predictor values of the same length as the coef_vectors[1]
 using LSurvival
  using Random
@@ -754,19 +754,19 @@
  mnz = sum(z)/length(z)
  # equivalent
  fit(PHSurv, [ft1,ft2], pred_profile=[mnx,mnz])
- risk_from_coxphmodels([ft1,ft2], pred_profile=[mnx,mnz])
source
LSurvival.survivaldata — Method

Loading example survival analysis datasets

using LSurvival, Plots # note Plots does not install by default
+ risk_from_coxphmodels([ft1,ft2], pred_profile=[mnx,mnz])
source
LSurvival.survivaldata — Method

Loading example survival analysis datasets

using LSurvival, Plots # note Plots does not install by default
 heartdata, heartmeta = survivaldata("heart")
 ft = coxph(@formula(Surv(start, stop, event)~surgery), heartdata);
 # plot baseline cumulative hazard
 basehazplot(ft)
 # plot Schoenfeld residuals
-coxdx(ft)
source
RecipesBase.apply_recipe — Method

Recipe for aalen-johansen risk curve

    using Plots, LSurvival
     res = z, x, outt, d, event, weights = LSurvival.dgm_comprisk(MersenneTwister(123123), 100)
     int = zeros(length(d)) # no late entry
     
         c = fit(AJSurv, int, outt, event)
         #risk2 = aalen_johansen(int, outt, event)
-        plot(c)
source
RecipesBase.apply_recipe — Method

Plotting a kaplan meier curve

    using Plots, LSurvival
 dat4 = (
     id = [1, 1, 2, 2, 2, 3, 4, 5, 5, 6],
     enter = [1, 2, 5, 4, 6, 7, 3, 6, 8, 0],
@@ -778,7 +778,7 @@
 )
 R = LSurvivalResp(dat4.enter, dat4.exit, dat4.status)
     k = kaplan_meier(dat4.enter, dat4.exit, dat4.status)
-    plot(k)
source
RecipesBase.apply_recipe — Method
function name(::Type{T}) where {T}
     #https://stackoverflow.com/questions/70043313/get-simple-name-of-type-in-julia
     isempty(T.parameters) ? T : T.name.wrapper
 end
@@ -807,7 +807,7 @@
 # Cumulative incidence/risk function
 aftdist(fte, type="risk", label="X=0")
 aftdist!(fte, type="risk", covlevels=[1.0, 2.0], color="red", label="X=1")
-
source
RecipesBase.apply_recipe — Method

Plotting baseline hazard for a Cox model

using Plots, LSurvival
 dat2 = (
     enter = [1, 2, 5, 2, 1, 7, 3, 4, 8, 8],
     exit = [2, 3, 6, 7, 8, 9, 9, 9, 14, 17],
@@ -818,7 +818,7 @@
 ftb = coxph(@formula(Surv(enter, exit, status)~x), dat2, ties="breslow", maxiter=0)
 
 plot(fte, label="Efron")
-plot!(ftb, label="Breslow")
source
RecipesBase.apply_recipe — Method

 using Plots, LSurvival
 dat2 = (
     enter = [1, 2, 5, 2, 1, 7, 3, 4, 8, 8],
@@ -829,7 +829,7 @@
 fte = coxph(@formula(Surv(enter, exit, status)~x), dat2)
 
 coxdx(fte)
-
source
RecipesBase.apply_recipe — Method
using Plots, LSurvival
 dat2 = (
     enter = [1, 2, 5, 2, 1, 7, 3, 4, 8, 8],
     exit = [2, 3, 6, 7, 8, 9, 9, 9, 14, 17],
@@ -840,7 +840,7 @@
 
 coxinfluence(fte, type="jackknife", par=1)
 coxinfluence!(fte, type="dfbeta", color=:red, par=1)
-
source
RecipesBase.apply_recipe — Method

Plotting baseline hazard for a Cox model

using Plots, LSurvival
 dat4 = (
     id = [1, 1, 2, 2, 2, 3, 4, 5, 5, 6],
     enter = [1, 2, 5, 4, 6, 7, 3, 6, 8, 0],
@@ -853,7 +853,7 @@
 
 k = kaplan_meier(dat4.enter, dat4.exit, dat4.status)
 
-lognlogplot(k)
source
RecipesBase.apply_recipe — Method

Plotting LSurvivalResp objects (outcomes in cox models, kaplan meier curves, parametric survival models)

Recipe for plotting time-to-event outcomes

using Plots, LSurvival
+lognlogplot(k)
source
RecipesBase.apply_recipe — Method

Plotting LSurvivalResp objects (outcomes in cox models, kaplan meier curves, parametric survival models)

Recipe for plotting time-to-event outcomes

using Plots, LSurvival
 
 dat4 = (
     id = [1, 1, 2, 2, 2, 3, 4, 5, 5, 6],
@@ -865,7 +865,7 @@
     w = [0, 0, 0, 0, 0, 1, 1, 1, 1, 0],
 )
 R = LSurvivalResp(dat4.enter, dat4.exit, dat4.status)
-plot([[R.enter[i], R.exit[i]] for i in eachindex(R.enter)], [[i, i] for i in values(R.id)])
source
RecipesBase.apply_recipe — Method

Recipe for cox-model based risk curves

    using Plots, LSurvival, Random, StatsBase
+plot([[R.enter[i], R.exit[i]] for i in eachindex(R.enter)], [[i, i] for i in values(R.id)])
source
RecipesBase.apply_recipe — Method

Recipe for cox-model based risk curves

    using Plots, LSurvival, Random, StatsBase
     res = z, x, outt, d, event, wts = LSurvival.dgm_comprisk(MersenneTwister(123123), 100)
     X = hcat(z, x)
     int = zeros(length(d)) # no late entry
@@ -873,13 +873,13 @@
     ft2 = fit(PHModel, X, int, outt, d .* (event .== 2), wts=wts)
     c = risk_from_coxphmodels([ft1, ft2], pred_profile = mean(X, dims=1))
     
-    plot(c)
source
StatsAPI.confint — Method

Greenwood's formula for variance and confidence intervals of a Aalen-Johansen risk function

Signatures:

StatsBase.stderror(m::AJSurv)
+    plot(c)
source
StatsAPI.confint — Method

Greenwood's formula for variance and confidence intervals of a Aalen-Johansen risk function

Signatures:

StatsBase.stderror(m::AJSurv)
 
 StatsBase.confint(m:AJSurv; level=0.95, method="normal")

Keyword arguments

  • method
    • "normal" normality-based confidence intervals
    • "lognlog" log(-log(S(t))) based confidence intervals
res = z, x, outt, d, event, wts = LSurvival.dgm_comprisk(MersenneTwister(123123), 100)
 int = zeros(length(d)) # no late entry
 m = fit(AJSurv, int, outt, event)
 stderror(m)
-confint(m, level=0.95)
source
StatsAPI.confint — Method

Greenwood's formula for variance and confidence intervals of a Kaplan-Meier survival curve

Signatures:

StatsBase.stderror(m::KMSurv)
+confint(m, level=0.95)
source
StatsAPI.confint — Method

Greenwood's formula for variance and confidence intervals of a Kaplan-Meier survival curve

Signatures:

StatsBase.stderror(m::KMSurv)
 
 StatsBase.confint(m:KMSurv; level=0.95, method="normal")

Keyword arguments

method:

  • "normal" normality-based confidence intervals
  • "lognlog" log(-log(S(t))) based confidence intervals
using LSurvival
 using Random
@@ -889,7 +889,7 @@
 mw = fit(KMSurv, enter, t, d, wts=wt)
 stderror(m)
 confint(m, method="normal")
-confint(m, method="lognlog") # log-log transformation
source
StatsAPI.confint — Method
using LSurvival
  dat1= (
    time = [1,1,6,6,8,9],
    status = [1,0,1,1,0,1],
@@ -914,7 +914,7 @@
  confint(ft2)
 
  # robust variance
- confint(ft2, type="robust")
source
StatsAPI.fit! — Method

Kaplan-Meier estimator for cumulative conditional risk

Signatures

StatsBase.fit!(m::T; kwargs...) where {T<:AbstractNPSurv}
+ confint(ft2, type="robust")
source
StatsAPI.fit! — Method

Kaplan-Meier estimator for cumulative conditional risk

Signatures

StatsBase.fit!(m::T; kwargs...) where {T<:AbstractNPSurv}
 
 kaplan_meier(enter::AbstractVector, exit::AbstractVector, y::AbstractVector,
    ; kwargs...)

Keyword arguments

  • wts::Vector{<:Real} = similar(enter, 0); vector of case weights (or zero length vector) for each observation
  • id::Vector{<:AbstractLSurvivalID} = [ID(i) for i in eachindex(y)]; Vector of AbstractSurvID objects denoting observations that form a single unit (used in bootstrap and jackknife methods)
  • atol = 0.00000001; absolute tolerance for defining tied event times
  • censval = 0; value of the outcome to be considered a censored event
  • keepy = true; keep the outcome vector after fitting (may save memory with large datasets)
  • eps = 0.00000001; deprecated (replaced by atol)
using LSurvival
@@ -952,7 +952,7 @@
  mnz = sum(z)/length(z)
  # equivalent
  fit(PHSurv, [ft1,ft2], pred_profile=[mnx,mnz])
- risk_from_coxphmodels([ft1,ft2], pred_profile=[mnx,mnz])
source
StatsAPI.fit — Method

Survival curve estimation using multiple cox models

Signatures

  risk_from_coxphmodels(fitlist::Vector{T}, args...; kwargs...) where {T<:PHModel}
+ risk_from_coxphmodels([ft1,ft2], pred_profile=[mnx,mnz])
source
StatsAPI.fit — Method

Survival curve estimation using multiple cox models

Signatures

  risk_from_coxphmodels(fitlist::Vector{T}, args...; kwargs...) where {T<:PHModel}
 
   fit(::Type{M}, fitlist::Vector{T}, ; fitargs...) where {M<:PHSurv,T<:PHModel}

Optional keywords

  • coef_vectors = nothing(default) or vector of coefficient vectors from the cox models [will default to the coefficients from fitlist models]
  • pred_profile = nothing(default) or vector of specific predictor values of the same length as the coef_vectors[1]
 using LSurvival
  using Random
@@ -972,7 +972,7 @@
  mnz = sum(z)/length(z)
  # equivalent
  fit(PHSurv, [ft1,ft2], pred_profile=[mnx,mnz])
- risk_from_coxphmodels([ft1,ft2], pred_profile=[mnx,mnz])
source
StatsAPI.fit — Method

Fit method for AbstractPH objects (Cox models)

Keyword arguments (used here, and passed on to internal structs)

  • ties "breslow" or "efron" (default)

  • wts observation weights

  • ties "breslow" or "efron" (default)

  • offset not currently used at all

  • fitargs arguments passed to other structs, which include

    • id cluster or individual level ID (defaults to a unique value for each row of data) see note below on ID
    • contrasts StatsModel style contrasts (dicts) that can be used for variable transformations/indicator variable creation (e.g. https://juliastats.org/StatsModels.jl/stable/contrasts/)

  • Arguments passed onto fitting routine:

    • eps (default: Float64 = 1e-9) tolerance for declaring convergence. Model is determined to be converged when relative change in log-partial likelihood is < eps .
    • getbasehaz (default: true): estimate baseline hazard
    • start (default: nothing) nothing, or vector of floats corresponding to initial values for parameters. Note that this defaults to a vector of zeros when set to nothing, and setting to other values invalidates some of the test statistics reported by default with coxph.
    • keepx (default: true) logical. Keep design matrix in AbstractPH object output (set to false for slight computational gain).
    • keepy (default: true)logical. Keep outcome in AbstractPH object output (set to false for slight computational gain).
    • bootstrap_sample (default: false) Fit the model to a bootstrap sample of the data (not generally used by end-users, but provides some convenience toward bootstrap variance estimation).
    • bootstrap_rng (default: Random.MersenneTwister()) Random number seed used when drawing a bootstrap sample of the data (not generally used by end-users, but provides some convenience toward bootstrap variance estimation).

    Signatures

  fit(::Type{M},
+ risk_from_coxphmodels([ft1,ft2], pred_profile=[mnx,mnz])
source
StatsAPI.fit — Method

Fit method for AbstractPH objects (Cox models)

Keyword arguments (used here, and passed on to internal structs)

  • ties "breslow" or "efron" (default)

  • wts observation weights

  • ties "breslow" or "efron" (default)

  • offset not currently used at all

  • fitargs arguments passed to other structs, which include

    • id cluster or individual level ID (defaults to a unique value for each row of data) see note below on ID
    • contrasts StatsModel style contrasts (dicts) that can be used for variable transformations/indicator variable creation (e.g. https://juliastats.org/StatsModels.jl/stable/contrasts/)

  • Arguments passed onto fitting routine:

    • eps (default: Float64 = 1e-9) tolerance for declaring convergence. Model is determined to be converged when relative change in log-partial likelihood is < eps .
    • getbasehaz (default: true): estimate baseline hazard
    • start (default: nothing) nothing, or vector of floats corresponding to initial values for parameters. Note that this defaults to a vector of zeros when set to nothing, and setting to other values invalidates some of the test statistics reported by default with coxph.
    • keepx (default: true) logical. Keep design matrix in AbstractPH object output (set to false for slight computational gain).
    • keepy (default: true)logical. Keep outcome in AbstractPH object output (set to false for slight computational gain).
    • bootstrap_sample (default: false) Fit the model to a bootstrap sample of the data (not generally used by end-users, but provides some convenience toward bootstrap variance estimation).
    • bootstrap_rng (default: Random.MersenneTwister()) Random number seed used when drawing a bootstrap sample of the data (not generally used by end-users, but provides some convenience toward bootstrap variance estimation).

    Signatures

  fit(::Type{M},
   X::AbstractMatrix,#{<:FP},
   enter::AbstractVector{<:Real},
   exit::AbstractVector{<:Real},
@@ -1041,7 +1041,7 @@
  stderror(ft2w, type="robust")      # robust variance (INCORRECT)
  
  ft2w # the coefficient table now shows incorrect confidence intervals and test statistics
-
source
StatsAPI.fit — Method

Aalen-Johansen estimator for cumulative cause-specific risk (in the presence of competing events)

Signatures

 StatsBase.fit!(m::T; kwargs...) where {T<:AbstractNPSurv}
+
source
StatsAPI.fit — Method

Aalen-Johansen estimator for cumulative cause-specific risk (in the presence of competing events)

Signatures

 StatsBase.fit!(m::T; kwargs...) where {T<:AbstractNPSurv}
 
  aalen_johansen(enter::AbstractVector, exit::AbstractVector, y::AbstractVector,
    ; kwargs...)
@@ -1051,7 +1051,7 @@
 enter = zeros(length(t));
    # event variable is coded 0[referent],1,2
 m = fit(AJSurv, enter, t, event)
-mw = fit(AJSurv, enter, t, event, wts=wt)

or, equivalently:

aalen_johansen(enter, t, event, wts=wt)
source
StatsAPI.fit — Method

Kaplan-Meier estimator for cumulative conditional risk

Signatures

StatsBase.fit!(m::T; kwargs...) where {T<:AbstractNPSurv}
+mw = fit(AJSurv, enter, t, event, wts=wt)

or, equivalently:

aalen_johansen(enter, t, event, wts=wt)
source
StatsAPI.fit — Method

Kaplan-Meier estimator for cumulative conditional risk

Signatures

StatsBase.fit!(m::T; kwargs...) where {T<:AbstractNPSurv}
 
 kaplan_meier(enter::AbstractVector, exit::AbstractVector, y::AbstractVector,
    ; kwargs...)

Keyword arguments

  • wts::Vector{<:Real} = similar(enter, 0); vector of case weights (or zero length vector) for each observation
  • id::Vector{<:AbstractLSurvivalID} = [ID(i) for i in eachindex(y)]; Vector of AbstractSurvID objects denoting observations that form a single unit (used in bootstrap and jackknife methods)
  • atol = 0.00000001; absolute tolerance for defining tied event times
  • censval = 0; value of the outcome to be considered a censored event
  • keepy = true; keep the outcome vector after fitting (may save memory with large datasets)
  • eps = 0.00000001; deprecated (replaced by atol)
using LSurvival
@@ -1059,7 +1059,7 @@
 z,x,t,d, event,wt = LSurvival.dgm_comprisk(MersenneTwister(1212), 1000);
 enter = zeros(length(t));
 m = fit(KMSurv, enter, t, d)
-mw = fit(KMSurv, enter, t, d, wts=wt)

or, equivalently:

kaplan_meier(enter, t, d, wts=wt)
source
StatsAPI.loglikelihood — Method

Maximum log partial likelihood for a fitted AbstractPH model Efron or Breslow (depending on the ties` parameter)

source
StatsAPI.nullloglikelihood — Method

Null log-partial likelihood for a fitted AbstractPH model Efron or Breslow (depending on the ties` parameter)

Note: this is just the log partial likelihood at the initial values of the model, which default to 0. If initial values are non-null, then this function no longer validly returns the null log-partial likelihood.

source
StatsAPI.nullloglikelihood — Method

Null log-partial likelihood for a fitted PSModel model

Note: this is just the log partial likelihood at the initial values of the model, which default to 0. If initial values are non-null, then this function no longer validly returns the null log-partial likelihood.

source
StatsAPI.residuals — Method

####################################################################

Cox proportional hazards model residuals:

Signature

  residuals(m::M; type = "martingale") where {M<:PHModel}

where type is one of

  • martingale

  • schoenfeld

  • score

  • dfbeta

  • jackknife

  • dfbetas (scaled dfbeta)

  • scaled_schoenfeld or schoenfelds (scaled Schoenfeld)

    Residuals from the residuals function are designed to exactly emulate those from the survival package in R. Currently, they are validated for single observation data (e.g. one data row per individual).

    ####################################################################

    Martingale residuals: Observed versus expected

  # example from https://cran.r-project.org/web/packages/survival/vignettes/validate.pdf
+mw = fit(KMSurv, enter, t, d, wts=wt)

or, equivalently:

kaplan_meier(enter, t, d, wts=wt)
source
StatsAPI.loglikelihood — Method

Maximum log partial likelihood for a fitted AbstractPH model Efron or Breslow (depending on the ties` parameter)

source
StatsAPI.nullloglikelihood — Method

Null log-partial likelihood for a fitted AbstractPH model Efron or Breslow (depending on the ties` parameter)

Note: this is just the log partial likelihood at the initial values of the model, which default to 0. If initial values are non-null, then this function no longer validly returns the null log-partial likelihood.

source
StatsAPI.nullloglikelihood — Method

Null log-partial likelihood for a fitted PSModel model

Note: this is just the log partial likelihood at the initial values of the model, which default to 0. If initial values are non-null, then this function no longer validly returns the null log-partial likelihood.

source
StatsAPI.residuals — Method

####################################################################

Cox proportional hazards model residuals:

Signature

  residuals(m::M; type = "martingale") where {M<:PHModel}

where type is one of

  • martingale

  • schoenfeld

  • score

  • dfbeta

  • jackknife

  • dfbetas (scaled dfbeta)

  • scaled_schoenfeld or schoenfelds (scaled Schoenfeld)

    Residuals from the residuals function are designed to exactly emulate those from the survival package in R. Currently, they are validated for single observation data (e.g. one data row per individual).

    ####################################################################

    Martingale residuals: Observed versus expected

  # example from https://cran.r-project.org/web/packages/survival/vignettes/validate.pdf
   # by Terry Therneau
 
   dat1 = (
@@ -1146,13 +1146,13 @@
   #jackknife(ft)
   residuals(ft, type="jackknife")
 
-
source
StatsAPI.stderror — Method

Greenwood's formula for variance and confidence intervals of a Aalen-Johansen risk function

Signatures:

StatsBase.stderror(m::AJSurv)
+
source
StatsAPI.stderror — Method

Greenwood's formula for variance and confidence intervals of a Aalen-Johansen risk function

Signatures:

StatsBase.stderror(m::AJSurv)
 
 StatsBase.confint(m:AJSurv; level=0.95, method="normal")

Keyword arguments

  • method
    • "normal" normality-based confidence intervals
    • "lognlog" log(-log(S(t))) based confidence intervals
res = z, x, outt, d, event, wts = LSurvival.dgm_comprisk(MersenneTwister(123123), 100)
 int = zeros(length(d)) # no late entry
 m = fit(AJSurv, int, outt, event)
 stderror(m)
-confint(m, level=0.95)
source
StatsAPI.stderror — Method

Greenwood's formula for variance and confidence intervals of a Kaplan-Meier survival curve

Signatures:

StatsBase.stderror(m::KMSurv)
+confint(m, level=0.95)
source
StatsAPI.stderror — Method

Greenwood's formula for variance and confidence intervals of a Kaplan-Meier survival curve

Signatures:

StatsBase.stderror(m::KMSurv)
 
 StatsBase.confint(m:KMSurv; level=0.95, method="normal")

Keyword arguments

method:

  • "normal" normality-based confidence intervals
  • "lognlog" log(-log(S(t))) based confidence intervals
using LSurvival
 using Random
@@ -1162,7 +1162,7 @@
 mw = fit(KMSurv, enter, t, d, wts=wt)
 stderror(m)
 confint(m, method="normal")
-confint(m, method="lognlog") # log-log transformation
source
StatsAPI.vcov — Method

Covariance matrix for Cox proportional hazards models

Keyword arguments

  • type nothing or "robust": determines whether model based or robust (dfbeta based) variance is returned.

    See ?residuals for info on dfbeta residuals

using LSurvival
+confint(m, method="lognlog") # log-log transformation
source
StatsAPI.vcov — Method

Covariance matrix for Cox proportional hazards models

Keyword arguments

  • type nothing or "robust": determines whether model based or robust (dfbeta based) variance is returned.

    See ?residuals for info on dfbeta residuals

using LSurvival
 dat1 = (
   time = [1,1,6,6,8,9],
   status = [1,0,1,1,0,1],
@@ -1195,4 +1195,4 @@
 vcov(ft2w, type="robust")    # robust variance (INCORRECT)
 nobs(ft2w)
 
-ft2w
source

Implementation details and further help

+ft2wsource

Implementation details and further help

diff --git a/dev/nonparametric/index.html b/dev/nonparametric/index.html index 3195a14..012e900 100644 --- a/dev/nonparametric/index.html +++ b/dev/nonparametric/index.html @@ -78,4 +78,4 @@ Number of events (j=1.0): 46.8371 Number of events (j=2.0): 58.4336 Number of unique event times: 106

Plotting marginal cause-specific risks

plot(res_aj)
-savefig("aj.svg")

Aalen-Johansen

+savefig("aj.svg")

Aalen-Johansen

diff --git a/dev/parametric/index.html b/dev/parametric/index.html index ae97a48..cc0db9c 100644 --- a/dev/parametric/index.html +++ b/dev/parametric/index.html @@ -180,4 +180,4 @@ Log-likelihood (full): -745.685 Log-likelihood (Intercept only): -963.049 LRT p-value (X^2=434.73, df=1): 0 -Solver iterations: 66

Notes: analytic gradients and Hessian matrixes are not available for this distribution, so the solver uses finite differencing, which can make this model considerably slower to fit than alternative models.

Log-logistic

In progress

Gompertz

In progress

+Solver iterations: 66

Notes: analytic gradients and Hessian matrixes are not available for this distribution, so the solver uses finite differencing, which can make this model considerably slower to fit than alternative models.

Log-logistic

In progress

Gompertz

In progress

diff --git a/dev/search/index.html b/dev/search/index.html index 94c6003..03eda3d 100644 --- a/dev/search/index.html +++ b/dev/search/index.html @@ -1,2 +1,2 @@ -- · LSurvival: survival analysis for left-truncated, right-censored outcomes
GitHub
+- · LSurvival: survival analysis for left-truncated, right-censored outcomes
GitHub