I’ve looked a bit into this. I’ll go through some steps that might be generally useful whenever the model fitting fails. Hopefully it’s then easier for other people to benefit from them.
As mentioned above, it’s always a good idea to check the deviance value in the starting values. I think we’ll add a check for that to the fit method at some point and we’ll generally try to make the debugging process simpler. For now, the starting value should be checked manually with
julia> deviance(model3,dfP,param, Pumas.FOCEI())
┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/avuk1/src/integrator_interface.jl:319
┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/avuk1/src/integrator_interface.jl:319
┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/avuk1/src/integrator_interface.jl:319
┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/avuk1/src/integrator_interface.jl:319
┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/avuk1/src/integrator_interface.jl:319
┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/avuk1/src/integrator_interface.jl:319
NaN
Whenever, the starting value is either Inf of NaN it won’t be possible to fit the model. To figure out which of the subjects that are causing problems, we can use the (unexported) findinfluential function, i.e.
julia> Pumas.findinfluential(model3, dfP, param, Pumas.FOCEI())
┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/avuk1/src/integrator_interface.jl:319
┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/avuk1/src/integrator_interface.jl:319
┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/avuk1/src/integrator_interface.jl:319
┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/avuk1/src/integrator_interface.jl:319
┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/avuk1/src/integrator_interface.jl:319
┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/avuk1/src/integrator_interface.jl:319
5-element Array{Tuple{String,Float64},1}:
("12", NaN)
("14", Inf)
("15", Inf)
("9", 118.83800947482402)
("4", 91.37243210337728)
So it looks like there are some problems with subjects 12, 14, and 15. The problems can happen as part of the evaluation of the conditional (on the random effects) likelihood function or when computing the derivatives with respect to the random effects in order to compute the Laplace approximation of the marginal likelihood.
To investigate this, we can evaluate just the conditional likelihood for the three subjects at zero random effects, i.e.
julia> map(["12", "14", "15"]) do id
(id, Pumas.conditional_nll(model3, filter(i -> i.id == id, dfP)[1], param, (η=zeros(4),)))
end
3-element Array{Tuple{String,Float64},1}:
("12", 1.0186482312664653e35)
("14", 43.07465583007613)
("15", 43.07348006826408)
The values of the conditional (log) likelihood function look reasonable for subjects 14, and 15 but is extreme for subject 12.
In such a case, it’s useful to inspect data for possible problems. Indeed, in this case it looks like there might be a problem with the timings. For subject 12, there is a positive concentration at t = 0.0 but the first dose happens only at t=0.0625. This is the reason for the very large likelihood value and is most likely an error in the data.
For subjects 14, and 15 things are more involved and I think we might have an issue in Pumas. By default, we use a hybrid ODE solver that will automatically switch to a stiff solver when necessary. For subjects 14 and 15 it looks like a stiff solver is indeed needed. We can try to evaluate the conditional likelihood with a non-stiff solver
julia> Pumas.conditional_nll(model3, filter(i -> i.id=="14", dfP)[1], param, (η=zeros(4),), alg=Vern7())
┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/avuk1/src/integrator_interface.jl:319
Inf
julia> Pumas.conditional_nll(model3, filter(i -> i.id=="15", dfP)[1], param, (η=zeros(4),), alg=Vern7())
┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/avuk1/src/integrator_interface.jl:319
Inf
It looks like we are correctly swtiching to the stiff solver when computing just the conditional likelihood but when we compute the deviance, the swiching to the stiff solver seems to fail. We can instead try to evaluate the marginal likelihood (same as the deviance) but forcing a stiff solver
julia> deviance(model3, filter(i -> i.id=="14", dfP)[1], param, Pumas.FOCEI(), alg=Rodas5())
65.93266393946297
julia> deviance(model3, filter(i -> i.id=="15", dfP)[1], param, Pumas.FOCEI(), alg=Rodas5())
65.93031242137423
and indeed it worked just fine.