# Corssvalidation: BySubject(marginal = LaplaceI())

Hello,

I would like to know if there is a difference in crossvalidation results or situations where I should set BySubject(marginal = LaplaceI()) versusBySubject(marginal = FOCE()) versus BySubject(marginal = nothing

Thanks

This part is explained in Bayesian Workflow Β· Pumas. The difference is how we marginalise the random effects or if we do at all. Setting marginal to nothing means the random effects will be set to 0 (assuming itβs Gaussian) and we will be computing the conditional likelihood of the fixed effects given the left out data and 0 random effects. Setting marginal to LaplaceI() or FOCE() will marginalise the random effects so we will be computing the marginal likelihood of the fixed effects given the left out data. This is generally preferred when comparing models because the marginal likelihood evaluates the entire distribution on the random effects not just the βmeanβ of 0 so it can identify a bad \Omega for example. But both can be useful to inspect because setting the random effects to 0 when making predictions for a new subject is common practice so that likelihood is important.

1 Like

Thanks @mohamed82008 It is much clearer for me now. I tried to run the cross validation task using the following code:

BAYESIAN_one_cmt_add_ecmo = @model begin
@param begin
#tvKTV1 ~ LogNormal(log(0.4), 1.0)
#tvKTV2 ~ LogNormal(log(0.1), 1.0)
#tvecmo ~ Normal(-0.7, 1.0)
#tvCL1 ~ LogNormal(log(0.9), 1.0)
tvCL2 ~ LogNormal(log(0.5), 1.0)
#tvcircuit ~ LogNormal(log(1.2), 1.0)
#tvVc1 ~ LogNormal(log(1.8), 0.1)
tvVc2 ~ LogNormal(log(0.5), 1.0)
#Ο_prop_1 ~ truncated(Cauchy(0.18, 5), 0.0, Inf)
Ο_add_2 ~ truncated(Cauchy(10000, 5), 0.0, Inf)
Ο_add_3 ~ truncated(Cauchy(10000, 5), 0.0, Inf)
#Ο_add ~ truncated(Cauchy(3500, 5), 0.0, Inf)
#ΟΒ²KTV ~ truncated(Normal(1.5, 10), 0.0, Inf)
Ο ~ MvLogNormal(ones(2), Diagonal([1.0,1.0]))

end

@random begin
Ξ· ~ MvNormal(Ο)
end
@covariates WT ECMO_DAYS IS_CIRCUIT_CHANGE site T
@pre begin

# PK parameters

siteeffect = if site==2
else
end
#siteeffect_2 = site==2 ? Ο_add : 0.0

#tvKTV = site == 1 ? tvKTV1 : tvKTV2*(ECMO_DAYS/6)^tvecmo
#KTV    = tvKTV*exp(Ξ·KTV)
#TVCL = site == 1 ? tvCL1 * (WT/10)^0.75 : tvCL2* (WT/10)^0.75  * tvcircuit^IS_CIRCUIT_CHANGE
CL     = tvCL2 * exp(Ξ·[1])#*(1-exp(-KTV*t))
#TVVc = site == 1 ? tvVc1 * (WT/10)^1 : tvVc2 * (WT/10)^1
Vc     = tvVc2 * exp(Ξ·[2])
end

@dynamics Central1
@derived begin
CP = @. Central / Vc
CONC = @. Normal(CP, sqrt(siteeffect))
end
end

one_cmt_add = truncate(one_cmt_add_ecmo   ; burnin = 100)

split_method = LeaveFutureK(; K=1, minimum = 475) # Total number of subjects is 678 so would like to keep #70% (475) as training dataset
split_by = BySubject()

cv_method = PSISCrossvalidation(;
split_method=split_method,
split_by=split_by,
)
cv_1cpt = crossvalidate(one_cmt_add, cv_method);

However I am hitting the following error:

cv_1cpt = crossvalidate(one_cmt_add, cv_method);
ERROR: ArgumentError: Invalid input for `log_ratios` (array is empty).
Stacktrace:
[1] _check_input_validity_psis(log_ratios::Array{Float64, 3})
@ ParetoSmooth /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/ParetoSmooth/A6x5U/src/ImportanceSampling.jl:325
[2] psis(log_ratios::Array{Float64, 3}; source::String, r_eff::Vector{Float64}, calc_ess::Bool, skip_checks::Bool)
@ ParetoSmooth /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/ParetoSmooth/A6x5U/src/ImportanceSampling.jl:137
@ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/bayes/crossvalidation/crossvalidation.jl:45
@ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/bayes/crossvalidation/crossvalidation.jl:86
@ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/bayes/crossvalidation/crossvalidation.jl:457
[6] top-level scope
@ ~/data/code/ECMO_UMMC/ECMO_BAYESIAN.jl:564

Another error when I use the following code:

BAYESIAN_one_cmt_prop_ecmo = @model begin
@param begin
ΞΈ ~ Constrained(
MvNormal(
[0.5,0.5],
[100 0.0
0.0 100]
),
lower = [0.0,0.0],
upper = [Inf, Inf],
init  = [0.5,0.5])
Οβ β Gamma(1,0.08)
Οβ β Gamma(1,0.08)
end
@random begin
Ξ·β ~ Normal(0, sqrt(Ξ©β))
Ξ·β ~ Normal(0, sqrt(Ξ©β))
end
@covariates WT ECMO_DAYS IS_CIRCUIT_CHANGE site T
@pre begin

# PK parameters

siteeffect = if site==2
Οβ
else
Οβ
end
#siteeffect_2 = site==2 ? Ο_add : 0.0

#tvKTV = site == 1 ? tvKTV1 : tvKTV2*(ECMO_DAYS/6)^tvecmo
#KTV    = tvKTV*exp(Ξ·KTV)
#TVCL = site == 1 ? tvCL1 * (WT/10)^0.75 : tvCL2* (WT/10)^0.75  * tvcircuit^IS_CIRCUIT_CHANGE
CL     = ΞΈ[1] * exp(Ξ·β[1])#*(1-exp(-KTV*t))
#TVVc = site == 1 ? tvVc1 * (WT/10)^1 : tvVc2 * (WT/10)^1
Vc     = ΞΈ[2] * exp(Ξ·β[1])
end
@dynamics Central1

@derived begin
CP = @. Central / Vc
CONC = @. Normal(CP, sqrt(CP^2*siteeffect))

end
end

one_cmt_prop_ecmo   = @time fit(BAYESIAN_one_cmt_prop_ecmo, estimation, init_params(BAYESIAN_one_cmt_prop_ecmo), Pumas.BayesMCMC(nsamples=6000,

one_cmt_prop = truncate(one_cmt_prop_ecmo, burnin = 100)

split_method = LeaveFutureK(; K=1, minimum = 475)
split_by = BySubject()
cv_method = PSISCrossvalidation(;
split_method=split_method,
split_by=split_by,
)

cv_1cpt_2 = crossvalidate(one_cmt_prop, cv_method);

I get the following error

cv_1cpt_2 = crossvalidate(one_cmt_prop, cv_method);
ERROR: DomainError with NaN:
Normal: the condition Ο >= zero(Ο) is not satisfied.
Stacktrace:
[1] #366
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Distributions/7iOJp/src/univariate/continuous/normal.jl:37 [inlined]
[2] check_args
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Distributions/7iOJp/src/utils.jl:89 [inlined]
[3] #Normal#365
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Distributions/7iOJp/src/univariate/continuous/normal.jl:37 [inlined]
[4] Normal
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Distributions/7iOJp/src/univariate/continuous/normal.jl:37 [inlined]
[7] getindex
[8] macro expansion
[9] macro expansion
@ ./simdloop.jl:77 [inlined]
[10] copyto!
[11] copyto!
[12] copy
[13] materialize
[14] (::Serialization.__deserialized_types__.var"##1771")(_pre::Pumas.Pre{Nothing, Serialization.__deserialized_types__.var"##1769"{Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}, Float64, Float64, Pumas.Promote{Float64}, Float64, Float64}}, _sol::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.Pre{Nothing, Serialization.__deserialized_types__.var"##1769"{Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}, Float64, Float64, Pumas.Promote{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.Pre{Nothing, Serialization.__deserialized_types__.var"##1769"{Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}, Float64, Float64, Pumas.Promote{Float64}, Float64, Float64}}}}, _obstimes::Vector{Float64}, _subject::Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}, _param::NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ), Tuple{Vector{Float64}, Float64, Float64, Float64, Float64}}, _random::NamedTuple{(:Ξ·β, :Ξ·β), Tuple{Float64, Float64}})
@ Main /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/dsl/model_macro.jl:1543
[15] #_derived#317
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/models/model_api.jl:1280 [inlined]
[16] #conditional_nll#496
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/likelihoods.jl:862 [inlined]
[17] _penalized_conditional_nll(model::PumasModel{(ΞΈ = 2, Ξ©β = 1, Ξ©β = 1, Οβ = 1, Οβ = 1), 2, (:Central,), ParamSet{NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ), Tuple{Constrained{FullNormal, VectorDomain{Vector{Float64}, Vector{TransformVariables.Infinity{true}}, Vector{Float64}}}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}}}}, Serialization.__deserialized_types__.var"##1766", Pumas.TimeDispatcher{Serialization.__deserialized_types__.var"##1767", Serialization.__deserialized_types__.var"##1768"}, Nothing, Serialization.__deserialized_types__.var"##1770", Central1, Serialization.__deserialized_types__.var"##1771", Serialization.__deserialized_types__.var"##1772", Nothing}, subject::Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}, param::NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ), Tuple{Vector{Float64}, Float64, Float64, Float64, Float64}}, vrandeffsorth::Vector{Float64}, diffeq_options::NamedTuple{(), Tuple{}})
@ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/likelihoods.jl:920
[18] #498
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/likelihoods.jl:998 [inlined]
[19] (::NLSolversBase.var"#95#101"{NLSolversBase.InplaceObjective{Nothing, Nothing, Pumas.var"#498#499"{PumasModel{(ΞΈ = 2, Ξ©β = 1, Ξ©β = 1, Οβ = 1, Οβ = 1), 2, (:Central,), ParamSet{NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ), Tuple{Constrained{FullNormal, VectorDomain{Vector{Float64}, Vector{TransformVariables.Infinity{true}}, Vector{Float64}}}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}}}}, Serialization.__deserialized_types__.var"##1766", Pumas.TimeDispatcher{Serialization.__deserialized_types__.var"##1767", Serialization.__deserialized_types__.var"##1768"}, Nothing, Serialization.__deserialized_types__.var"##1770", Central1, Serialization.__deserialized_types__.var"##1771", Serialization.__deserialized_types__.var"##1772", Nothing}, NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ), Tuple{Vector{Float64}, Float64, Float64, Float64, Float64}}, LaplaceI, NamedTuple{(), Tuple{}}, TransformVariables.TransformTuple{NamedTuple{(:Ξ·β, :Ξ·β), Tuple{Pumas.NormalTransform{Normal{Float64}}, Pumas.NormalTransform{Normal{Float64}}}}}, ParamSet{NamedTuple{(:Ξ·β, :Ξ·β), Tuple{Normal{Float64}, Normal{Float64}}}}, Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, Nothing, Nothing}, Float64})(x::Vector{Float64})
@ NLSolversBase /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/NLSolversBase/cfJrN/src/objective_types/incomplete.jl:112
[20] value!!(obj::NLSolversBase.TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}, x::Vector{Float64})
@ NLSolversBase /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/NLSolversBase/cfJrN/src/interface.jl:9
[21] value!
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/NLSolversBase/cfJrN/src/interface.jl:28 [inlined]
[22] update_state!(d::NLSolversBase.TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}, state::Optim.NewtonTrustRegionState{Vector{Float64}, Float64, Vector{Float64}}, method::Optim.NewtonTrustRegion{Float64})
@ Optim /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Optim/Zq1jM/src/multivariate/solvers/second_order/newton_trust_region.jl:344
[23] optimize(d::NLSolversBase.TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, method::Optim.NewtonTrustRegion{Float64}, options::Optim.Options{Float64, Nothing}, state::Optim.NewtonTrustRegionState{Vector{Float64}, Float64, Vector{Float64}})
@ Optim /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Optim/Zq1jM/src/multivariate/optimize/optimize.jl:54
[24] optimize
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Optim/Zq1jM/src/multivariate/optimize/optimize.jl:36 [inlined]
[25] _orth_empirical_bayes!(vrandeffsorth::Vector{Float64}, m::PumasModel{(ΞΈ = 2, Ξ©β = 1, Ξ©β = 1, Οβ = 1, Οβ = 1), 2, (:Central,), ParamSet{NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ), Tuple{Constrained{FullNormal, VectorDomain{Vector{Float64}, Vector{TransformVariables.Infinity{true}}, Vector{Float64}}}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}}}}, Serialization.__deserialized_types__.var"##1766", Pumas.TimeDispatcher{Serialization.__deserialized_types__.var"##1767", Serialization.__deserialized_types__.var"##1768"}, Nothing, Serialization.__deserialized_types__.var"##1770", Central1, Serialization.__deserialized_types__.var"##1771", Serialization.__deserialized_types__.var"##1772", Nothing}, _subject::Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}, param::NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ), Tuple{Vector{Float64}, Float64, Float64, Float64, Float64}}, approx::LaplaceI, diffeq_options::NamedTuple{(), Tuple{}})
@ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/likelihoods.jl:1007
[26] _orth_empirical_bayes(m::PumasModel{(ΞΈ = 2, Ξ©β = 1, Ξ©β = 1, Οβ = 1, Οβ = 1), 2, (:Central,), ParamSet{NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ), Tuple{Constrained{FullNormal, VectorDomain{Vector{Float64}, Vector{TransformVariables.Infinity{true}}, Vector{Float64}}}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}}}}, Serialization.__deserialized_types__.var"##1766", Pumas.TimeDispatcher{Serialization.__deserialized_types__.var"##1767", Serialization.__deserialized_types__.var"##1768"}, Nothing, Serialization.__deserialized_types__.var"##1770", Central1, Serialization.__deserialized_types__.var"##1771", Serialization.__deserialized_types__.var"##1772", Nothing}, subject::Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}, param::NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ), Tuple{Vector{Float64}, Float64, Float64, Float64, Float64}}, approx::LaplaceI, diffeq_options::NamedTuple{(), Tuple{}})
@ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/likelihoods.jl:938
[27] #marginal_nll#504
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/likelihoods.jl:1208 [inlined]
[28] (::Pumas.var"#985#987"{Vector{Float64}, BySubject{LaplaceI}, Pumas.var"#984#986"{TransformVariables.TransformTuple{NamedTuple{(:Ξ·β, :Ξ·β), Tuple{Pumas.NormalTransform{Normal{Float64}}, Pumas.NormalTransform{Normal{Float64}}}}}}, NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ), Tuple{Vector{Float64}, Float64, Float64, Float64, Float64}}, NamedTuple{(), Tuple{}}, Vector{Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, PumasModel{(ΞΈ = 2, Ξ©β = 1, Ξ©β = 1, Οβ = 1, Οβ = 1), 2, (:Central,), ParamSet{NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ), Tuple{Constrained{FullNormal, VectorDomain{Vector{Float64}, Vector{TransformVariables.Infinity{true}}, Vector{Float64}}}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}}}}, Serialization.__deserialized_types__.var"##1766", Pumas.TimeDispatcher{Serialization.__deserialized_types__.var"##1767", Serialization.__deserialized_types__.var"##1768"}, Nothing, Serialization.__deserialized_types__.var"##1770", Central1, Serialization.__deserialized_types__.var"##1771", Serialization.__deserialized_types__.var"##1772", Nothing}})(si::Int64)
@ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/bayes/crossvalidation/datasplit.jl:489
[29] iterate
@ ./generator.jl:47 [inlined]
[30] collect_to!
@ ./array.jl:782 [inlined]
[31] collect_to_with_first!
@ ./array.jl:760 [inlined]
[32] _collect(c::UnitRange{Int64}, itr::Base.Generator{UnitRange{Int64}, Pumas.var"#985#987"{Vector{Float64}, BySubject{LaplaceI}, Pumas.var"#984#986"{TransformVariables.TransformTuple{NamedTuple{(:Ξ·β, :Ξ·β), Tuple{Pumas.NormalTransform{Normal{Float64}}, Pumas.NormalTransform{Normal{Float64}}}}}}, NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ), Tuple{Vector{Float64}, Float64, Float64, Float64, Float64}}, NamedTuple{(), Tuple{}}, Vector{Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, PumasModel{(ΞΈ = 2, Ξ©β = 1, Ξ©β = 1, Οβ = 1, Οβ = 1), 2, (:Central,), ParamSet{NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ), Tuple{Constrained{FullNormal, VectorDomain{Vector{Float64}, Vector{TransformVariables.Infinity{true}}, Vector{Float64}}}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}}}}, Serialization.__deserialized_types__.var"##1766", Pumas.TimeDispatcher{Serialization.__deserialized_types__.var"##1767", Serialization.__deserialized_types__.var"##1768"}, Nothing, Serialization.__deserialized_types__.var"##1770", Central1, Serialization.__deserialized_types__.var"##1771", Serialization.__deserialized_types__.var"##1772", Nothing}}}, #unused#::Base.EltypeUnknown, isz::Base.HasShape{1})
@ Base ./array.jl:754
[33] collect_similar(cont::UnitRange{Int64}, itr::Base.Generator{UnitRange{Int64}, Pumas.var"#985#987"{Vector{Float64}, BySubject{LaplaceI}, Pumas.var"#984#986"{TransformVariables.TransformTuple{NamedTuple{(:Ξ·β, :Ξ·β), Tuple{Pumas.NormalTransform{Normal{Float64}}, Pumas.NormalTransform{Normal{Float64}}}}}}, NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ), Tuple{Vector{Float64}, Float64, Float64, Float64, Float64}}, NamedTuple{(), Tuple{}}, Vector{Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, PumasModel{(ΞΈ = 2, Ξ©β = 1, Ξ©β = 1, Οβ = 1, Οβ = 1), 2, (:Central,), ParamSet{NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ), Tuple{Constrained{FullNormal, VectorDomain{Vector{Float64}, Vector{TransformVariables.Infinity{true}}, Vector{Float64}}}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}}}}, Serialization.__deserialized_types__.var"##1766", Pumas.TimeDispatcher{Serialization.__deserialized_types__.var"##1767", Serialization.__deserialized_types__.var"##1768"}, Nothing, Serialization.__deserialized_types__.var"##1770", Central1, Serialization.__deserialized_types__.var"##1771", Serialization.__deserialized_types__.var"##1772", Nothing}}})
@ Base ./array.jl:653
[34] map(f::Function, A::UnitRange{Int64})
@ Base ./abstractarray.jl:2867
@ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/bayes/crossvalidation/datasplit.jl:479
[36] #965
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/bayes/crossvalidation/datasplit.jl:419 [inlined]
[37] iterate
@ ./generator.jl:47 [inlined]
@ Base ./array.jl:782
@ Base ./array.jl:760

