Pumas_be errors

1

Hello,

I’m running simulations over several scenarios, and I’m getting errors when running pumas_be on one specific scenario. I’m unable to identify the cause of the errors.

I’m running the following lines of code on the data I’ve pasted at the bottom:

cmax_output = pumas_be(test, endpoint = :cmax)
auc_output = pumas_be(test, endpoint = :auc)
auc0_28_output = pumas_be(test, endpoint = :auc0_28)
auc7_28_output = pumas_be(test, endpoint = :auc7_28)

auc_output and auc7_28_output run, but I get the following errors for cmax_output and auc0_28 output:

julia> cmax_output = pumas_be(test5, endpoint = :cmax)
ERROR: ArgumentError: reducing over an empty collection is not allowed
Stacktrace:
  [1] _empty_reduce_error()
    @ Base ./reduce.jl:301
  [2] mapreduce_empty(f::Function, op::Base.BottomRF{typeof(Base.add_sum)}, T::Type)
    @ Base ./reduce.jl:344
  [3] reduce_empty(op::Base.MappingRF{Bioequivalence.var"#6#11", Base.BottomRF{typeof(Base.add_sum)}}, #unused#::Type{Float64})
    @ Base ./reduce.jl:331
  [4] reduce_empty(op::Base.FilteringRF{Bioequivalence.var"#7#12", Base.MappingRF{Bioequivalence.var"#6#11", Base.BottomRF{typeof(Base.add_sum)}}}, #unused#::Type{Float64})
    @ Base ./reduce.jl:332
  [5] reduce_empty_iter
    @ ./reduce.jl:357 [inlined]
  [6] reduce_empty_iter
    @ ./reduce.jl:356 [inlined]
  [7] foldl_impl
    @ ./reduce.jl:49 [inlined]
  [8] mapfoldl_impl
    @ ./reduce.jl:44 [inlined]
  [9] #mapfoldl#244
    @ ./reduce.jl:162 [inlined]
 [10] mapfoldl
    @ ./reduce.jl:162 [inlined]
 [11] #mapreduce#248
    @ ./reduce.jl:289 [inlined]
 [12] mapreduce
    @ ./reduce.jl:289 [inlined]
 [13] #sum#251
    @ ./reduce.jl:503 [inlined]
 [14] sum
    @ ./reduce.jl:503 [inlined]
 [15] #sum#252
    @ ./reduce.jl:532 [inlined]
 [16] sum(a::Base.Generator{Base.Iterators.Filter{Bioequivalence.var"#7#12", Vector{Float64}}, Bioequivalence.var"#6#11"})
    @ Base ./reduce.jl:532
 [17] satterthwaite(model::MixedModels.LinearMixedModel{Float64}, L::Matrix{Float64})
    @ Bioequivalence /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Bioequivalence/AEgIQ/src/satterthwaite.jl:95
 [18] walds_tests(model::MixedModels.LinearMixedModel{Float64})
    @ Bioequivalence /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Bioequivalence/AEgIQ/src/Bioequivalence.jl:299
 [19] pumas_be(data::DataFrame; endpoint::Symbol, logtransformed::Bool, reference_scale::Float64, cv_max::Float64, id::String, sequence::String, period::String, nonparametric::Bool, reml::Bool)
    @ Bioequivalence /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Bioequivalence/AEgIQ/src/Bioequivalence.jl:1025
 [20] top-level scope
    @ ~/data/code/padagis_mibe/modeling/partial replicate design simulations (parallel absorption - zo, Erlang) 3.jl:676

julia> auc0_28_output = pumas_be(test5, endpoint = :auc0_28)
ERROR: DomainError with -13.175565702861608:
FDist: the condition ν2 > zero(ν2) is not satisfied.
Stacktrace:
  [1] #269
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Distributions/7iOJp/src/univariate/continuous/fdist.jl:30 [inlined]
  [2] check_args
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Distributions/7iOJp/src/utils.jl:89 [inlined]
  [3] _
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Distributions/7iOJp/src/univariate/continuous/fdist.jl:30 [inlined]
  [4] #FDist#272
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Distributions/7iOJp/src/univariate/continuous/fdist.jl:35 [inlined]
  [5] #FDist#274
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Distributions/7iOJp/src/univariate/continuous/fdist.jl:37 [inlined]
  [6] FDist(ν1::Int64, ν2::Float64)
    @ Distributions /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Distributions/7iOJp/src/univariate/continuous/fdist.jl:37
  [7] satterthwaite(model::MixedModels.LinearMixedModel{Float64}, L::Matrix{Float64})
    @ Bioequivalence /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Bioequivalence/AEgIQ/src/satterthwaite.jl:98
  [8] walds_tests(model::MixedModels.LinearMixedModel{Float64})
    @ Bioequivalence /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Bioequivalence/AEgIQ/src/Bioequivalence.jl:299
  [9] pumas_be(data::DataFrame; endpoint::Symbol, logtransformed::Bool, reference_scale::Float64, cv_max::Float64, id::String, sequence::String, period::String, nonparametric::Bool, reml::Bool)
    @ Bioequivalence /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Bioequivalence/AEgIQ/src/Bioequivalence.jl:1025
 [10] top-level scope
    @ ~/data/code/padagis_mibe/modeling/partial replicate design simulations (parallel absorption - zo, Erlang) 3.jl:678

There are only 4 subjects with an RTT/TRR replicate design, but a scenario with 40% WSV ran without problems. This scenario has no WSV. Could I get some help identifying why I’m getting these errors?

The data causing the errors is below:

id	rep	period	cmax	auc	auc0_28	auc7_28	Nsubj	input_tr_ratio	BOV	sequence	formulation	p_extrap_auc0_28
1	1	1	26.95198647	11.05709908	10.19941708	2.135242623	4	1	0.1	RTT	R	7.756844729
1	1	2	38.42458639	12.4261716	11.40120656	2.407952729	4	1	0.1	RTT	T	8.248437836
1	1	3	33.437928	12.62193615	11.69189361	2.548529996	4	1	0.1	RTT	T	7.368461747
2	1	1	25.95525597	12.24250729	11.41183829	3.997105946	4	1	0.1	RTT	R	6.785121543
2	1	2	27.55211256	13.32263612	11.51540158	4.025919439	4	1	0.1	RTT	T	13.56514227
2	1	3	22.27539908	16.13506787	13.79141672	5.1853434	4	1	0.1	RTT	T	14.52520168
3	1	1	19.61582162	6.753974148	6.527068701	2.333494383	4	1	0.1	TRR	T	3.359584169
3	1	2	21.31397657	7.995812583	7.951315931	2.655606377	4	1	0.1	TRR	R	0.556499435
3	1	3	20.20298612	9.744124406	9.640481691	4.010845008	4	1	0.1	TRR	R	1.06364318
4	1	1	27.4860378	14.62849581	13.5428536	5.179847392	4	1	0.1	TRR	T	7.42142066
4	1	2	31.37180564	15.16105866	12.51296622	4.044817363	4	1	0.1	TRR	R	17.46640854
4	1	3	24.53625692	15.47109261	13.63232752	4.709957563	4	1	0.1	TRR	R	11.88516637
2

Hi Donald,

Thanks for reporting this issue. I’m able to reproduce it. The dataset exercises a code path related to the calculation of the Satterthwaite degrees of freedom that hasn’t yet been exercised by our current tests. I have prepared a fix which will be available in the upcoming release.

3

Hi @andreasnoack,

I’m now able to get pumas_be to run, but now I get an error with assess_be when the criterion is either FDA or FDA_NarrowTherapeuticIndex.

Here is the pumas_be output:

image
image1041×291 7.27 KB

And here is the error:

ERROR: MethodError: no method matching Chisq(::Missing)
Closest candidates are:
  Chisq(::Integer; check_args) at /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Distributions/YQQXX/src/univariate/continuous/chisq.jl:33
  Chisq(::Real; check_args) at /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Distributions/YQQXX/src/univariate/continuous/chisq.jl:28
Stacktrace:
  [1] update_rsabe_theta!(obj::ReferenceScaledAverageBioequivalance{StatsModels.TableRegressionModel{GLM.LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}}}, 𝜃::Float64)
    @ Bioequivalence /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Bioequivalence/XifjK/src/Bioequivalence.jl:903
  [2] assess_be(criterion::AverageBioequivalenceWithExpandingLimits, result::BioequivalenceStudy{NamedTuple{(:treatment, :sequence, :period), Tuple{DataFrame, DataFrame, DataFrame}}, NamedTuple{(:RRT, :RTR, :TRR), Tuple{Int64, Int64, Int64}}, MixedModels.LinearMixedModel{Float64}, NamedTuple{(:within_subject_variability, :rsabe), Tuple{DataFrame, ReferenceScaledAverageBioequivalance{StatsModels.TableRegressionModel{GLM.LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}}}}}, Float64})
    @ Bioequivalence /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Bioequivalence/XifjK/src/Bioequivalence.jl:1844
  [3] _broadcast_getindex_evalf
    @ ./broadcast.jl:670 [inlined]
  [4] _broadcast_getindex
    @ ./broadcast.jl:643 [inlined]
  [5] getindex
    @ ./broadcast.jl:597 [inlined]
  [6] copyto_nonleaf!(dest::Vector{BitVector}, bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(assess_be), Tuple{Base.Broadcast.Extruded{Vector{Bioequivalence.BioequivalenceCriterion}, Tuple{Bool}, Tuple{Int64}}, Base.RefValue{BioequivalenceStudy{NamedTuple{(:treatment, :sequence, :period), Tuple{DataFrame, DataFrame, DataFrame}}, NamedTuple{(:RRT, :RTR, :TRR), Tuple{Int64, Int64, Int64}}, MixedModels.LinearMixedModel{Float64}, NamedTuple{(:within_subject_variability, :rsabe), Tuple{DataFrame, ReferenceScaledAverageBioequivalance{StatsModels.TableRegressionModel{GLM.LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}}}}}, Float64}}}}, iter::Base.OneTo{Int64}, state::Int64, count::Int64)
    @ Base.Broadcast ./broadcast.jl:1055
  [7] copy
    @ ./broadcast.jl:907 [inlined]
  [8] materialize
    @ ./broadcast.jl:860 [inlined]
  [9] assess_be(criteria::Vector{Bioequivalence.BioequivalenceCriterion}, result::BioequivalenceStudy{NamedTuple{(:treatment, :sequence, :period), Tuple{DataFrame, DataFrame, DataFrame}}, NamedTuple{(:RRT, :RTR, :TRR), Tuple{Int64, Int64, Int64}}, MixedModels.LinearMixedModel{Float64}, NamedTuple{(:within_subject_variability, :rsabe), Tuple{DataFrame, ReferenceScaledAverageBioequivalance{StatsModels.TableRegressionModel{GLM.LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}}}}}, Float64})
    @ Bioequivalence /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Bioequivalence/XifjK/src/Bioequivalence.jl:1943
 [10] top-level scope
    @ REPL[13]:1

