Simulating with vcov

Exactly what @vijay said! :slight_smile:

EDIT: I would say that Pumas has a clear advantage in ease of use compared to NONMEM.

Hello,

@mohamed82008 is there a way to send a manually derived vcov matrix, array/matrix type, to the simobs function? I see that the vcov(fit) return type is Pumas.CovarianceEstimate…

Yes, you can pass the second argument to simobs as a matrix instead (preferrably wrapped with Symmetric). Example:

simobs(model_fit, Symmetric(vcov(model_fit).Σ), samples = 10)

Note that vcov(model_fit).Σ is just the underlying covariance matrix estimate.

2 Likes

Hi, @mohamed82008

Do I understand this correctly? The above line of code simobs (model_fit, Symmetric(vcov(model_fit).Σ), samples = 10) mean that 10 subjects will be created from vcov(model_fit) with the variability and the SE. The created subjects will be simulated using the model_fit.

Almost. It’s 10 samples that will be generated. Each sample is an entire simulated population. The data structure you get at the end is a vector of vectors where if we have 3 subjects and 10 samples, the length of the vector of vectors will be 10 (number of samples) and each element of the vector of vectors will be a vector of length 3 (number of subjects) corresponding to the simulated population for a particular sample.

1 Like

@mohamed82008 I am attempting to use this feature again, and I’m running into problems.

