Gamma distribution time-to-event model

How-to
1

Hi - for fitting a TTE model in R with gamma distribution, I would do:

flexsurvreg(Surv(time, event) ~ cov1 + cov2 , data=df, dist="gamma")

Could someone please let me know how to do this in Pumas? Thanks.

2

If you follow the Time to Event tutorials and implement the Gamma distribution in the @pre and @vars blocks that should be enough.

3

Thank you, Jose. I apologize, I should’ve been clearer in my ask. My question is stemming from the fact there is no analytical form to the hazard function in this case (like Exponential, Weibull, and Gompertz as you have shown in the tutorials). Can you please let me know how the model code would look like if we assume the failure time to follow a gamma distribution?

4

@storopoli I was able to perform some simulations using the gamma function available in the SpecialFunctions package but I cannot get the fit to work. I am attaching the code for simulations, fit, and error from fit at the bottom. Please let me know if it is possible to make the fitting work. Thanks.

Simulation code:

using SpecialFunctions

gamma_model = @model begin
    @param begin
        scale1  ∈ RealDomain(lower=0)
        shape1  ∈ RealDomain(lower=0)
    end
    @pre begin
        scale_typ  = scale1
        shape_typ  = shape1

        pdf_gamma = (t^(shape_typ-1.) * exp(-t/scale_typ)) / (gamma(shape_typ) * scale_typ^shape_typ)
        cdf_gamma = gamma(shape_typ, t/scale_typ) - gamma(shape_typ)
        haz       = pdf_gamma / (1. - cdf_gamma)
    end
    @vars begin
        λ   = haz
    end
    @dynamics begin
        Λ' = λ
    end
    @derived begin
        event ~ @. TimeToEvent(λ, Λ)
    end
end

dt_sim = DataFrame(id = 1, time = 10., event=missing)
pop4tte = read_pumas(dt_sim, observations = [:event], event_data=false)
gamma_params = (; scale1 = 0.15, shape1 = 1.2,)
gamma_sims = simobstte(
    gamma_model, pop4tte[1], gamma_params, NamedTuple(); 
    minT = 0, nT = 100, repeated=false, maxT=500)
gamma_sims_df = DataFrame(gamma_sims)

Fit code:

num_subj = 10
dt_fit = DataFrame(id = 1:10, time = rand(Gamma(0.2, 10.), num_subj), event = float.(rand(Bernoulli(0.4), num_subj)))
pop4tte = read_pumas(dt_fit, observations = [:event], event_data=false)
gamma_params = (; scale1 = 0.15, shape1 = 10.2,)
gamma_fit = fit(gamma_model, pop4tte, gamma_params, Pumas.NaivePooled())

Fitting error:

