Hello Pumas community!
I am attempting to simulate an ER model using logistic regression. The simulation runs when simulating a single dose regimen, but gives the following error when adding an additional dose. The error is posted at the end due to the length of the stacktrace.
Parameters:
params_sim = (
# PK
tvcl = 18,
tvvc = 93,
# Logistic Regression - Endpoint
Interc = -0.82,
BetaAUC = 0.4,
# RUV
σ_prop = 0.01,
σ_add = 0.01,
)
Model:
mdl_auc_sim = @model begin
@param begin
# PK
tvcl ∈ RealDomain()
tvvc ∈ RealDomain()
# Logistic Regression - Endpoint
Interc ∈ RealDomain()
BetaAUC ∈ RealDomain()
# RUV
σ_prop ∈ RealDomain()
σ_add ∈ RealDomain()
end
@covariates regimen wt
@pre begin
CL = tvcl * (wt/70)^0.75
Vc = tvvc * (wt/70)^1
_BetaAUC = BetaAUC
_Interc = Interc
end
@vars begin
cp := Central / Vc
end
@dynamics begin
Central' = -CL/Vc * Central
end
@derived begin
dv ~ @. Normal(cp, sqrt((σ_prop * cp)^2 + σ_add^2))
nca_cp := @. @nca cp # compute NCA on cp
AUC = NCA.auc(nca_cp, method=:linuplogdown)
cmax = NCA.cmax(nca_cp)
linear_pred = AUC * _BetaAUC
PAINRELIEF ~ @. Bernoulli(logistic(_Interc + linear_pred))
end
end
Problematic Dosing Regimen and Population for Simulation:
dose_10bidbolus = DosageRegimen(10; time = 0, ii=12, addl=1, route=NCA.IVBolus)
pop_10bid = map(101:200) do i
Subject(id = i,
events = dose_10bidbolus,
covariates = (wt = rand(55:80), regimen = "10mg BID IV Bolus"),
observations = (; dv = nothing,
PAINRELIEF = nothing)
)
end
pop_10bid_obs = DataFrame(simobs(mdl_auc_sim, pop_10bid, params_sim, obstimes = 0:2:24))
The following error occurs on the last line with the simobs step:
ERROR: TaskFailedException
nested task error: MethodError: no method matching *(::Vector{Float64}, ::Vector{Float64})
The function `*` exists, but no method is defined for this combination of argument types.
Closest candidates are:
*(::Any, ::Any, ::Any, ::Any...)
@ Base operators.jl:596
*(::MutableArithmetics.Zero, ::Any)
@ MutableArithmetics C:\Users\rayga\.julia\packages\MutableArithmetics\tNSBd\src\rewrite.jl:61
*(::Any, ::ChainRulesCore.NoTangent)
@ ChainRulesCore C:\Users\rayga\.julia\packages\ChainRulesCore\XAgYn\src\tangent_arithmetic.jl:65
...
Stacktrace:
[1] (::var"#15#21")(_pre::Returns{…}, _sol::Pumas.PKPDAnalyticalSolution{…}, _obstimes::StepRange{…}, _subject::Subject{…}, _param::@NamedTuple{…}, _random::@NamedTuple{})
@ Main .\none:2214
[2] DerivedObj
@ .\none:2245 [inlined]
[3] _derived
@ .\none:3173 [inlined]
[4] _simobs(model::PumasModel{…}, subject::Subject{…}, param::@NamedTuple{…}, randeffs::Nothing; obstimes::StepRange{…}, ensemblealg::EnsembleSerial, diffeq_options::@NamedTuple{…}, rng::Xoshiro, simulate_error::Val{…}, isfor::Val{…})
@ Pumas .\none:2263
[5] _simobs
@ .\none:2198 [inlined]
[6] _simobs_worker!(sims::Vector{…}, model::PumasModel{…}, population::Vector{…}, vparam::StaticArraysCore.SVector{…}, vrandeffs::Nothing, rng::Xoshiro, seeds::Vector{…}, ::EnsembleSerial; diffeq_options::@NamedTuple{…}, isfor::Val{…}, simulate_error::Val{…}, obstimes::StepRange{…}, return_model::Val{…}, iterstart::Int64, iterstop::Int64, outstart::Int64)
@ Pumas .\none:2918
[7] _simobs_worker!
@ .\none:2881 [inlined]
[8] macro expansion
@ .\none:2986 [inlined]
[9] (::Pumas.var"#610#threadsfor_fun#436"{Pumas.var"#610#threadsfor_fun#435#437"{…}})(tid::Int64; onethread::Bool)
@ Pumas .\threadingconstructs.jl:253
[10] #610#threadsfor_fun
@ .\threadingconstructs.jl:220 [inlined]
[11] (::Base.Threads.var"#1#2"{Pumas.var"#610#threadsfor_fun#436"{Pumas.var"#610#threadsfor_fun#435#437"{…}}, Int64})()
@ Base.Threads .\threadingconstructs.jl:154
...and 1 more exception.
Stacktrace:
[1] threading_run(fun::Pumas.var"#610#threadsfor_fun#436"{Pumas.var"#610#threadsfor_fun#435#437"{…}}, static::Bool)
@ Base.Threads .\threadingconstructs.jl:173
[2] macro expansion
@ .\threadingconstructs.jl:190 [inlined]
[3] #_simobs_worker!#434
@ .\none:2983 [inlined]
[4] _simobs_worker!
@ .\none:2957 [inlined]
[5] #_simobs!#429
@ .\none:2740 [inlined]
[6] _simobs!
@ .\none:2725 [inlined]
[7] _simobs(model::PumasModel{…}, population::Vector{…}, vparam::StaticArraysCore.SVector{…}, vrandeffs::Nothing; rng::TaskLocalRNG, diffeq_options::@NamedTuple{…}, isfor::Val{…}, simulate_error::Val{…}, obstimes::StepRange{…}, ensemblealg::EnsembleThreads)
@ Pumas .\none:2705
[8] _simobs
@ .\none:2639 [inlined]
[9] #simobs#422
@ .\none:2545 [inlined]
[10] simobs
@ .\none:2530 [inlined]
[11] top-level scope
@ c:\Users\rayga\OneDrive\Desktop\Poster with PumasAI Post-Fellowship\mwe.jl:66
Some type information was truncated. Use `show(err)` to see complete types.
This is using Pumas version 2.7.0.
I am not sure what is causing this issue as the single-dose regimen simulation executed and a similar layout was used to create the populations. I tested different methods of creating the population, using either a , or a ; in the dosing regimen for better separation, using .* vs * for the linear_pred = AUC * _BetaAUC step in @derived, but these have not resolved the issue. Other ideas are welcome; please let me know if any additional information is needed.
If needed, the working single-dose regimen is included:
dose_10qdbolus = DosageRegimen(10, time = 0, route=NCA.IVBolus)
pop_10qd = map(1:100) do i
Subject(id = i,
events = dose_10qdbolus,
covariates = (wt = rand(55:80), regimen = "10mg QD IV Bolus"),
observations = (; dv = nothing,
PAINRELIEF = nothing)
)
end
pop_10qd_obs = DataFrame(simobs(mdl_auc_sim, pop_10qd, params_sim, obstimes = 0:2:24))