Weibull Model Fit Gives NaN Gradient

Estimation
1

Hi - I am trying to fit a weibull model on some data, but it fails. If I remove the exponent on t and fit an exponential function it does not fail. I have looked at the simulated data from both models and it looks the way it should, so model mis-specification is not a problem. Is this NaN gradient a bug or am I doing something wrong?
Code is attached below. I have tested this code on Pumas Desktop 2.0 and also Pumas 2.2.1 on the cloud and the result is same in both.

Code with exponent - doesn't work : 
using Pumas, PumasUtilities

weibull_model = @model begin
    @param begin
        PEmax ∈ RealDomain()
        TD ∈ RealDomain(lower = 0)
        POW ∈ RealDomain()
        Ω   ∈ PDiagDomain(3)
        σ_add ∈ RealDomain(lower=0.001)
    end
    @random η ~ MvNormal(Ω)
    @covariates BASE_SCORE
    @pre begin
        PEmax_i = PEmax + η[1]
        TD_i    = TD * exp(η[2])
        POW_i   = POW + η[3]
        bsl_i   = BASE_SCORE
    end
    @vars begin
        score_i_t := bsl_i * (1 - PEmax_i * (1 - exp( - (t/TD_i)^POW_i )))
    end
    @derived begin
        score_out ~ @. Normal(score_i_t, σ_add)
    end
end

weibull_params = (PEmax = 0.45, TD = 4, POW = 0.1, Ω = Diagonal(fill(0.09, 3)), σ_add = 4)
weibull_sims = simobs(
    weibull_model, 
    [Subject(id = i, covariates = (BASE_SCORE = 50, )) for i in 1:20],
    weibull_params,
    obstimes = [0,1,2,4,6,8,10]
)
weibull_refit_pop = Subject.(weibull_sims)
weibull_refit = fit(weibull_model, weibull_refit_pop, weibull_params, Pumas.FOCE())
Code without exponent - works : 
using Pumas, PumasUtilities

weibull_model = @model begin
    @param begin
        PEmax ∈ RealDomain()
        TD ∈ RealDomain(lower = 0)
        Ω   ∈ PDiagDomain(2)
        σ_add ∈ RealDomain(lower=0.001)
    end
    @random η ~ MvNormal(Ω)
    @covariates BASE_SCORE
    @pre begin
        PEmax_i = PEmax + η[1]
        TD_i    = TD * exp(η[2])
        bsl_i   = BASE_SCORE
    end
    @vars begin
        score_i_t := bsl_i * (1 - PEmax_i * (1 - exp( - (t/TD_i) )))
    end
    @derived begin
        score_out ~ @. Normal(score_i_t, σ_add)
    end
end

weibull_params = (PEmax = 0.45, TD = 4, Ω = Diagonal(fill(0.09, 2)), σ_add = 4)
weibull_sims = simobs(
    weibull_model, 
    [Subject(id = i, covariates = (BASE_SCORE = 50, )) for i in 1:20],
    weibull_params,
    obstimes = [0,1,2,4,6,8,10]
)
weibull_refit_pop = Subject.(weibull_sims)
weibull_refit = fit(weibull_model, weibull_refit_pop, weibull_params, Pumas.FOCE())
2

Hi,

The way that you’ve specified your model, your POW_i can go negative. You’re also evaluating the model at t=0 which means that you’ll occasionally be evaluating 0^POW_i and that’ll return Inf - resulting in an infinitely bad loglikelihood. You might either ensure that POW_i stays positive or that you don’t have 0 in your obstimes.

Also, I have not myself used the @vars block without having @dynamics the way that you do. I would just have placed

score_i_t = bsl_i * (1 - PEmax_i * (1 - exp( - (t/TD_i)^POW_i )))

right in the @pre block. I’m not at all sure that your solution is problematic, I just reacted to it and thought I’d mention this in case it might help.

1 Like
3

I had the same problem as you and in fact I wanted also to account for both negative and positive trends ( to account for the worsening and improvement in the outcome). To run the estimation task and to circumvent the t=0 issue I started the equation at a very low value as following:

@vars begin
      Placebo := @. Pmax*(1-exp(-((t/TD+1e-5)^POW)))
      Drug := ((Emax*average)/(EC50+average))
      Delay := @. (1 - exp(-(Kt*t)))
  end
  @derived begin
    PANSS := @. baseline * (1 - Placebo) * (1 - Drug*Delay)
    base := baseline
    Observed ~ @. Normal(PANSS,σ_add)
  end
end

Hope this can be helpful for your case as well.

1 Like
4

Thank you, Niklas and Ahmed for your insights and solutions. I tested your approach where I constrained the POW to be greater than 0 and it worked, because I allowed the PEmax to be less than 0, I was able to capture both improving and worsening trends.
Also, I think the simulations were not working the way I expected when I included the equation in the @pre block but worked rather when I put it in the @vars block - which is the primary reason to put the equation in the @vars block