Or when I when I change the cross validation for the above code to

split_method_1 = LeaveFutureK(; K=10, minimum = 475)
split_by_1 = BySubject(;marginal = nothing)

cv_method_2 = PSISCrossvalidation(;
split_method=split_method_1,
split_by=split_by_1,
)

cv_1cpt_2 = crossvalidate(one_cmt_prop, cv_method_2);

I receive the following error

cv_1cpt_2 = crossvalidate(one_cmt_prop, cv_method_2);
ERROR: ArgumentError: Invalid input for `log_ratios` (array is empty).
Stacktrace:
[1] _check_input_validity_psis(log_ratios::Array{Float64, 3})
@ ParetoSmooth /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/ParetoSmooth/A6x5U/src/ImportanceSampling.jl:325
[2] psis(log_ratios::Array{Float64, 3}; source::String, r_eff::Vector{Float64}, calc_ess::Bool, skip_checks::Bool)
@ ParetoSmooth /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/ParetoSmooth/A6x5U/src/ImportanceSampling.jl:137
@ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/bayes/crossvalidation/crossvalidation.jl:45
@ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/bayes/crossvalidation/crossvalidation.jl:86
@ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/bayes/crossvalidation/crossvalidation.jl:457
[6] top-level scope
@ ~/data/code/ECMO_UMMC/ECMO_BAYESIAN.jl:577

