Cholesky factorization failed - Simulation for Population Not Working

1

Hi I am trying to run Exercise PK-27 from the tutorials. For a population I am not able to run the simulation, it shows the following error. But the same code works fine for a single subject. Please help …!

sim_pop3_sub = simobs(pk_19, pop3_sub, param, obstimes=0.1:1:300)
ERROR: PosDefException: matrix is not positive definite; Cholesky factorization failed.
Stacktrace:
  [1] checkpositivedefinite
    @ /opt/julia-1.7.2/share/julia/stdlib/v1.7/LinearAlgebra/src/factorization.jl:18 [inlined]
  [2] cholesky!(A::Hermitian{Float64, Matrix{Float64}}, ::Val{false}; check::Bool)
    @ LinearAlgebra /opt/julia-1.7.2/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:266
  [3] cholesky!(A::Matrix{Float64}, ::Val{false}; check::Bool)
    @ LinearAlgebra /opt/julia-1.7.2/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:298
  [4] #cholesky#143
    @ /opt/julia-1.7.2/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:394 [inlined]
  [5] cholesky (repeats 2 times)
    @ /opt/julia-1.7.2/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:394 [inlined]
  [6] PDMat
    @ /build/_work/PumasSystemImages/PumasSystemImages/.julia/packages/PDMats/ZW0lN/src/pdmat.jl:19 [inlined]
  [7] MvNormal(μ::FillArrays.Zeros{Float64, 1, Tuple{Base.OneTo{Int64}}}, Σ::Matrix{Float64})
    @ Distributions /build/_work/PumasSystemImages/PumasSystemImages/.julia/packages/Distributions/7iOJp/src/multivariate/mvnormal.jl:201
  [8] MvNormal
    @ /build/_work/PumasSystemImages/PumasSystemImages/.julia/packages/Distributions/7iOJp/src/multivariate/mvnormal.jl:218 [inlined]
  [9] #5
    @ /build/_work/PumasSystemImages/PumasSystemImages/.julia/packages/Pumas/89e4W/src/dsl/model_macro.jl:334 [inlined]
 [10] sample_randeffs
    @ /build/_work/PumasSystemImages/PumasSystemImages/.julia/packages/Pumas/89e4W/src/models/model_api.jl:627 [inlined]
 [11] #371
    @ ./none:0 [inlined]
 [12] iterate
    @ ./generator.jl:47 [inlined]
 [13] collect(itr::Base.Generator{Vector{NamedTuple{(:tvvc, :tvvp, :tvq, :tvvmax, :tvkm, :tvkme, :tvvme, :Ω, :σ_prop_cp, :σ_prop_met), Tuple{Float64, Float64, Float64, Float64, Float64, Float64, Float64, Matrix{Float64}, Float64, Float64}}}, Pumas.var"#371#373"{TaskLocalRNG, PumasModel{(tvvc = 1, tvvp = 1, tvq = 1, tvvmax = 1, tvkm = 1, tvkme = 1, tvvme = 1, Ω = 7, σ_prop_cp = 1, σ_prop_met = 1), 7, (:Central, :Peripheral, :Metabolite), ParamSet{NamedTuple{(:tvvc, :tvvp, :tvq, :tvvmax, :tvkm, :tvkme, :tvvme, :Ω, :σ_prop_cp, :σ_prop_met), Tuple{RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}}}}, var"#5#12", var"#6#13", Nothing, var"#7#14", ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, ODEFunction{false, ModelingToolkit.ODEFunctionClosure{var"#8#15", var"#9#16"}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, var"#10#17", var"#11#18", ModelingToolkit.ODESystem}}})
    @ Base ./array.jl:724
 [14] _simobs(::PumasModel{(tvvc = 1, tvvp = 1, tvq = 1, tvvmax = 1, tvkm = 1, tvkme = 1, tvvme = 1, Ω = 7, σ_prop_cp = 1, σ_prop_met = 1), 7, (:Central, :Peripheral, :Metabolite), ParamSet{NamedTuple{(:tvvc, :tvvp, :tvq, :tvvmax, :tvkm, :tvkme, :tvvme, :Ω, :σ_prop_cp, :σ_prop_met), Tuple{RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}}}}, var"#5#12", var"#6#13", Nothing, var"#7#14", ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, ODEFunction{false, ModelingToolkit.ODEFunctionClosure{var"#8#15", var"#9#16"}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, var"#10#17", var"#11#18", ModelingToolkit.ODESystem}, ::Vector{Subject{NamedTuple{(:cp, :met), Tuple{Vector{Missing}, Vector{Missing}}}, Pumas.ConstantCovar{NamedTuple{(), Tuple{}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Nothing}}, ::Pumas.RepeatedVector{NamedTuple{(:tvvc, :tvvp, :tvq, :tvvmax, :tvkm, :tvkme, :tvvme, :Ω, :σ_prop_cp, :σ_prop_met), Tuple{Float64, Float64, Float64, Float64, Float64, Float64, Float64, Diagonal{Float64, Vector{Float64}}, Float64, Float64}}, NamedTuple{(:tvvc, :tvvp, :tvq, :tvvmax, :tvkm, :tvkme, :tvvme, :Ω, :σ_prop_cp, :σ_prop_met), Tuple{Float64, Float64, Float64, Float64, Float64, Float64, Float64, Diagonal{Float64, Vector{Float64}}, Float64, Float64}}}, ::Nothing; obstimes::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, diffeq_options::NamedTuple{(), Tuple{}}, ensemblealg::EnsembleSerial, isfor_derived::Val{false}, rng::TaskLocalRNG, simulate_error::Val{true})
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/.julia/packages/Pumas/89e4W/src/models/model_api.jl:2054
 [15] simobs(::PumasModel{(tvvc = 1, tvvp = 1, tvq = 1, tvvmax = 1, tvkm = 1, tvkme = 1, tvvme = 1, Ω = 7, σ_prop_cp = 1, σ_prop_met = 1), 7, (:Central, :Peripheral, :Metabolite), ParamSet{NamedTuple{(:tvvc, :tvvp, :tvq, :tvvmax, :tvkm, :tvkme, :tvvme, :Ω, :σ_prop_cp, :σ_prop_met), Tuple{RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}}}}, var"#5#12", var"#6#13", Nothing, var"#7#14", ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, ODEFunction{false, ModelingToolkit.ODEFunctionClosure{var"#8#15", var"#9#16"}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, var"#10#17", var"#11#18", ModelingToolkit.ODESystem}, ::Vector{Subject{NamedTuple{(:cp, :met), Tuple{Vector{Missing}, Vector{Missing}}}, Pumas.ConstantCovar{NamedTuple{(), Tuple{}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Nothing}}, ::NamedTuple{(:tvvc, :tvvp, :tvq, :tvvmax, :tvkm, :tvkme, :tvvme, :Ω, :σ_prop_cp, :σ_prop_met), Tuple{Float64, Float64, Float64, Float64, Float64, Float64, Float64, Diagonal{Float64, Vector{Float64}}, Float64, Float64}}, ::Nothing; rng::TaskLocalRNG, kwargs::Base.Pairs{Symbol, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Tuple{Symbol}, NamedTuple{(:obstimes,), Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}})
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/.julia/packages/Pumas/89e4W/src/models/model_api.jl:1677
 [16] top-level scope
    @ ~/data/code/PBPK_Voriconazole.jl:98
2

I can tell you that it is because the variance elements are all zero, but I have to look into why it used to work and why it doesn’t not work now.

A work around is to either specify positive diagonal elements or to input the zero random effects if that is what you want. To do the former, do something like this

param = ( tvcl        = 0.001,
          tvkon       = 0.096,
          tvkoff      = 0.001,
          tvvp        = 0.100,
          tvq         = 0.003,
          tvkin       = 0.11,
          tvkout      = 0.0089,
          tvkerl      = 0.003,
          tvvc        = 0.05,
          Ω           = Diagonal([0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15]),
          σ²_prop_cp  = 0.02,
          σ²_prop_rec = 0.012,
          σ²_prop_com = 0.015)

To do the latter, do something like this instead (then you don’t have to update the omega as above):

sim_pop4 = simobs(pk_27, pop4_sub, param, init_randeffs(pk_27, param, pop4_sub); obstimes = 0.1:1:500)

This will set all the random effects to 0 for gaussian random effects.

3

That thing works. Thanks @pkofod