julia> gamma_fit = fit(gamma_model, pop4tte, gamma_params, Pumas.NaivePooled())
[ Info: Checking the initial parameter values.
ERROR: MethodError: no method matching floatmax(::ForwardDiff.Dual{ForwardDiff.Tag{Pumas.Tag, Float64}, Float64, 2})
Closest candidates are:
  floatmax() at /opt/julia-1.7.3/share/julia/base/float.jl:903
  floatmax(::T) where T<:AbstractFloat at /opt/julia-1.7.3/share/julia/base/float.jl:900
  floatmax(::Type{T}) where T<:FixedPointNumbers.FixedPoint at /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/FixedPointNumbers/HAGk2/src/FixedPointNumbers.jl:119
  ...
Stacktrace:
  [1] En_cf_gamma(ν::ForwardDiff.Dual{ForwardDiff.Tag{Pumas.Tag, Float64}, Float64, 2}, z::ForwardDiff.Dual{ForwardDiff.Tag{Pumas.Tag, Float64}, Float64, 2}, n::Int64)
    @ SpecialFunctions /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/SpecialFunctions/hefUc/src/expint.jl:215
  [2] En_cf(ν::ForwardDiff.Dual{ForwardDiff.Tag{Pumas.Tag, Float64}, Float64, 2}, z::ForwardDiff.Dual{ForwardDiff.Tag{Pumas.Tag, Float64}, Float64, 2}, niter::Int64)
    @ SpecialFunctions /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/SpecialFunctions/hefUc/src/expint.jl:248
  [3] _expint(ν::ForwardDiff.Dual{ForwardDiff.Tag{Pumas.Tag, Float64}, Float64, 2}, z::ForwardDiff.Dual{ForwardDiff.Tag{Pumas.Tag, Float64}, Float64, 2}, niter::Int64, ::Val{false})
    @ SpecialFunctions /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/SpecialFunctions/hefUc/src/expint.jl:437
  [4] expint (repeats 2 times)
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/SpecialFunctions/hefUc/src/expint.jl:516 [inlined]
  [5] _gamma(a::ForwardDiff.Dual{ForwardDiff.Tag{Pumas.Tag, Float64}, Float64, 2}, x::ForwardDiff.Dual{ForwardDiff.Tag{Pumas.Tag, Float64}, Float64, 2})
    @ SpecialFunctions 
...
5

Right not, I think it might not be possible to fit such a model in Pumas. The problem is that the derivatives of the gamma pdf are incorrect

julia> using ForwardDiff

julia> using Distributions

julia> ForwardDiff.gradient(t -> pdf(Gamma(t[1], t[2]), t[3]), [3.1, 2.2, 0.0])
3-element Vector{Float64}:
 NaN
 NaN
 NaN

julia> ForwardDiff.gradient(t -> pdf(Gamma(t[1], t[2]), t[3]), [3.1, 2.2, 1e-100])
3-element Vector{Float64}:
 -9.163068198811177e-210
 -5.56513941603919e-212
  8.293852936161633e-112

and not yet supported at all for the cdf

julia> ForwardDiff.gradient(t -> cdf(Gamma(t[1], t[2]), t[3]), [3.1, 2.2, 1e-100])
ERROR: MethodError: no method matching _gamma_inc(::ForwardDiff.Dual{ForwardDiff.Tag{var"#19#20", Float64}, Float64, 3}, ::ForwardDiff.Dual{ForwardDiff.Tag{var"#19#20", Float64}, Float64, 3}, ::Int64)
Closest candidates are:
  _gamma_inc(::Float64, ::Float64, ::Integer) at ~/.julia/packages/SpecialFunctions/hefUc/src/gamma_inc.jl:860
  _gamma_inc(::BigFloat, ::BigFloat, ::Integer) at ~/.julia/packages/SpecialFunctions/hefUc/src/gamma_inc.jl:950
  _gamma_inc(::Float32, ::Float32, ::Integer) at ~/.julia/packages/SpecialFunctions/hefUc/src/gamma_inc.jl:956
  ...
Stacktrace:
  [1] gamma_inc(a::ForwardDiff.Dual{ForwardDiff.Tag{var"#19#20", Float64}, Float64, 3}, x::Float64, ind::Int64)
    @ SpecialFunctions ~/.julia/packages/SpecialFunctions/hefUc/src/gamma_inc.jl:858
  [2] gamma_inc(a::ForwardDiff.Dual{ForwardDiff.Tag{var"#19#20", Float64}, Float64, 3}, d::ForwardDiff.Dual{ForwardDiff.Tag{var"#19#20", Float64}, Float64, 3}, ind::Int64)
    @ ForwardDiff ~/.julia/packages/ForwardDiff/QdStj/src/dual.jl:773
  [3] gamma_inc(a::ForwardDiff.Dual{ForwardDiff.Tag{var"#19#20", Float64}, Float64, 3}, x::ForwardDiff.Dual{ForwardDiff.Tag{var"#19#20", Float64}, Float64, 3})
    @ SpecialFunctions ~/.julia/packages/SpecialFunctions/hefUc/src/gamma_inc.jl:858
  [4] gammacdf(k::ForwardDiff.Dual{ForwardDiff.Tag{var"#19#20", Float64}, Float64, 3}, θ::ForwardDiff.Dual{ForwardDiff.Tag{var"#19#20", Float64}, Float64, 3}, x::ForwardDiff.Dual{ForwardDiff.Tag{var"#19#20", Float64}, Float64, 3})
    @ StatsFuns ~/.julia/packages/StatsFuns/beHKc/src/distrs/gamma.jl:34
  [5] cdf
    @ ~/.julia/packages/Distributions/bQ6Gj/src/univariates.jl:636 [inlined]
  [6] #19
    @ ./REPL[58]:1 [inlined]
  [7] vector_mode_dual_eval!
    @ ~/.julia/packages/ForwardDiff/QdStj/src/apiutils.jl:37 [inlined]
  [8] vector_mode_gradient(f::var"#19#20", x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{var"#19#20", Float64}, Float64, 3, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#19#20", Float64}, Float64, 3}}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/QdStj/src/gradient.jl:106
  [9] gradient(f::Function, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{var"#19#20", Float64}, Float64, 3, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#19#20", Float64}, Float64, 3}}}, ::Val{true})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/QdStj/src/gradient.jl:0
 [10] gradient(f::Function, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{var"#19#20", Float64}, Float64, 3, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#19#20", Float64}, Float64, 3}}}) (repeats 2 times)
    @ ForwardDiff ~/.julia/packages/ForwardDiff/QdStj/src/gradient.jl:17
 [11] top-level scope
    @ REPL[58]:1

We are working adding support for computing these gradients but, meanwhile, you’ll probably have to try with other functional forms.

1 Like