Switch to LeaveK instead of LeaveFutureK and remove the minimum. LeaveFutureK and BySubject assume you have more subjects than the value of minimum = 475. Because it will incrementally add subjects and evaluate the likelihood of the next one starting with 475 subjects in your case. LeaveK will just do leave-K-out cross-validation.

Also generally you would like K to be small with LeaveK. In LeaveK, K is the number of subjects left out in each run so that depends on the number of subjects but K = 1 is a common choice. Alternatively, you can use KFold instead which will figure out how many subjects to leave out each time to do K folds. K in this case has a different meaning which is the number of cross-validation folds.

Thanks @mohamed82008 for the reply. I tried your suggestion

split_method_1 = LeaveK(; K=1)
split_by_1 = BySubject(;)
cv_method_2 = PSISCrossvalidation(;
split_method=split_method_1,
split_by=split_by_1,
)

cv_1cpt = crossvalidate(one_cmt_add, cv_method_2);

But I am getting this error

cv_1cpt = crossvalidate(one_cmt_add, cv_method_2);
ERROR: DomainError with Invalid input for `log_ratios` (contains NaN  or inf values).:

Stacktrace:
[1] _check_input_validity_psis(log_ratios::Array{Float64, 3})
@ ParetoSmooth /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/ParetoSmooth/A6x5U/src/ImportanceSampling.jl:323
[2] psis(log_ratios::Array{Float64, 3}; source::String, r_eff::Vector{Float64}, calc_ess::Bool, skip_checks::Bool)
@ ParetoSmooth /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/ParetoSmooth/A6x5U/src/ImportanceSampling.jl:137
@ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/bayes/crossvalidation/crossvalidation.jl:45
@ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/bayes/crossvalidation/crossvalidation.jl:86
@ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/bayes/crossvalidation/crossvalidation.jl:462
[6] top-level scope
@ ~/data/code/ECMO_UMMC/ECMO_BAYESIAN.jl:576

Can I see the model? And does this happen with simulated data as well?

