Error for BIC, AIC calculation from a Bayesian output

Estimation bayesian
1

Hello,
I am trying to calculate BIC and AIC using the following code

split_method= LeaveK(K = 1)

split_by= ByObservation(allsubjects = true)

res_t = Pumas.truncate(res; burnin=100)

pl = loglikelihood( res_t ; split_method , split_by)

Pumas.aic(pl)
Pumas.aicc(pl)
Pumas.bic(pl)
Pumas.waic(pl)

The loglikelihood function runs ok with the following output

pl = loglikelihood( res_t ; split_method , split_by)
[ Info: Computing in-sample loglikelihoods.
1×5900×4 Array{Float64, 3}:
[:, :, 1] =
 -3372.21  -3372.21  -3325.67  -3325.67  -3315.43  -3329.12  -3320.48  -3320.48  -3334.26  -3325.41  -3344.78  …  -3334.55  -3325.02  -3366.97  -3388.51  -3393.29  -3433.4  -3399.69  -3363.97  -3377.28  -3371.55

[:, :, 2] =
 -3379.97  -3379.97  -3346.1  -3346.1  -3346.72  -3357.6  -3357.6  -3345.27  -3338.86  -3344.66  -3306.59  …  -3350.65  -3375.88  -3355.34  -3360.2  -3390.24  -3348.77  -3337.15  -3383.26  -3344.77  -3310.35

[:, :, 3] =
 -3376.41  -3378.06  -3365.15  -3380.94  -3365.8  -3389.8  -3389.8  -3369.57  -3388.12  -3388.12  -3371.54  …  -3330.06  -3407.95  -3455.92  -3449.55  -3407.8  -3470.26  -3470.26  -3497.54  -3487.15  -3503.74

[:, :, 4] =
 -3331.21  -3331.21  -3333.46  -3333.46  -3316.36  -3366.62  -3366.62  -3405.8  -3374.23  -3369.58  -3369.58  …  -3360.76  -3350.42  -3364.6  -3380.73  -3409.23  -3358.75  -3330.57  -3336.66  -3306.26  -3322.87

However I am getting the following error with Pumas.aic(), Pumas.aicc(), Pumas.bic() or Pumas.waic()

julia> aic(pl)
ERROR: MethodError: no method matching aic(::Array{Float64, 3})
Closest candidates are:
  aic(::StatsAPI.StatisticalModel) at /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/StatsAPI/y7Ydc/src/statisticalmodel.jl:189
  aic(::Pumas.AbstractFittedPumasModel) at /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/diagnostics.jl:1287
  aic(::Pumas.PointwiseLogLikelihood) at /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/bayes/crossvalidation/crossvalidation.jl:98
Stacktrace:
 [1] top-level scope
   @ ~/data/code/ECMO_UMMC/NON_ECMO.jl:127

Please let me know how to fix this error. Thanks

2

Can you try Pumas.cv_pointwise_loglikelihood instead of loglikelihood?

1 Like
3

This works but for the following combination of split_method and split_by it gives errors. example of cases when it does not work:

split_method= KFold(K = 5)

split_by= ByObservation(allsubjects = true)


pl = Pumas.cv_pointwise_loglikelihood(res_t,split_method,split_by)

I get the following error

pl = Pumas.cv_pointwise_loglikelihood(res_t,split_method,split_by)
ERROR: BoundsError: attempt to access 8×1 Matrix{Float64} at index [7:9, 1:1]
Stacktrace:
  [1] throw_boundserror(A::Matrix{Float64}, I::Tuple{UnitRange{Int64}, Base.Slice{Base.OneTo{Int64}}})
    @ Base ./abstractarray.jl:691
  [2] checkbounds
    @ ./abstractarray.jl:656 [inlined]
  [3] view
    @ ./subarray.jl:177 [inlined]
  [4] #973
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/bayes/crossvalidation/datasplit.jl:446 [inlined]
  [5] iterate
    @ ./generator.jl:47 [inlined]
  [6] collect_to!(dest::Vector{Float64}, itr::Base.Generator{Vector{UnitRange{Int64}}, Pumas.var"#973#980"{Matrix{Float64}}}, offs::Int64, st::Int64)
    @ Base ./array.jl:782
  [7] collect_to_with_first!
    @ ./array.jl:760 [inlined]
  [8] _collect(c::Vector{UnitRange{Int64}}, itr::Base.Generator{Vector{UnitRange{Int64}}, Pumas.var"#973#980"{Matrix{Float64}}}, #unused#::Base.EltypeUnknown, isz::Base.HasShape{1})
    @ Base ./array.jl:754
  [9] collect_similar
    @ ./array.jl:653 [inlined]
 [10] map
    @ ./abstractarray.jl:2867 [inlined]
 [11] #972
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/bayes/crossvalidation/datasplit.jl:445 [inlined]
 [12] mapreduce_impl(f::Pumas.var"#972#979"{Vector{UnitRange{Int64}}}, op::typeof(+), A::Vector{Matrix{Float64}}, ifirst::Int64, ilast::Int64, blksize::Int64)
    @ Base ./reduce.jl:244
 [13] mapreduce_impl
    @ ./reduce.jl:259 [inlined]
 [14] _mapreduce(f::Pumas.var"#972#979"{Vector{UnitRange{Int64}}}, op::typeof(+), #unused#::IndexLinear, A::Vector{Matrix{Float64}})
    @ Base ./reduce.jl:417
 [15] _mapreduce_dim(f::Function, op::Function, #unused#::Base._InitialValue, A::Vector{Matrix{Float64}}, #unused#::Colon)
    @ Base ./reducedim.jl:330
 [16] #mapreduce#731
    @ ./reducedim.jl:322 [inlined]
 [17] mapreduce(f::Function, op::Function, A::Vector{Matrix{Float64}})
    @ Base ./reducedim.jl:322
 [18] (::Pumas.var"#971#978"{Pumas.BayesMCMCResults{BayesMCMC{AdvancedHMC.GeneralisedNoUTurn{Float64}, Float64, Int64, Int64, Int64, Int64, Int64, Float64, EnsembleThreads, Bool, Bool, Bool, Bool, TaskLocalRNG, NamedTuple{(:abstol, :reltol, :alg), Tuple{Float64, Float64, Rodas5P{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}}, Bool, Int64, NamedTuple{(), Tuple{}}}, Vector{Pumas.ThreadedBayesLogDensity{PumasModel{(θ = 5, Ω₁ = 1, Ω₂ = 1, Ω₃ = 1, Ω₄ = 1, σ₁ = 1), 4, (: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"#127#137", Pumas.TimeDispatcher{Serialization.__deserialized_types__.var"#128#138", Serialization.__deserialized_types__.var"#130#140"}, Nothing, Serialization.__deserialized_types__.var"#132#142", ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, ODEFunction{false, ModelingToolkit.ODEFunctionClosure{Serialization.__deserialized_types__.var"#133#143", Serialization.__deserialized_types__.var"#134#144"}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Serialization.__deserialized_types__.var"#135#145", Serialization.__deserialized_types__.var"#136#146", ModelingToolkit.ODESystem}, NamedTuple{(:θ, :Ω₁, :Ω₂, :Ω₃, :Ω₄, :σ₁), Tuple{Vector{Float64}, Float64, Float64, Float64, Float64, Float64}}, 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}}}, Vector{Vector{Float64}}, Vector{Vector{Float64}}, Vector{ForwardDiff.GradientConfig{ForwardDiff.Tag{Pumas.Tag, Float64}, Float64, 7, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Pumas.Tag, Float64}, Float64, 7}}}}, Vector{DiffResults.MutableDiffResult{1, Float64, Tuple{Vector{Float64}}}}, NamedTuple{(:abstol, :reltol, :alg), Tuple{Float64, Float64, Rodas5P{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}}}}, Vector{Vector{Vector{Float64}}}, Vector{AbstractMCMC.SamplingStats}, Nothing}, ByObservation{true}, Vector{UnitRange{Int64}}})(vaug::Vector{Float64})
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/bayes/crossvalidation/datasplit.jl:444
 [19] iterate
    @ ./generator.jl:47 [inlined]
 [20] _collect(c::Vector{Vector{Float64}}, itr::Base.Generator{Vector{Vector{Float64}}, Pumas.var"#971#978"{Pumas.BayesMCMCResults{BayesMCMC{AdvancedHMC.GeneralisedNoUTurn{Float64}, Float64, Int64, Int64, Int64, Int64, Int64, Float64, EnsembleThreads, Bool, Bool, Bool, Bool, TaskLocalRNG, NamedTuple{(:abstol, :reltol, :alg), Tuple{Float64, Float64, Rodas5P{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}}, Bool, Int64, NamedTuple{(), Tuple{}}}, Vector{Pumas.ThreadedBayesLogDensity{PumasModel{(θ = 5, Ω₁ = 1, Ω₂ = 1, Ω₃ = 1, Ω₄ = 1, σ₁ = 1), 4, (: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"#127#137", Pumas.TimeDispatcher{Serialization.__deserialized_types__.var"#128#138", Serialization.__deserialized_types__.var"#130#140"}, Nothing, Serialization.__deserialized_types__.var"#132#142", ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, ODEFunction{false, ModelingToolkit.ODEFunctionClosure{Serialization.__deserialized_types__.var"#133#143", Serialization.__deserialized_types__.var"#134#144"}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Serialization.__deserialized_types__.var"#135#145", Serialization.__deserialized_types__.var"#136#146", ModelingToolkit.ODESystem}, NamedTuple{(:θ, :Ω₁, :Ω₂, :Ω₃, :Ω₄, :σ₁), Tuple{Vector{Float64}, Float64, Float64, Float64, Float64, Float64}}, 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}}}, Vector{Vector{Float64}}, Vector{Vector{Float64}}, Vector{ForwardDiff.GradientConfig{ForwardDiff.Tag{Pumas.Tag, Float64}, Float64, 7, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Pumas.Tag, Float64}, Float64, 7}}}}, Vector{DiffResults.MutableDiffResult{1, Float64, Tuple{Vector{Float64}}}}, NamedTuple{(:abstol, :reltol, :alg), Tuple{Float64, Float64, Rodas5P{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}}}}, Vector{Vector{Vector{Float64}}}, Vector{AbstractMCMC.SamplingStats}, Nothing}, ByObservation{true}, Vector{UnitRange{Int64}}}}, #unused#::Base.EltypeUnknown, isz::Base.HasShape{1})

Another case is

split_method_3= LeaveFutureK(K = 2)

split_by= ByObservation(allsubjects = true)
pl = Pumas.cv_pointwise_loglikelihood(res_t,split_method_3,split_by)

It gives the following output

julia> split_method_3= LeaveFutureK(K = 2)
LeaveFutureK(2, 2, false)

julia> pl = Pumas.cv_pointwise_loglikelihood(res_t,split_method_3,split_by)
Pointwise log likelihood
  Type: LeaveFutureK(2, 2, false) - ByObservation - allsubjects = true
  Number of groups: 0
  Number of samples per chain: 1900
  Number of chains: 4

but when I run Pumas.bic(pl) it gives the following error:

Pumas.bic(pl)
ERROR: ArgumentError: reducing over an empty collection is not allowed
Stacktrace:
  [1] _empty_reduce_error()
    @ Base ./reduce.jl:301
  [2] reduce_empty(op::Function, #unused#::Type{Float64})
    @ Base ./reduce.jl:311
  [3] mapreduce_empty(#unused#::typeof(identity), op::Function, T::Type)
    @ Base ./reduce.jl:345
  [4] reduce_empty(op::Base.MappingRF{typeof(identity), typeof(max)}, #unused#::Type{Float64})
    @ Base ./reduce.jl:331
  [5] reduce_empty_iter
    @ ./reduce.jl:357 [inlined]
  [6] mapreduce_empty_iter(f::Function, op::Function, itr::Array{Float64, 3}, ItrEltype::Base.HasEltype)
    @ Base ./reduce.jl:353
  [7] _mapreduce(f::typeof(identity), op::typeof(max), #unused#::IndexLinear, A::Array{Float64, 3})
    @ Base ./reduce.jl:402
  [8] _mapreduce_dim
    @ ./reducedim.jl:330 [inlined]
  [9] #mapreduce#731
    @ ./reducedim.jl:322 [inlined]
 [10] mapreduce
    @ ./reducedim.jl:322 [inlined]
 [11] #_maximum#749
    @ ./reducedim.jl:894 [inlined]
 [12] _maximum
    @ ./reducedim.jl:894 [inlined]
 [13] #_maximum#748
    @ ./reducedim.jl:893 [inlined]
 [14] _maximum
    @ ./reducedim.jl:893 [inlined]
 [15] #maximum#746
    @ ./reducedim.jl:889 [inlined]
 [16] maximum
    @ ./reducedim.jl:889 [inlined]
 [17] bic(pl::Pumas.PointwiseLogLikelihood{Array{Float64, 3}, Pumas.SplitData{UnitRange{Int64}, Vector{UnitRange{Int64}}, Vector{UnitRange{Int64}}, UnitRange{Int64}, LeaveFutureK, ByObservation{true}, Vector{UnitRange{Int64}}, Vector{UnitRange{Int64}}}, Pumas.BayesMCMCResults{BayesMCMC{AdvancedHMC.GeneralisedNoUTurn{Float64}, Float64, Int64, Int64, Int64, Int64, Int64, Float64, EnsembleThreads, Bool, Bool, Bool, Bool, TaskLocalRNG, NamedTuple{(), Tuple{}}, Bool, Int64, NamedTuple{(), Tuple{}}}, Vector{Pumas.ThreadedBayesLogDensity{PumasModel{(θ = 4, Ω₁ = 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"#35#42", Pumas.TimeDispatcher{Serialization.__deserialized_types__.var"#36#43", Serialization.__deserialized_types__.var"#37#44"}, Nothing, Serialization.__deserialized_types__.var"#39#46", Central1, Serialization.__deserialized_types__.var"#40#47", Serialization.__deserialized_types__.var"#41#48", Nothing}, NamedTuple{(:θ, :Ω₁, :Ω₂, :σ₁, :σ₂, :σ₃), Tuple{Vector{Float64}, Float64, Float64, Float64, Float64, Float64}}, 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}}}, Vector{Vector{Float64}}, Vector{Vector{Float64}}, Vector{ForwardDiff.GradientConfig{ForwardDiff.Tag{Pumas.Tag, Float64}, Float64, 6, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Pumas.Tag, Float64}, Float64, 6}}}}, Vector{DiffResults.MutableDiffResult{1, Float64, Tuple{Vector{Float64}}}}, NamedTuple{(), Tuple{}}}}, Vector{Vector{Vector{Float64}}}, Vector{AbstractMCMC.SamplingStats}, Nothing}})
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/bayes/crossvalidation/crossvalidation.jl:117
 [18] top-level scope
    @ ~/data/code/ECMO_UMMC/Total_BAYESIAN.jl:578
4

Do you have unequal number of observations per subject?

5

Yes
median 12 ranging from 2 to 35

6

There is a known bug in ByObservation when the number of observations is different per subject. Please use BySubjectfor now.

1 Like