The error is the same for FDA_NarrowTherapeuticIndex.

Also, how do I access the assess_be results? For example, for pumas_be outputs, I can do cmax_ouput.result.LB.

4

Thanks for reporting this. I can reproduce the issue for the RRT_TRT_TRR design. We’ll look into it.

Regarding your second question then will also look into improving the documention. Meanwhile, the simplest solution might be for you to browse the returned structs in the variable explorer window in VSCode. Then you can hopefully find the fields you are interested in.

5

Hi @donaldlee3, thanks for reporting the issue.
I wanted to add that the FDA_NarrowTherapeuticIndex criteria includes the WithinSubjectVariabilityRatio criterion which is only applicable for fully replicated crossover designs. For NTI, partially replicated designs are not accepted as one cannot estimate the within subject variability of the test reference to assess the ratio. For designs such as RRT|RTR|TRR using the FDA will apply the criteria for non-NTI which include when applicable reference scaling for highly variable drugs.

In terms of accessing the result of the assessment, we will update the documentation.
Using,

assessment = assess_be(FDA, result)

Allows to examine each criterion.

assessment.criteria # provides a vector of each criterion

Through

assessment.results

It provides a BitMatrix where each row is a comparison (e.g. R|T, R|S) and the columns are indicators on whether passed the criterion.