BAYESIAN_one_cmt_add_ecmo = @model begin
@param begin
#tvKTV1 ~ LogNormal(log(0.4), 1.0)
#tvKTV2 ~ LogNormal(log(0.1), 1.0)
#tvecmo ~ Normal(-0.7, 1.0)
#tvCL1 ~ LogNormal(log(0.9), 1.0)
tvCL2 ~ LogNormal(log(0.5), 1.0)
#tvcircuit ~ LogNormal(log(1.2), 1.0)
#tvVc1 ~ LogNormal(log(1.8), 0.1)
tvVc2 ~ LogNormal(log(0.5), 1.0)
#Ο_prop_1 ~ truncated(Cauchy(0.18, 5), 0.0, Inf)
Ο_add_2 ~ truncated(Cauchy(10000, 5), 0.0, Inf)
Ο_add_3 ~ truncated(Cauchy(10000, 5), 0.0, Inf)
#Ο_add ~ truncated(Cauchy(3500, 5), 0.0, Inf)
#ΟΒ²KTV ~ truncated(Normal(1.5, 10), 0.0, Inf)
Ο ~ MvLogNormal(ones(2), Diagonal([1.0,1.0]))

end

@random begin
Ξ· ~ MvNormal(Ο)
end
@covariates WT ECMO_DAYS IS_CIRCUIT_CHANGE site T
@pre begin

# PK parameters

siteeffect = if site==2
else
end
#siteeffect_2 = site==2 ? Ο_add : 0.0

#tvKTV = site == 1 ? tvKTV1 : tvKTV2*(ECMO_DAYS/6)^tvecmo
#KTV    = tvKTV*exp(Ξ·KTV)
#TVCL = site == 1 ? tvCL1 * (WT/10)^0.75 : tvCL2* (WT/10)^0.75  * tvcircuit^IS_CIRCUIT_CHANGE
CL     = tvCL2 * exp(Ξ·[1])#*(1-exp(-KTV*t))
#TVVc = site == 1 ? tvVc1 * (WT/10)^1 : tvVc2 * (WT/10)^1
Vc     = tvVc2 * exp(Ξ·[2])
end

@dynamics Central1
@derived begin
CP = @. Central / Vc
CONC = @. Normal(CP, sqrt(siteeffect))
end
end

I havenβt tried any simulation yet but the vpc function works fine ( it depends on simulation output)