julia>     test = simobs(i_fit, [pop], vcov(i_fit); samples=2)
ERROR: MethodError: no method matching simobs(::Pumas.FittedPumasModel{PumasModel{(tvtr_ratio = 1, tvCL = 1, tvVc = 1, tvDuration = 1, ωCL_BSV = 1, ωVc_BSV = 1, ωduration_BSV = 1, σ_add = 1, σ_prop = 1), 3, (:Central,), ParamSet{NamedTuple{(:tvtr_ratio, :tvCL, :tvVc, :tvDuration, :ωCL_BSV, :ωVc_BSV, :ωduration_BSV, :σ_add, :σ_prop), NTuple{9, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}, var"#5#14", Pumas.TimeDispatcher{var"#6#15", var"#7#16"}, Pumas.DCPChecker{Pumas.TimeDispatcher{var"#8#17", var"#9#18"}}, var"#11#20", Central1, var"#12#21", var"#13#22", Nothing}, Vector{Subject{NamedTuple{(:dv,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantCovar{NamedTuple{(:Nsubj_PBE, :rep_PBE, :occ, :period, :sequence, :seq_n, :formulation, :isT, :input_tr_ratio), Tuple{Int64, Int64, Int64, Int64, String, Int64, String, Int64, Float64}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, Optim.MultivariateOptimizationResults{Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}, Vector{Float64}, Float64, Float64, Vector{Optim.OptimizationState{Float64, Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}}}, Bool, NamedTuple{(:f_limit_reached, :g_limit_reached, :h_limit_reached, :time_limit, :callback, :f_increased), NTuple{6, Bool}}}, FOCE, Vector{Vector{Float64}}, NamedTuple{(:optimize_fn, :constantcoef, :omegas, :ensemblealg, :diffeq_options), Tuple{Pumas.DefaultOptimizeFN{Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, NamedTuple{(), Tuple{}}, Tuple{}, EnsembleSerial, NamedTuple{(:abstol, :reltol, :alg), Tuple{Float64, Float64, OrdinaryDiffEq.CompositeAlgorithm{Tuple{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5{1, true, LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.AutoSwitch{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5{1, true, LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}}}}}, ParamSet{NamedTuple{(:tvtr_ratio, :tvCL, :tvVc, :tvDuration, :ωCL_BSV, :ωVc_BSV, :ωduration_BSV, :σ_add, :σ_prop), NTuple{9, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}, Pumas.OptimState{NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, Optim.BFGSState{Vector{Float64}, Matrix{Float64}, Float64, Vector{Float64}}}}, ::Vector{Vector{Subject{NamedTuple{(:dv,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantCovar{NamedTuple{(:Nsubj_PBE, :rep_PBE, :occ, :period, :sequence, :seq_n, :formulation, :isT, :input_tr_ratio), Tuple{Int64, Int64, Int64, Int64, String, Int64, String, Int64, Float64}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}}, ::Symmetric{Float64, Matrix{Float64}}; samples=2)
Closest candidates are:
  simobs(::Pumas.FittedPumasModel, ::AbstractVector{T} where T<:Subject, ::AbstractMatrix) at /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/9AGx1/src/estimation/inference.jl:700 got unsupported keyword argument "samples"
  simobs(::Pumas.FittedPumasModel, ::AbstractVector{T} where T<:Subject, ::AbstractMatrix, ::Union{Nothing, AbstractVector{<:NamedTuple}}; samples, rng) at /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/9AGx1/src/estimation/inference.jl:700
  simobs(::Any, ::AbstractVector{T} where T<:Subject, ::Any, ::NamedTuple, ::Any...; kwargs...) at /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/9AGx1/src/models/model_api.jl:2129
  ...
Stacktrace:
 [1] top-level scope
   @ ~/data/code/complex_generics/modeling/Depo-SubQ Provera 104/Learn-Apply Paradigm Manuscript v4/FBE_MPA_SD_pl_BE_fxns copy.jl:253

Removing the pop for testing purposes, causes the following error:

julia>     test = simobs(i_fit, vcov(i_fit); samples=2)
ERROR: PosDefException: matrix is not positive definite; Cholesky factorization failed.
Stacktrace:
 [1] checkpositivedefinite
   @ /opt/julia-1.8.5/share/julia/stdlib/v1.8/LinearAlgebra/src/factorization.jl:18 [inlined]
 [2] cholesky!(A::Symmetric{Float64, Matrix{Float64}}, ::LinearAlgebra.NoPivot; check::Bool)
   @ LinearAlgebra /opt/julia-1.8.5/share/julia/stdlib/v1.8/LinearAlgebra/src/cholesky.jl:270
 [3] #cholesky#162
   @ /opt/julia-1.8.5/share/julia/stdlib/v1.8/LinearAlgebra/src/cholesky.jl:402 [inlined]
 [4] cholesky (repeats 2 times)
   @ /opt/julia-1.8.5/share/julia/stdlib/v1.8/LinearAlgebra/src/cholesky.jl:402 [inlined]
 [5] PDMat
   @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/PDMats/CbBv1/src/pdmat.jl:19 [inlined]
 [6] MvNormal(μ::Vector{Float64}, Σ::Symmetric{Float64, Matrix{Float64}})
   @ Distributions /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Distributions/uREy8/src/multivariate/mvnormal.jl:201
 [7] simobs(fpm::Pumas.FittedPumasModel{PumasModel{(tvtr_ratio = 1, tvCL = 1, tvVc = 1, tvDuration = 1, ωCL_BSV = 1, ωVc_BSV = 1, ωduration_BSV = 1, σ_add = 1, σ_prop = 1), 3, (:Central,), ParamSet{NamedTuple{(:tvtr_ratio, :tvCL, :tvVc, :tvDuration, :ωCL_BSV, :ωVc_BSV, :ωduration_BSV, :σ_add, :σ_prop), NTuple{9, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}, var"#5#14", Pumas.TimeDispatcher{var"#6#15", var"#7#16"}, Pumas.DCPChecker{Pumas.TimeDispatcher{var"#8#17", var"#9#18"}}, var"#11#20", Central1, var"#12#21", var"#13#22", Nothing}, Vector{Subject{NamedTuple{(:dv,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantCovar{NamedTuple{(:Nsubj_PBE, :rep_PBE, :occ, :period, :sequence, :seq_n, :formulation, :isT, :input_tr_ratio), Tuple{Int64, Int64, Int64, Int64, String, Int64, String, Int64, Float64}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, Optim.MultivariateOptimizationResults{Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}, Vector{Float64}, Float64, Float64, Vector{Optim.OptimizationState{Float64, Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}}}, Bool, NamedTuple{(:f_limit_reached, :g_limit_reached, :h_limit_reached, :time_limit, :callback, :f_increased), NTuple{6, Bool}}}, FOCE, Vector{Vector{Float64}}, NamedTuple{(:optimize_fn, :constantcoef, :omegas, :ensemblealg, :diffeq_options), Tuple{Pumas.DefaultOptimizeFN{Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, NamedTuple{(), Tuple{}}, Tuple{}, EnsembleSerial, NamedTuple{(:abstol, :reltol, :alg), Tuple{Float64, Float64, OrdinaryDiffEq.CompositeAlgorithm{Tuple{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5{1, true, LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.AutoSwitch{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5{1, true, LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}}}}}, ParamSet{NamedTuple{(:tvtr_ratio, :tvCL, :tvVc, :tvDuration, :ωCL_BSV, :ωVc_BSV, :ωduration_BSV, :σ_add, :σ_prop), NTuple{9, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}, Pumas.OptimState{NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, Optim.BFGSState{Vector{Float64}, Matrix{Float64}, Float64, Vector{Float64}}}}, population::Vector{Subject{NamedTuple{(:dv,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantCovar{NamedTuple{(:Nsubj_PBE, :rep_PBE, :occ, :period, :sequence, :seq_n, :formulation, :isT, :input_tr_ratio), Tuple{Int64, Int64, Int64, Int64, String, Int64, String, Int64, Float64}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, vcov::Symmetric{Float64, Matrix{Float64}}, randeffs::Nothing; samples::Int64, rng::TaskLocalRNG)
   @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/9AGx1/src/estimation/inference.jl:711
 [8] #simobs#1116
   @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/9AGx1/src/estimation/inference.jl:698 [inlined]
 [9] top-level scope
   @ ~/data/code/complex_generics/modeling/Depo-SubQ Provera 104/Learn-Apply Paradigm Manuscript v4/FBE_MPA_SD_pl_BE_fxns copy.jl:253

Here is the vcov matrix:

julia> vcov(i_fit)
9×9 Symmetric{Float64, Matrix{Float64}}:
  0.00866509      14.5721          17.5201         0.103529      0.00020401     0.000269455   0.000133468   3.8551e-6     -0.00075309
 14.5721       24529.7          35863.9          177.324         0.194283       1.09845       0.228999      0.000509347   -1.25094
 17.5201       35863.9              2.00567e8  -7063.99        -66.301        853.466         1.22517       8.06038      -84.9449
  0.103529       177.324        -7063.99           2.00763      -0.0166879      0.0630772     0.00222508   -0.00116127    -0.00329978
  0.00020401       0.194283       -66.301         -0.0166879     0.000942773   -0.00413658   -2.54599e-5    3.63152e-5    -0.000103562
  0.000269455      1.09845        853.466          0.0630772    -0.00413658     0.0198991     0.000127268  -0.000129573    0.000100419
  0.000133468      0.228999         1.22517        0.00222508   -2.54599e-5     0.000127268   2.31051e-5   -1.1039e-6     -8.48535e-6
  3.8551e-6        0.000509347      8.06038       -0.00116127    3.63152e-5    -0.000129573  -1.1039e-6     1.97774e-6    -8.51019e-6
 -0.00075309      -1.25094        -84.9449        -0.00329978   -0.000103562    0.000100419  -8.48535e-6   -8.51019e-6     0.000114244

It looks like I should be using pop and not [pop], as is in the documentation. Is that correct?

I want to simulate with different number of subjects than are in the population used for the estimation.

pop is the correct use, where did you see [pop]?

What’s the output of isposdef(vcov(i_fit))? I think it will be false which means that your matrix is not positive definite.

pop is the correct use, where did you see [pop] ?

What’s the output of isposdef(vcov(i_fit)) ? I think it will be false which means that your matrix is not positive definite.

It is false. Could you educate me a bit regarding the positive definite requirement for the covariance matrix? It’s been a long time since I’ve taken linear algebra.

  1. I simulated a simple 1-comp model without correlations on any of the random effects. Is this why the vcov matrix is not positive definite? I’m curious as to why simulating with the cov matrix from such a simple model would fail.

  2. If I were to make a model block and use the above vcov matrix as the random effects (i.e., the mvnormal matrix), and use simobs to generate a distribution of PK parameters (population-level)…essentially doing some of what I believe Pumas may be doing under the hood when you give simobs a vcov matrix, would this also error due to the vcov matrix not being positive definite?

The [population::Population,] in the docs means that it is optional. Maybe we need to clarify that. CC: @storopoli.

On positive definiteness, a covariance matrix is supposed to be positive definite in Pumas. Generally speaking, a covariance matrix is not positive definite if any of the following conditions is satisfied:

  • One of the parameters has a variance of 0.
  • Two parameters are perfectly correlated (correlation of 1) or anti-correlated (correlation of -1). This can happen when your model is not identifiable wrt your 2 parameters. A model can be non-identifiable because of the model structure (structural non-identifiability) or because of lack of data (practical non-identifiability).
  • There is multi-collinearity between your parameters. This is a generalization of the second condition but for more than 2 parameters.

The positive definiteness condition must be satisfied for any covariance matrix you use in Pumas, e.g. when defining MvNormal priors on population or subject-specific parameters or when sampling from such distributions.

In Pumas, vcov tries to estimate of the covariance matrix of the sampling distribution of your parameters to estimate the uncertainty around your maximum likelihood estimates. This covariance matrix estimation process assumes that the maximum likelihood parameters obtained in the fit must have been truly maximal, i.e. the optimization succeeded. An example of this not happening is if the optimization algorithm converged to an inflection point instead of a point that locally maximizes the (log) likelihood. Another example is if the optimization failed to improve the log likelihood because the line search algorithm used in BFGS failed. This can happen when there are numerical errors in the gradient or the objective function evaluation. There are many potential sources of numerical errors in NLME, e.g:

  • FO/FOCE/Laplace introduce some numerical error when estimating the marginal likelihood and its gradient especially if the inner optimization performed to get the empirical Bayes estimates was not perfect,
  • ODE solvers can introduce some numerical errors,
  • Dividing by a variable that is very close to 0 or taking the exponential of a variable that is very large in your model can also introduce numerical errors. This can happen when your parameters are poorly scaled, e.g. some parameters are very large in magnitude and others are very small.

When the fit is not successful, often the vcov estimate of the covariance matrix of your parameters’ sampling distribution will fail or will give you a matrix that is not positive definite. This can trigger an error when trying to sample from an MvNormal distribution using this matrix as a covariance matrix which is what simobs does under the hood when used in the way you did.

Without seeing your model and data, it will be difficult to diagnose why you are getting the error you are getting.

1 Like

This is documented in the Julia documentation.

1 Like

For some reason, vcov started giving positive definite matrices.

I now am getting an error with the following:

i_FBE_sims = simobs(fits_PBE[i], pop, vcovs_PBE[i]*uncertainty_factor; samples=Nreps_FBE, rng=rng)

I edited the population to have more subjects than were used in the fit.

What does this error indicate?

ERROR: BoundsError: attempt to access 0-element reinterpret(LabelledArrays.SLArray{Tuple{1}, Float64, 1, 1, (:Central,)}, ::Vector{StaticArraysCore.SVector{1, Float64}}) at index [1]
Stacktrace:
  [1] throw_boundserror(A::Base.ReinterpretArray{LabelledArrays.SLArray{Tuple{1}, Float64, 1, 1, (:Central,)}, 1, StaticArraysCore.SVector{1, Float64}, Vector{StaticArraysCore.SVector{1, Float64}}, false}, I::Tuple{Int64})
    @ Base ./abstractarray.jl:703
  [2] checkbounds
    @ ./abstractarray.jl:668 [inlined]
  [3] _getindex_ra
    @ ./reinterpretarray.jl:377 [inlined]
  [4] getindex
    @ ./reinterpretarray.jl:343 [inlined]
  [5] first(a::Base.ReinterpretArray{LabelledArrays.SLArray{Tuple{1}, Float64, 1, 1, (:Central,)}, 1, StaticArraysCore.SVector{1, Float64}, Vector{StaticArraysCore.SVector{1, Float64}}, false})
    @ Base ./abstractarray.jl:404
  [6] tovecs(v::Base.ReinterpretArray{LabelledArrays.SLArray{Tuple{1}, Float64, 1, 1, (:Central,)}, 1, StaticArraysCore.SVector{1, Float64}, Vector{StaticArraysCore.SVector{1, Float64}}, false})
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/9AGx1/src/models/model_api.jl:1715
  [7] tovecs
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/9AGx1/src/models/model_api.jl:1733 [inlined]
  [8] simobs(::PumasModel{(tvtr_ratio = 1, tvCL = 1, tvVc = 1, tvDuration = 1, ωCL_BSV = 1, ωVc_BSV = 1, σ_prop = 1), 2, (:Central,), ParamSet{NamedTuple{(:tvtr_ratio, :tvCL, :tvVc, :tvDuration, :ωCL_BSV, :ωVc_BSV, :σ_prop), NTuple{7, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}, var"#30#39", Pumas.TimeDispatcher{var"#31#40", var"#32#41"}, Pumas.DCPChecker{Pumas.TimeDispatcher{var"#33#42", var"#34#43"}}, var"#36#45", Central1, var"#37#46", var"#38#47", Nothing}, ::Subject{NamedTuple{(:dv,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantCovar{NamedTuple{(:rep_PBE, :Nreps_PBE, :rep_FBE, :Nreps_FBE, :input_tr_ratio, :Nsubj_PBE, :Nsubj_FBE, :sequence, :seq_n, :period, :occ, :formulation, :isT), Tuple{Int64, Int64, Int64, Int64, Float64, Int64, Int64, String, Int64, Int64, Int64, String, Int64}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}, ::NamedTuple{(:tvtr_ratio, :tvCL, :tvVc, :tvDuration, :ωCL_BSV, :ωVc_BSV, :σ_prop), NTuple{7, Float64}}, ::Nothing; obstimes::Nothing, ensemblealg::EnsembleSerial, diffeq_options::NamedTuple{(:alg, :callback, :abstol, :reltol), Tuple{OrdinaryDiffEq.CompositeAlgorithm{Tuple{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5{1, true, LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.AutoSwitch{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5{1, true, LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}, Nothing, Float64, Float64}}, rng::TaskLocalRNG, simulate_error::Val{true}, isfor_derived::Val{false}, return_model::Val{true})
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/9AGx1/src/models/model_api.jl:1539
  [9] _simobs!(::Vector{Vector{Pumas.SimulatedObservations{PumasModel{(tvtr_ratio = 1, tvCL = 1, tvVc = 1, tvDuration = 1, ωCL_BSV = 1, ωVc_BSV = 1, σ_prop = 1), 2, (:Central,), ParamSet{NamedTuple{(:tvtr_ratio, :tvCL, :tvVc, :tvDuration, :ωCL_BSV, :ωVc_BSV, :σ_prop), NTuple{7, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}, var"#30#39", Pumas.TimeDispatcher{var"#31#40", var"#32#41"}, Pumas.DCPChecker{Pumas.TimeDispatcher{var"#33#42", var"#34#43"}}, var"#36#45", Central1, var"#37#46", var"#38#47", Nothing}, Subject{NamedTuple{(:dv,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantCovar{NamedTuple{(:rep_PBE, :Nreps_PBE, :rep_FBE, :Nreps_FBE, :input_tr_ratio, :Nsubj_PBE, :Nsubj_FBE, :sequence, :seq_n, :period, :occ, :formulation, :isT), Tuple{Int64, Int64, Int64, Int64, Float64, Int64, Int64, String, Int64, Int64, Int64, String, Int64}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Symbol}}, Vector{Float64}}, Vector{Float64}, NamedTuple{(:tvtr_ratio, :tvCL, :tvVc, :tvDuration, :ωCL_BSV, :ωVc_BSV, :σ_prop), NTuple{7, Float64}}, NamedTuple{(:ηCL, :ηVc), Tuple{Float64, Float64}}, NamedTuple{(:rep_PBE, :Nreps_PBE, :rep_FBE, :Nreps_FBE, :input_tr_ratio, :Nsubj_PBE, :Nsubj_FBE, :sequence, :seq_n, :period, :occ, :formulation, :isT), Tuple{Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{String}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{String}, Vector{Int64}}}, NamedTuple{(:callback, :continuity, :saveat, :alg, :abstol, :reltol), Tuple{Nothing, Symbol, Vector{Float64}, OrdinaryDiffEq.CompositeAlgorithm{Tuple{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5{1, true, LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.AutoSwitch{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5{1, true, LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}, Float64, Float64}}, Pumas.PKPDAnalyticalSolution{StaticArraysCore.SVector{1, Float64}, 2, Vector{StaticArraysCore.SVector{1, Float64}}, Vector{Float64}, Vector{StaticArraysCore.SVector{1, Float64}}, Vector{StaticArraysCore.SVector{1, Float64}}, Pumas.Returns{NamedTuple{(:tr_ratio, :CL, :Vc), Tuple{Float64, Float64, Float64}}}, Pumas.AnalyticalPKPDProblem{StaticArraysCore.SVector{1, Float64}, Float64, false, Central1, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}, Pumas.Returns{NamedTuple{(:tr_ratio, :CL, :Vc), Tuple{Float64, Float64, Float64}}}}}, NamedTuple{(:duration, :bioav), Tuple{Vector{NamedTuple{(:Central,), Tuple{Float64}}}, Vector{NamedTuple{(:Central,), Tuple{Float64}}}}}, NamedTuple{(:tr_ratio, :CL, :Vc), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}}, NamedTuple{(:dv,), Tuple{Vector{Float64}}}}}}, ::PumasModel{(tvtr_ratio = 1, tvCL = 1, tvVc = 1, tvDuration = 1, ωCL_BSV = 1, ωVc_BSV = 1, σ_prop = 1), 2, (:Central,), ParamSet{NamedTuple{(:tvtr_ratio, :tvCL, :tvVc, :tvDuration, :ωCL_BSV, :ωVc_BSV, :σ_prop), NTuple{7, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}, var"#30#39", Pumas.TimeDispatcher{var"#31#40", var"#32#41"}, Pumas.DCPChecker{Pumas.TimeDispatcher{var"#33#42", var"#34#43"}}, var"#36#45", Central1, var"#37#46", var"#38#47", Nothing}, ::Vector{Subject{NamedTuple{(:dv,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantCovar{NamedTuple{(:rep_PBE, :Nreps_PBE, :rep_FBE, :Nreps_FBE, :input_tr_ratio, :Nsubj_PBE, :Nsubj_FBE, :sequence, :seq_n, :period, :occ, :formulation, :isT), Tuple{Int64, Int64, Int64, Int64, Float64, Int64, Int64, String, Int64, Int64, Int64, String, Int64}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, ::Vector{NamedTuple{(:tvtr_ratio, :tvCL, :tvVc, :tvDuration, :ωCL_BSV, :ωVc_BSV, :σ_prop), NTuple{7, Float64}}}, ::Nothing, ::TaskLocalRNG; diffeq_options::NamedTuple{(:alg, :callback, :abstol, :reltol), Tuple{OrdinaryDiffEq.CompositeAlgorithm{Tuple{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5{1, true, LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.AutoSwitch{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5{1, true, LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}, Nothing, Float64, Float64}}, isfor_derived::Val{false}, simulate_error::Val{true}, obstimes::Nothing, iterstart::Int64, iterstop::Int64, return_model::Val{true})
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/9AGx1/src/models/model_api.jl:1819
 [10] _simobs(::PumasModel{(tvtr_ratio = 1, tvCL = 1, tvVc = 1, tvDuration = 1, ωCL_BSV = 1, ωVc_BSV = 1, σ_prop = 1), 2, (:Central,), ParamSet{NamedTuple{(:tvtr_ratio, :tvCL, :tvVc, :tvDuration, :ωCL_BSV, :ωVc_BSV, :σ_prop), NTuple{7, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}, var"#30#39", Pumas.TimeDispatcher{var"#31#40", var"#32#41"}, Pumas.DCPChecker{Pumas.TimeDispatcher{var"#33#42", var"#34#43"}}, var"#36#45", Central1, var"#37#46", var"#38#47", Nothing}, ::Vector{Subject{NamedTuple{(:dv,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantCovar{NamedTuple{(:rep_PBE, :Nreps_PBE, :rep_FBE, :Nreps_FBE, :input_tr_ratio, :Nsubj_PBE, :Nsubj_FBE, :sequence, :seq_n, :period, :occ, :formulation, :isT), Tuple{Int64, Int64, Int64, Int64, Float64, Int64, Int64, String, Int64, Int64, Int64, String, Int64}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, ::Vector{NamedTuple{(:tvtr_ratio, :tvCL, :tvVc, :tvDuration, :ωCL_BSV, :ωVc_BSV, :σ_prop), NTuple{7, Float64}}}, ::Nothing; rng::TaskLocalRNG, diffeq_options::NamedTuple{(:abstol, :reltol, :alg), Tuple{Float64, Float64, OrdinaryDiffEq.CompositeAlgorithm{Tuple{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5{1, true, LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.AutoSwitch{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5{1, true, LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}}}, isfor_derived::Val{false}, simulate_error::Val{true}, obstimes::Nothing, ensemblealg::EnsembleSerial)
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/9AGx1/src/models/model_api.jl:2102
 [11] simobs(::PumasModel{(tvtr_ratio = 1, tvCL = 1, tvVc = 1, tvDuration = 1, ωCL_BSV = 1, ωVc_BSV = 1, σ_prop = 1), 2, (:Central,), ParamSet{NamedTuple{(:tvtr_ratio, :tvCL, :tvVc, :tvDuration, :ωCL_BSV, :ωVc_BSV, :σ_prop), NTuple{7, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}, var"#30#39", Pumas.TimeDispatcher{var"#31#40", var"#32#41"}, Pumas.DCPChecker{Pumas.TimeDispatcher{var"#33#42", var"#34#43"}}, var"#36#45", Central1, var"#37#46", var"#38#47", Nothing}, ::Vector{Subject{NamedTuple{(:dv,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantCovar{NamedTuple{(:rep_PBE, :Nreps_PBE, :rep_FBE, :Nreps_FBE, :input_tr_ratio, :Nsubj_PBE, :Nsubj_FBE, :sequence, :seq_n, :period, :occ, :formulation, :isT), Tuple{Int64, Int64, Int64, Int64, Float64, Int64, Int64, String, Int64, Int64, Int64, String, Int64}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, ::Vector{NamedTuple{(:tvtr_ratio, :tvCL, :tvVc, :tvDuration, :ωCL_BSV, :ωVc_BSV, :σ_prop), NTuple{7, Float64}}}, ::Nothing; rng::TaskLocalRNG, kwargs::Base.Pairs{Symbol, NamedTuple{(:abstol, :reltol, :alg), Tuple{Float64, Float64, OrdinaryDiffEq.CompositeAlgorithm{Tuple{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5{1, true, LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.AutoSwitch{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5{1, true, LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}}}, Tuple{Symbol}, NamedTuple{(:diffeq_options,), Tuple{NamedTuple{(:abstol, :reltol, :alg), Tuple{Float64, Float64, OrdinaryDiffEq.CompositeAlgorithm{Tuple{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5{1, true, LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.AutoSwitch{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5{1, true, LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}}}}}})
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/9AGx1/src/models/model_api.jl:2008
 [12] simobs(fpm::Pumas.FittedPumasModel{PumasModel{(tvtr_ratio = 1, tvCL = 1, tvVc = 1, tvDuration = 1, ωCL_BSV = 1, ωVc_BSV = 1, σ_prop = 1), 2, (:Central,), ParamSet{NamedTuple{(:tvtr_ratio, :tvCL, :tvVc, :tvDuration, :ωCL_BSV, :ωVc_BSV, :σ_prop), NTuple{7, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}, var"#30#39", Pumas.TimeDispatcher{var"#31#40", var"#32#41"}, Pumas.DCPChecker{Pumas.TimeDispatcher{var"#33#42", var"#34#43"}}, var"#36#45", Central1, var"#37#46", var"#38#47", Nothing}, Vector{Subject{NamedTuple{(:dv,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantCovar{NamedTuple{(:rep_PBE, :Nreps_PBE, :rep_FBE, :Nreps_FBE, :input_tr_ratio, :Nsubj_PBE, :Nsubj_FBE, :sequence, :seq_n, :period, :occ, :formulation, :isT), Tuple{Int64, Int64, Int64, Int64, Float64, Int64, Int64, String, Int64, Int64, Int64, String, Int64}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, Optim.MultivariateOptimizationResults{Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}, Vector{Float64}, Float64, Float64, Vector{Optim.OptimizationState{Float64, Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}}}, Bool, NamedTuple{(:f_limit_reached, :g_limit_reached, :h_limit_reached, :time_limit, :callback, :f_increased), NTuple{6, Bool}}}, FOCE, Vector{Vector{Float64}}, NamedTuple{(:optimize_fn, :constantcoef, :omegas, :ensemblealg, :diffeq_options), Tuple{Pumas.DefaultOptimizeFN{Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, NamedTuple{(), Tuple{}}, Tuple{}, EnsembleSerial, NamedTuple{(:abstol, :reltol, :alg), Tuple{Float64, Float64, OrdinaryDiffEq.CompositeAlgorithm{Tuple{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5{1, true, LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.AutoSwitch{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5{1, true, LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}}}}}, ParamSet{NamedTuple{(:tvtr_ratio, :tvCL, :tvVc, :tvDuration, :ωCL_BSV, :ωVc_BSV, :σ_prop), NTuple{7, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}, Pumas.OptimState{NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, Optim.BFGSState{Vector{Float64}, Matrix{Float64}, Float64, Vector{Float64}}}}, population::Vector{Subject{NamedTuple{(:dv,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantCovar{NamedTuple{(:rep_PBE, :Nreps_PBE, :rep_FBE, :Nreps_FBE, :input_tr_ratio, :Nsubj_PBE, :Nsubj_FBE, :sequence, :seq_n, :period, :occ, :formulation, :isT), Tuple{Int64, Int64, Int64, Int64, Float64, Int64, Int64, String, Int64, Int64, Int64, String, Int64}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, vcov::Symmetric{Float64, Matrix{Float64}}, randeffs::Nothing; samples::Int64, rng::TaskLocalRNG)
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/9AGx1/src/estimation/inference.jl:713
 [13] top-level scope
    @ ~/data/code/complex_generics/modeling/Depo-SubQ Provera 104/Learn-Apply Paradigm Manuscript v4/FBE_MPA_SD_pl_BE_fxns copy.jl:266

I have serendipitously figured out the problem. This is the error you get when you simulate with a Pumas population that only has dosing times and no obstimes.

Is adding in obstimes possible with simobs(fit, pop, vcov(fit); samples=x)?

I tried the following which do not work:

simobs(fit, pop, vcov(fit), obstimes=obstimes; samples=x)
simobs(fit, pop, vcov(fit); samples=x, obstimes=obstimes)

I was able to manually add obstimes to the Pumas population I created with some data manipulation, but being able to add in obstimes from simobs would be much more convenient for users when simulating with parameter uncertainty.

How are you creating your population? If you are using the Subject in a map call to create it, then you can assign the times to that directly. Look for ?Subject to see the help

1 Like

Thank you, @vijay.

I wasn’t aware you could add times to Subject.