Hi @mohamed82008, Just as a follow up I tried the PSISCrossvalidation in another model ( just changed the priors and the error model as following:

BAYESIAN_one_cmt_prop_ecmo = @model begin
@param begin
ΞΈ ~ Constrained(
MvNormal(
[0.5,0.5],
[100 0.0
0.0 100]
),
lower = [0.0,0.0],
upper = [Inf, Inf],
init  = [0.5,0.5])
Οβ β Gamma(1,0.08)
Οβ β Gamma(1,0.08)
end
@random begin
Ξ·β ~ Normal(0, sqrt(Ξ©β))
Ξ·β ~ Normal(0, sqrt(Ξ©β))
end
@covariates WT ECMO_DAYS IS_CIRCUIT_CHANGE site T
@pre begin

# PK parameters

siteeffect = if site==2
Οβ
else
Οβ
end
#siteeffect_2 = site==2 ? Ο_add : 0.0

#tvKTV = site == 1 ? tvKTV1 : tvKTV2*(ECMO_DAYS/6)^tvecmo
#KTV    = tvKTV*exp(Ξ·KTV)
#TVCL = site == 1 ? tvCL1 * (WT/10)^0.75 : tvCL2* (WT/10)^0.75  * tvcircuit^IS_CIRCUIT_CHANGE
CL     = ΞΈ[1] * exp(Ξ·β[1])#*(1-exp(-KTV*t))
#TVVc = site == 1 ? tvVc1 * (WT/10)^1 : tvVc2 * (WT/10)^1
Vc     = ΞΈ[2] * exp(Ξ·β[1])
end
@dynamics Central1

@derived begin
CP = @. Central / Vc
CONC = @. Normal(CP, sqrt(CP^2*siteeffect))

end
end

one_cmt_prop_ecmo   = @time fit(BAYESIAN_one_cmt_prop_ecmo, estimation, init_params(BAYESIAN_one_cmt_prop_ecmo), Pumas.BayesMCMC(nsamples=6000,

one_cmt_prop = truncate(one_cmt_prop_ecmo , burnin = 100)

split_method_1 = LeaveK(; K=2)

split_by_1 = BySubject(marginal = nothing)

cv_method_2 = PSISCrossvalidation(;
split_method=split_method_1,
split_by=split_by_1,
)

cv_1cpt_2 = crossvalidate(one_cmt_prop, cv_method_2);

However, I am getting the following warning

cv_1cpt_2 = crossvalidate(one_cmt_prop, cv_method_2);
β Warning: All the data points were removed as part of the filtering of high and NaN Pareto shape parameter values. This is a sign that the Pareto smoothed importance sampling algorithm for crossvalidation failed. Please consider using a different data splitting method or a different crossvalidation method.
β @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/bayes/crossvalidation/crossvalidation.jl:488

and the elpd function produce

julia> elpd(cv_1cpt_2)
NaN

Please use the ExactCrossvalidation method when PSISCrossvalidation fails. PSIS is sensitive to problems with a few data points. I would still like to dig into the first error so if you can reproduce it with a simulated data set and post the script here, that would be helpful.

2 Likes

Thanks so much @mohamed82008, I will follow-up with the simulated data crossvalidation

Hi @mohamed82008 I simulated data and run PSISCrossvalidation using a simulated dataset as following

df_simulation = @chain df begin
@rsubset :evid .==1
end

observations=[:CONC],
id=:ID,
time=:TIME,
evid=:evid,
mdv = :mdv,
amt=:AMT,
rate = :RATE,
cmt = :cmt,
covariates = [:WT]
)

simulated_data = DataFrame([simobs(mod, modfit.data[i], fparam,
obstimes = minimum(modfit.data[i].time):2:maximum(modfit.data[i].time))
for i in 1:length(modfit.data)])

simulated_data = @chain simulated_data begin
@rsubset :evid β [2]
@rtransform :cmt = 1
end

observations=[:CONC],
id=:id,
time=:time,
evid=:evid,
amt=:amt,
rate = :rate,
cmt = :cmt,
covariates = [:WT]
)

BAYESIAN_one_cmt_prop_ecmo = @model begin
@param begin
ΞΈ ~ Constrained(
MvNormal(
[0.5,0.5],
[100 0.0
0.0 100]
),
lower = [0.0,0.0],
upper = [Inf, Inf],
init  = [0.5,0.5])
Οβ β Gamma(1,0.01)
end
@random begin
Ξ·β ~ Normal(0, sqrt(Ξ©β))
Ξ·β ~ Normal(0, sqrt(Ξ©β))
end
@covariates WT
@pre begin

CL     = ΞΈ[1] * exp(Ξ·β[1])* (WT/10)^0.75
Vc     = ΞΈ[2] * exp(Ξ·β[1])* (WT/10)^1
end
@dynamics Central1

@derived begin
CP = @. Central / Vc
CONC = @. Normal(CP, sqrt(CP^2*Οβ))

end
end

one_cmt_prop_ecmo   = @time fit(BAYESIAN_one_cmt_prop_ecmo, estimation, init_params(BAYESIAN_one_cmt_prop_ecmo), Pumas.BayesMCMC(nsamples=2000,

res_t = Pumas.truncate(one_cmt_prop_ecmo; burnin=1000)

split_method_1 = LeaveK(; K=1)
split_by_1 = BySubject(marginal = nothing)

cv_method_2 = PSISCrossvalidation(;
split_method=split_method_1,
split_by=split_by_1,
)

cv_1cpt = crossvalidate(res_t, cv_method_2);

However I am still getting the error

cv_1cpt = crossvalidate(one_cmt_add, cv_method);
β Warning: All the data points were removed as part of the filtering of high and NaN Pareto shape parameter values. This is a sign that the Pareto smoothed importance sampling algorithm for crossvalidation failed. Please consider using a different data splitting method or a different crossvalidation method.
β @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/bayes/crossvalidation/crossvalidation.jl:488

Also When I use ExactCrossvalidation it takes excessive amount of time ( 3 days and no output yet). Given that I have around 700 subjects and the way the ExactCrossvalidation works (running Bayesian estimation for each validation round), I think this method will take an excessive amount of time. Is there a way to increase the speed of the function ?

Thanks

If PSISCrossvalidation shows this warning, it means that its result cannot be trusted. But it could also point towards failed inference. The way PSIS works with LeaveK(K = 1) and BySubject for example is by re-weighting the MCMC samples after removing each subject. The weights are then smoothed which is an important step I will revisit later. If you only have a few data points (subjects in this case), the posterior distribution resulting from removing even a single data point may be significantly different from the posterior distribution when using all of the data points. In the weighting scheme, this means that most weights will be very close to 0 and the effective sample size of the new weighted samples will be significantly lower.

In the weight smoothing part of the algorithm, the weights are fitted to a generalised Pareto distribution where the shape parameter k of the distribution gives you a diagnostic that tells you when the estimates coming from the re-weighted posterior cannot be trusted. A value of k > 0.7 is usually a sign that we shouldnβt trust this set of weights. Sometimes the smoothing step may also fail and result in NaNs. In Pumas, when removing a data point gives a high shape parameter k > 0.7 or NaN, we skip removing this data point and move on to the next one. If removing any of the points leads to a high shape parameter k, then you get the warning above. This typically happens if you have a few influential data points but given that you have 700 subjects, itβs surprising to see that PSIS is failing. This is why I think itβs possible your inference may have failed in some way. Did you check all of the diagnostics and predictive plots before doing cross-validation? Just as an experiment, try setting the pareto_shape_threshold keyword argument of PSISCrossvalidation to a high value instead of the default 0.7. See if this gives you an answer.

Regarding ExactCrossvalidation, you can use ensemblealg = EnsenbleDistributed() and many cores to crunch through the runs. You may also want to use KFold(K = 10) instead of LeaveK(K.= 1) because the latter will need to run 700 MCMC runs while the former will only need 10 MCMC runs. Let me know if any of these suggestions works for you.

1 Like

@ahmed.salem please try PSIS with marginal = LaplaceI() as well.

Thanks so much @mohamed82008 for responding. KFold crossvalidation worked faster by setting KFold(K = 5) which took around 6 hours to finish the crossvalidation ( I subset the dataset and used information from only 188 subjects). I think setting KFold(K = 10) will take longer time so I will try to set ensemblealg = EnsenbleDistributed()to speed up the process. Do you recommend specific vPCUs and GB memory while initiating the Pumas task that is optimal for ensemblealg = EnsenbleDistributed() ?

I tried conditional (by setting BySubject(marginal = nothing) and marginal crossvalidation (by setting BySubject(marginal = LaplaceI()). The elpd function works fine with the conditional crossvalidation, however I am receiving the following error from marginal crossvalidation BySubject(marginal = LaplaceI())

elpd(cv_1cpt_4)
[ Info: Computing out-of-sample loglikelihoods.
ERROR: DomainError with NaN:
Normal: the condition Ο >= zero(Ο) is not satisfied.
Stacktrace:
[1] #366
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Distributions/7iOJp/src/univariate/continuous/normal.jl:37 [inlined]
[2] check_args
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Distributions/7iOJp/src/utils.jl:89 [inlined]
[3] #Normal#365
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Distributions/7iOJp/src/univariate/continuous/normal.jl:37 [inlined]
[4] Normal
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Distributions/7iOJp/src/univariate/continuous/normal.jl:37 [inlined]
[7] getindex
[8] macro expansion
[9] macro expansion
@ ./simdloop.jl:77 [inlined]
[10] copyto!
[11] copyto!
[12] copy
[13] materialize
[14] (::Serialization.__deserialized_types__.var"#70#77")(_pre::Pumas.Pre{Nothing, Serialization.__deserialized_types__.var"#68#75"{Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}, Float64, Float64, Pumas.Promote{Float64}, Float64, Float64, Float64}}, _sol::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.Pre{Nothing, Serialization.__deserialized_types__.var"#68#75"{Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}, Float64, Float64, Pumas.Promote{Float64}, 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.Pre{Nothing, Serialization.__deserialized_types__.var"#68#75"{Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}, Float64, Float64, Pumas.Promote{Float64}, Float64, Float64, Float64}}}}, _obstimes::Vector{Float64}, _subject::Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}, _param::NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ, :Οβ), Tuple{Vector{Float64}, Float64, Float64, Float64, Float64, Float64}}, _random::NamedTuple{(:Ξ·β, :Ξ·β), Tuple{Float64, Float64}})
@ Main /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/dsl/model_macro.jl:1543
[15] #_derived#317
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/models/model_api.jl:1280 [inlined]
[16] #conditional_nll#496
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/likelihoods.jl:862 [inlined]
[17] _penalized_conditional_nll(model::PumasModel{(ΞΈ = 2, Ξ©β = 1, Ξ©β = 1, Οβ = 1, Οβ = 1, Οβ = 1), 2, (:Central,), ParamSet{NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ, :Οβ), Tuple{Constrained{FullNormal, VectorDomain{Vector{Float64}, Vector{TransformVariables.Infinity{true}}, Vector{Float64}}}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}}}}, Serialization.__deserialized_types__.var"#65#72", Pumas.TimeDispatcher{Serialization.__deserialized_types__.var"#66#73", Serialization.__deserialized_types__.var"#67#74"}, Nothing, Serialization.__deserialized_types__.var"#69#76", Central1, Serialization.__deserialized_types__.var"#70#77", Serialization.__deserialized_types__.var"#71#78", Nothing}, subject::Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}, param::NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ, :Οβ), Tuple{Vector{Float64}, Float64, Float64, Float64, Float64, Float64}}, vrandeffsorth::Vector{Float64}, diffeq_options::NamedTuple{(), Tuple{}})
@ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/likelihoods.jl:920
[18] #498
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/likelihoods.jl:998 [inlined]
[19] (::NLSolversBase.var"#95#101"{NLSolversBase.InplaceObjective{Nothing, Nothing, Pumas.var"#498#499"{PumasModel{(ΞΈ = 2, Ξ©β = 1, Ξ©β = 1, Οβ = 1, Οβ = 1, Οβ = 1), 2, (:Central,), ParamSet{NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ, :Οβ), Tuple{Constrained{FullNormal, VectorDomain{Vector{Float64}, Vector{TransformVariables.Infinity{true}}, Vector{Float64}}}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}}}}, Serialization.__deserialized_types__.var"#65#72", Pumas.TimeDispatcher{Serialization.__deserialized_types__.var"#66#73", Serialization.__deserialized_types__.var"#67#74"}, Nothing, Serialization.__deserialized_types__.var"#69#76", Central1, Serialization.__deserialized_types__.var"#70#77", Serialization.__deserialized_types__.var"#71#78", Nothing}, NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ, :Οβ), Tuple{Vector{Float64}, Float64, Float64, Float64, Float64, Float64}}, LaplaceI, NamedTuple{(), Tuple{}}, TransformVariables.TransformTuple{NamedTuple{(:Ξ·β, :Ξ·β), Tuple{Pumas.NormalTransform{Normal{Float64}}, Pumas.NormalTransform{Normal{Float64}}}}}, ParamSet{NamedTuple{(:Ξ·β, :Ξ·β), Tuple{Normal{Float64}, Normal{Float64}}}}, Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, Nothing, Nothing}, Float64})(x::Vector{Float64})
@ NLSolversBase /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/NLSolversBase/cfJrN/src/objective_types/incomplete.jl:112
[20] value!!(obj::NLSolversBase.TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}, x::Vector{Float64})
@ NLSolversBase /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/NLSolversBase/cfJrN/src/interface.jl:9
[21] value!
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/NLSolversBase/cfJrN/src/interface.jl:28 [inlined]
[22] update_state!(d::NLSolversBase.TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}, state::Optim.NewtonTrustRegionState{Vector{Float64}, Float64, Vector{Float64}}, method::Optim.NewtonTrustRegion{Float64})
@ Optim /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Optim/Zq1jM/src/multivariate/solvers/second_order/newton_trust_region.jl:344
[23] optimize(d::NLSolversBase.TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, method::Optim.NewtonTrustRegion{Float64}, options::Optim.Options{Float64, Nothing}, state::Optim.NewtonTrustRegionState{Vector{Float64}, Float64, Vector{Float64}})
@ Optim /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Optim/Zq1jM/src/multivariate/optimize/optimize.jl:54
[24] optimize
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Optim/Zq1jM/src/multivariate/optimize/optimize.jl:36 [inlined]
[25] _orth_empirical_bayes!(vrandeffsorth::Vector{Float64}, m::PumasModel{(ΞΈ = 2, Ξ©β = 1, Ξ©β = 1, Οβ = 1, Οβ = 1, Οβ = 1), 2, (:Central,), ParamSet{NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ, :Οβ), Tuple{Constrained{FullNormal, VectorDomain{Vector{Float64}, Vector{TransformVariables.Infinity{true}}, Vector{Float64}}}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}}}}, Serialization.__deserialized_types__.var"#65#72", Pumas.TimeDispatcher{Serialization.__deserialized_types__.var"#66#73", Serialization.__deserialized_types__.var"#67#74"}, Nothing, Serialization.__deserialized_types__.var"#69#76", Central1, Serialization.__deserialized_types__.var"#70#77", Serialization.__deserialized_types__.var"#71#78", Nothing}, _subject::Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}, param::NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ, :Οβ), Tuple{Vector{Float64}, Float64, Float64, Float64, Float64, Float64}}, approx::LaplaceI, diffeq_options::NamedTuple{(), Tuple{}})
@ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/likelihoods.jl:1007
[26] _orth_empirical_bayes(m::PumasModel{(ΞΈ = 2, Ξ©β = 1, Ξ©β = 1, Οβ = 1, Οβ = 1, Οβ = 1), 2, (:Central,), ParamSet{NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ, :Οβ), Tuple{Constrained{FullNormal, VectorDomain{Vector{Float64}, Vector{TransformVariables.Infinity{true}}, Vector{Float64}}}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}}}}, Serialization.__deserialized_types__.var"#65#72", Pumas.TimeDispatcher{Serialization.__deserialized_types__.var"#66#73", Serialization.__deserialized_types__.var"#67#74"}, Nothing, Serialization.__deserialized_types__.var"#69#76", Central1, Serialization.__deserialized_types__.var"#70#77", Serialization.__deserialized_types__.var"#71#78", Nothing}, subject::Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}, param::NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ, :Οβ), Tuple{Vector{Float64}, Float64, Float64, Float64, Float64, Float64}}, approx::LaplaceI, diffeq_options::NamedTuple{(), Tuple{}})
@ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/likelihoods.jl:938
[27] #marginal_nll#504
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/likelihoods.jl:1208 [inlined]
[28] (::Pumas.var"#1040#1048"{Vector{Float64}, Pumas.var"#1039#1047"{TransformVariables.TransformTuple{NamedTuple{(:Ξ·β, :Ξ·β), Tuple{Pumas.NormalTransform{Normal{Float64}}, Pumas.NormalTransform{Normal{Float64}}}}}}, Pumas.var"#1038#1046"{TransformVariables.TransformTuple{NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ, :Οβ), Tuple{Pumas.ElementArrayTransform{TransformVariables.ShiftedExp{true, Float64}, 1}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}}}}}, Vector{Float64}, NamedTuple{(), Tuple{}}, PumasModel{(ΞΈ = 2, Ξ©β = 1, Ξ©β = 1, Οβ = 1, Οβ = 1, Οβ = 1), 2, (:Central,), ParamSet{NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ, :Οβ), Tuple{Constrained{FullNormal, VectorDomain{Vector{Float64}, Vector{TransformVariables.Infinity{true}}, Vector{Float64}}}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}}}}, Serialization.__deserialized_types__.var"#65#72", Pumas.TimeDispatcher{Serialization.__deserialized_types__.var"#66#73", Serialization.__deserialized_types__.var"#67#74"}, Nothing, Serialization.__deserialized_types__.var"#69#76", Central1, Serialization.__deserialized_types__.var"#70#77", Serialization.__deserialized_types__.var"#71#78", Nothing}, Vector{Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, BySubject{LaplaceI}})(si::Int64)
@ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/bayes/crossvalidation/crossvalidation.jl:620
[29] macro expansion
@ ./reduce.jl:247 [inlined]
[30] macro expansion
@ ./simdloop.jl:77 [inlined]
[31] mapreduce_impl(f::Pumas.var"#1040#1048"{Vector{Float64}, Pumas.var"#1039#1047"{TransformVariables.TransformTuple{NamedTuple{(:Ξ·β, :Ξ·β), Tuple{Pumas.NormalTransform{Normal{Float64}}, Pumas.NormalTransform{Normal{Float64}}}}}}, Pumas.var"#1038#1046"{TransformVariables.TransformTuple{NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ, :Οβ), Tuple{Pumas.ElementArrayTransform{TransformVariables.ShiftedExp{true, Float64}, 1}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}}}}}, Vector{Float64}, NamedTuple{(), Tuple{}}, PumasModel{(ΞΈ = 2, Ξ©β = 1, Ξ©β = 1, Οβ = 1, Οβ = 1, Οβ = 1), 2, (:Central,), ParamSet{NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ, :Οβ), Tuple{Constrained{FullNormal, VectorDomain{Vector{Float64}, Vector{TransformVariables.Infinity{true}}, Vector{Float64}}}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}}}}, Serialization.__deserialized_types__.var"#65#72", Pumas.TimeDispatcher{Serialization.__deserialized_types__.var"#66#73", Serialization.__deserialized_types__.var"#67#74"}, Nothing, Serialization.__deserialized_types__.var"#69#76", Central1, Serialization.__deserialized_types__.var"#70#77", Serialization.__deserialized_types__.var"#71#78", Nothing}, Vector{Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, BySubject{LaplaceI}}, op::typeof(Base.add_sum), A::UnitRange{Int64}, ifirst::Int64, ilast::Int64, blksize::Int64)
@ Base ./reduce.jl:245
[32] mapreduce_impl
@ ./reduce.jl:259 [inlined]
[33] _mapreduce(f::Pumas.var"#1040#1048"{Vector{Float64}, Pumas.var"#1039#1047"{TransformVariables.TransformTuple{NamedTuple{(:Ξ·β, :Ξ·β), Tuple{Pumas.NormalTransform{Normal{Float64}}, Pumas.NormalTransform{Normal{Float64}}}}}}, Pumas.var"#1038#1046"{TransformVariables.TransformTuple{NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ, :Οβ), Tuple{Pumas.ElementArrayTransform{TransformVariables.ShiftedExp{true, Float64}, 1}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}}}}}, Vector{Float64}, NamedTuple{(), Tuple{}}, PumasModel{(ΞΈ = 2, Ξ©β = 1, Ξ©β = 1, Οβ = 1, Οβ = 1, Οβ = 1), 2, (:Central,), ParamSet{NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Οβ, :Οβ, :Οβ), Tuple{Constrained{FullNormal, VectorDomain{Vector{Float64}, Vector{TransformVariables.Infinity{true}}, Vector{Float64}}}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}}}}, Serialization.__deserialized_types__.var"#65#72", Pumas.TimeDispatcher{Serialization.__deserialized_types__.var"#66#73", Serialization.__deserialized_types__.var"#67#74"}, Nothing, Serialization.__deserialized_types__.var"#69#76", Central1, Serialization.__deserialized_types__.var"#70#77", Serialization.__deserialized_types__.var"#71#78", Nothing}, Vector{Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, BySubject{LaplaceI}}, op::typeof(Base.add_sum), #unused#::IndexLinear, A::UnitRange{Int64})
@ Base ./reduce.jl:417
[34] _mapreduce_dim(f::Function, op::Function, #unused#::Base._InitialValue, A::UnitRange{Int64}, #unused#::Colon)
@ Base ./reducedim.jl:330
[35] #mapreduce#731
@ ./reducedim.jl:322 [inlined]
[36] mapreduce
@ ./reducedim.jl:322 [inlined]
[37] #_sum#741
@ ./reducedim.jl:894 [inlined]
[38] _sum
@ ./reducedim.jl:894 [inlined]
[39] #sum#739
@ ./reducedim.jl:890 [inlined]
[40] sum
@ ./reducedim.jl:890 [inlined]
@ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/bayes/crossvalidation/crossvalidation.jl:609
[42] iterate
@ ./generator.jl:47 [inlined]

This error happens when I call elpd function using the crossvalidate output. But I do not receive any error while running the crossvalidate function`

Thanks @mohamed82008 for the suggestion. I looked into the diagnostic plots of the model that throw the above error when using the PSISCrossvalidate() function and it seems that there is a high autocorrelation for one of the parameters as shown in the figure

The error still appear after setting pareto_shape_thresholdto high values (I tried setting to 1 and 10). and the error appears with both marginal = LaplaceI()and marginal = nothing

cv_1cpt = crossvalidate(one_cmt_add, cv_method_3);
β Warning: All the data points were removed as part of the filtering of high and NaN Pareto shape parameter values. This is a sign that the Pareto smoothed importance sampling algorithm for crossvalidation failed. Please consider using a different data splitting method or a different crossvalidation method.
β @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/bayes/crossvalidation/crossvalidation.jl:488

FYI information this is the model that I used

BAYESIAN_one_cmt_add_ecmo = @model begin
@param begin
#tvKTV1 ~ LogNormal(log(0.4), 1.0)
#tvKTV2 ~ LogNormal(log(0.1), 1.0)
#tvecmo ~ Normal(-0.7, 1.0)
#tvCL1 ~ LogNormal(log(0.9), 1.0)
tvCL2 ~ LogNormal(log(0.5), 1.0)
#tvcircuit ~ LogNormal(log(1.2), 1.0)
#tvVc1 ~ LogNormal(log(1.8), 0.1)
tvVc2 ~ LogNormal(log(0.5), 1.0)
#Ο_prop_1 ~ truncated(Cauchy(0.18, 5), 0.0, Inf)
Ο_add_2 ~ truncated(Normal(10000, 100), 0.0, Inf)
Ο_add_3 ~ truncated(Normal(10000, 100), 0.0, Inf)
#Ο_add ~ truncated(Cauchy(3500, 5), 0.0, Inf)
#ΟΒ²KTV ~ truncated(Normal(1.5, 10), 0.0, Inf)
Ο ~ Constrained(
MvNormal(zeros(2), Diagonal([100.0,100.0])),
lower=zeros(2),
upper=fill(Inf, 2),
init=ones(2),
)
end

@random begin
Ξ· ~ MvNormal(Ο)
end
@covariates WT ECMO_DAYS IS_CIRCUIT_CHANGE site T
@pre begin

# PK parameters

siteeffect = if site==2
else
end
#siteeffect_2 = site==2 ? Ο_add : 0.0

#tvKTV = site == 1 ? tvKTV1 : tvKTV2*(ECMO_DAYS/6)^tvecmo
#KTV    = tvKTV*exp(Ξ·KTV)
#TVCL = site == 1 ? tvCL1 * (WT/10)^0.75 : tvCL2* (WT/10)^0.75  * tvcircuit^IS_CIRCUIT_CHANGE
CL     = tvCL2 * exp(Ξ·[1])#*(1-exp(-KTV*t))
#TVVc = site == 1 ? tvVc1 * (WT/10)^1 : tvVc2 * (WT/10)^1
Vc     = tvVc2 * exp(Ξ·[2])
end

@dynamics Central1
@derived begin
CP = @. Central / Vc
CONC = @. Normal(CP, sqrt(siteeffect))
end
end

I tried another model which could fit the data better as following:

BAYESIAN_one_cmt_comb_ecmo_not_all_wt_KTV = @model begin
@param begin
ΞΈ ~ Constrained(
MvNormal(
[0.1, 0.8,2.6],
[100 0.0 0.0;
0.0 100 0.0;
0.0 0.0 100]
),
lower = [0.0,0.0,0.0],
upper = [Inf, Inf,  Inf],
init  = [0.1,0.8,2.6])
Οβ β Gamma(1,0.08)
Οβ β Gamma(1,0.08)
Οβ β Gamma(1,3000)
#Οβ β Gamma(1,3000)
end
@random begin
Ξ·β ~ Normal(0, sqrt(Ξ©β))
Ξ·β ~ Normal(0, sqrt(Ξ©β))
Ξ·β ~ Normal(0, sqrt(Ξ©β))
end
@covariates WT ECMO_DAYS IS_CIRCUIT_CHANGE site T
@pre begin

# PK parameters

siteeffect = if site==2
Οβ
else
Οβ
end
siteeffect_2 = site==2 ? Οβ : 0.0

#tvKTV = site == 1 ? tvKTV1 : tvKTV2*(ECMO_DAYS/6)^tvecmo
KTV    = ΞΈ[1] *exp(Ξ·β[1])
#TVCL = site == 1 ? tvCL1 * (WT/10)^0.75 : tvCL2* (WT/10)^0.75  * tvcircuit^IS_CIRCUIT_CHANGE
CL     = ΞΈ[2] * exp(Ξ·β[1])* (WT/10)^0.75 *(1-exp(-KTV*t+1e-5))
#TVVc = site == 1 ? tvVc1 * (WT/10)^1 : tvVc2 * (WT/10)^1
Vc     = ΞΈ[3] * exp(Ξ·β[1])* (WT/10)^1
end
@dynamics Central1

@derived begin
CP = @. Central / Vc
CONC = @. Normal(CP, sqrt(CP^2*siteeffect + siteeffect_2))

end
end

All the diagnostic plots and statistics show signs of successful convergence

I tried to run

cv_method = PSISCrossvalidation(;
split_method=KFold(; K=5,shuffle = false),
split_by=BySubject(marginal = LaplaceI())
)

But I am getting the following error

cross = crossvalidate(res_t,cv_method)
ERROR: DomainError with NaN:
Normal: the condition Ο >= zero(Ο) is not satisfied.
Stacktrace:
[1] #366
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Distributions/7iOJp/src/univariate/continuous/normal.jl:37 [inlined]
[2] check_args
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Distributions/7iOJp/src/utils.jl:89 [inlined]
[3] #Normal#365
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Distributions/7iOJp/src/univariate/continuous/normal.jl:37 [inlined]
[4] Normal
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Distributions/7iOJp/src/univariate/continuous/normal.jl:37 [inlined]
[7] getindex
[8] macro expansion
[9] macro expansion
@ ./simdloop.jl:77 [inlined]
[10] copyto!
[11] copyto!
[12] copy
[13] materialize
[14] (::Serialization.__deserialized_types__.var"##427")(_pre::Pumas.Pre{Nothing, Serialization.__deserialized_types__.var"##425"{Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}, Float64, Pumas.Promote{Float64}, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64}}, _sol::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.Pre{Nothing, Serialization.__deserialized_types__.var"##425"{Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}, Float64, Pumas.Promote{Float64}, Float64, Float64, Float64, Float64, Float64, 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.Pre{Nothing, Serialization.__deserialized_types__.var"##425"{Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}, Float64, Pumas.Promote{Float64}, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64}}}}, _obstimes::Vector{Float64}, _subject::Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}, _param::NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Ξ©β, :Οβ, :Οβ, :Οβ), Tuple{Vector{Float64}, Float64, Float64, Float64, Float64, Float64, Float64}}, _random::NamedTuple{(:Ξ·β, :Ξ·β, :Ξ·β), Tuple{Float64, Float64, Float64}})
@ Main /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/dsl/model_macro.jl:1543
[15] #_derived#317
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/models/model_api.jl:1280 [inlined]
[16] #conditional_nll#496
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/likelihoods.jl:862 [inlined]
[17] _penalized_conditional_nll(model::PumasModel{(ΞΈ = 3, Ξ©β = 1, Ξ©β = 1, Ξ©β = 1, Οβ = 1, Οβ = 1, Οβ = 1), 3, (:Central,), ParamSet{NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Ξ©β, :Οβ, :Οβ, :Οβ), Tuple{Constrained{FullNormal, VectorDomain{Vector{Float64}, Vector{TransformVariables.Infinity{true}}, Vector{Float64}}}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}}}}, Serialization.__deserialized_types__.var"##422", Pumas.TimeDispatcher{Serialization.__deserialized_types__.var"##423", Serialization.__deserialized_types__.var"##424"}, Nothing, Serialization.__deserialized_types__.var"##426", Central1, Serialization.__deserialized_types__.var"##427", Serialization.__deserialized_types__.var"##428", Nothing}, subject::Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}, param::NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Ξ©β, :Οβ, :Οβ, :Οβ), Tuple{Vector{Float64}, Float64, Float64, Float64, Float64, Float64, Float64}}, vrandeffsorth::Vector{Float64}, diffeq_options::NamedTuple{(), Tuple{}})
@ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/likelihoods.jl:920
[18] #498
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/likelihoods.jl:998 [inlined]
[19] (::NLSolversBase.var"#95#101"{NLSolversBase.InplaceObjective{Nothing, Nothing, Pumas.var"#498#499"{PumasModel{(ΞΈ = 3, Ξ©β = 1, Ξ©β = 1, Ξ©β = 1, Οβ = 1, Οβ = 1, Οβ = 1), 3, (:Central,), ParamSet{NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Ξ©β, :Οβ, :Οβ, :Οβ), Tuple{Constrained{FullNormal, VectorDomain{Vector{Float64}, Vector{TransformVariables.Infinity{true}}, Vector{Float64}}}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}}}}, Serialization.__deserialized_types__.var"##422", Pumas.TimeDispatcher{Serialization.__deserialized_types__.var"##423", Serialization.__deserialized_types__.var"##424"}, Nothing, Serialization.__deserialized_types__.var"##426", Central1, Serialization.__deserialized_types__.var"##427", Serialization.__deserialized_types__.var"##428", Nothing}, NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Ξ©β, :Οβ, :Οβ, :Οβ), Tuple{Vector{Float64}, Float64, Float64, Float64, Float64, Float64, Float64}}, LaplaceI, NamedTuple{(), Tuple{}}, TransformVariables.TransformTuple{NamedTuple{(:Ξ·β, :Ξ·β, :Ξ·β), Tuple{Pumas.NormalTransform{Normal{Float64}}, Pumas.NormalTransform{Normal{Float64}}, Pumas.NormalTransform{Normal{Float64}}}}}, ParamSet{NamedTuple{(:Ξ·β, :Ξ·β, :Ξ·β), Tuple{Normal{Float64}, Normal{Float64}, Normal{Float64}}}}, Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, Nothing, Nothing}, Float64})(x::Vector{Float64})
@ NLSolversBase /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/NLSolversBase/cfJrN/src/objective_types/incomplete.jl:112
[20] value!!(obj::NLSolversBase.TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}, x::Vector{Float64})
@ NLSolversBase /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/NLSolversBase/cfJrN/src/interface.jl:9
[21] value!
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/NLSolversBase/cfJrN/src/interface.jl:28 [inlined]
[22] update_state!(d::NLSolversBase.TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}, state::Optim.NewtonTrustRegionState{Vector{Float64}, Float64, Vector{Float64}}, method::Optim.NewtonTrustRegion{Float64})
@ Optim /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Optim/Zq1jM/src/multivariate/solvers/second_order/newton_trust_region.jl:344
[23] optimize(d::NLSolversBase.TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, method::Optim.NewtonTrustRegion{Float64}, options::Optim.Options{Float64, Nothing}, state::Optim.NewtonTrustRegionState{Vector{Float64}, Float64, Vector{Float64}})
@ Optim /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Optim/Zq1jM/src/multivariate/optimize/optimize.jl:54
[24] optimize
@ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Optim/Zq1jM/src/multivariate/optimize/optimize.jl:36 [inlined]
[25] _orth_empirical_bayes!(vrandeffsorth::Vector{Float64}, m::PumasModel{(ΞΈ = 3, Ξ©β = 1, Ξ©β = 1, Ξ©β = 1, Οβ = 1, Οβ = 1, Οβ = 1), 3, (:Central,), ParamSet{NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Ξ©β, :Οβ, :Οβ, :Οβ), Tuple{Constrained{FullNormal, VectorDomain{Vector{Float64}, Vector{TransformVariables.Infinity{true}}, Vector{Float64}}}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}, Gamma{Float64}}}}, Serialization.__deserialized_types__.var"##422", Pumas.TimeDispatcher{Serialization.__deserialized_types__.var"##423", Serialization.__deserialized_types__.var"##424"}, Nothing, Serialization.__deserialized_types__.var"##426", Central1, Serialization.__deserialized_types__.var"##427", Serialization.__deserialized_types__.var"##428", Nothing}, _subject::Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE, :site, :T), Tuple{Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}, param::NamedTuple{(:ΞΈ, :Ξ©β, :Ξ©β, :Ξ©β, :Οβ, :Οβ, :Οβ), Tuple{Vector{Float64}, Float64, Float64, Float64, Float64, Float64, Float64}}, approx::LaplaceI, diffeq_options::NamedTuple{(), Tuple{}})

But when I try to set the following

cv_method = PSISCrossvalidation(;
split_method=KFold(; K=5,shuffle = false),
split_by=BySubject(marginal = nothing)
)

The function works fine but when I try to call elpd() to the output, I receive the following output

cv_method = PSISCrossvalidation(;
split_method=KFold(; K=5,shuffle = false),
split_by=BySubject(marginal = nothing)
)

cross = crossvalidate(res_t,cv_method)

elpd(cross)
NaN

I am not sure why elpd() function does not provide output

The error you are getting with Laplace is happening because the empirical Bayes estimate (EBE) calculation is failing. This is not your fault but the model may be incompatible with LaplaceI. Try BySubject(marginal = LLQuad(imaxiters=1000)) instead which will use a quadrature integration method instead of the Laplace method to estimate the marginal likelihoods.