Skip to content

Commit c4dadd4

Browse files
PharmCatPharmCat
authored
Merge pull request #16 from PharmCat/dev
Dev
2 parents 2383992 + 7a4a6ff commit c4dadd4

5 files changed

Lines changed: 127 additions & 95 deletions

File tree

‎Project.toml‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@ uuid = "00e1d80e-504f-4d7b-946f-0aea2e7abfe3"
33
keywords = ["bioequivalence"]
44
desc = "Bioequivalence"
55
authors = ["Vladimir Arnautov <mail@pharmcat.net>"]
6-
version = "0.2.10"
6+
version = "0.2.11"
77

88
[deps]
99
CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597"
@@ -16,14 +16,14 @@ MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
1616
StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
1717

1818
[compat]
19-
CategoricalArrays = "0.10"
19+
CategoricalArrays = "0.10, 1"
2020
DataFrames = "1"
21-
Distributions = "0.21, 0.22, 0.23, 0.24, 0.25"
21+
Distributions = "0.22, 0.23, 0.24, 0.25"
2222
GLM = "1.8, 1.9"
2323
Metida = "0.15, 0.16, 0.17"
2424
MetidaBase = "0.12, 0.13, 0.14"
2525
MixedModels = "4"
26-
StatsBase = "0.30, 0.31, 0.32, 0.33, 0.34"
26+
StatsBase = "0.31, 0.32, 0.33, 0.34"
2727
julia = "1.8"
2828

2929
[extras]

‎docs/Project.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -21,4 +21,4 @@ MixedModels = "4"
2121
PrettyTables = "1, 2"
2222
StatsBase = "≥0.33"
2323
Distributions = "0.24, 0.25"
24-
CategoricalArrays = "0.10"
24+
CategoricalArrays = "0.10, 1"

‎docs/src/api.md‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,11 @@ MetidaBioeq.estimate
1616
MetidaBioeq.result
1717
```
1818

19+
### MetidaBioeq.variance
20+
```@docs
21+
MetidaBioeq.variance
22+
```
23+
1924
### MetidaBioeq.makeseq
2025
```@docs
2126
MetidaBioeq.makeseq

‎src/bioequivalence.jl‎

Lines changed: 115 additions & 90 deletions
Original file line numberDiff line numberDiff line change
@@ -14,9 +14,12 @@ function nomissing(data, cols::AbstractVector)
1414
true
1515
end
1616

17-
function functional_term(f, arg_expr...)
18-
expr = Expr(:call, Symbol(f), arg_expr...)
19-
eval(:(@formula 0 ~ $expr)).rhs
17+
function functional_term(f, args...)
18+
FunctionTerm(
19+
f,
20+
term.(args),
21+
Expr(:call, Symbol(f), args...)
22+
)
2023
end
2124

2225
"""
@@ -336,6 +339,40 @@ function makeseq(data;
336339
map(x-> subjdict[x], getcol(data, subject))
337340
end
338341

342+
343+
function models_create(be, estimator, design = "auto")
344+
if design == "parallel"
345+
if be.logt
346+
return [FormulaTerm(term(i), term(be.formulation)) for i in be.vars]
347+
else
348+
return [FormulaTerm(functional_term(log, i), term(be.formulation)) for i in be.vars]
349+
end
350+
end
351+
352+
353+
if estimator == "glm"
354+
if be.logt
355+
return [FormulaTerm(term(i), term(be.formulation) + term(be.period) + term(be.sequence) + term(be.subject)) for i in be.vars]
356+
else
357+
return [FormulaTerm(functional_term(log, i), term(be.formulation) + term(be.period) + term(be.sequence) + term(be.subject)) for i in be.vars]
358+
end
359+
elseif estimator == "met"
360+
if be.logt
361+
return [FormulaTerm(term(i), term(be.formulation) + term(be.period) + term(be.sequence)) for i in be.vars]
362+
else
363+
return [FormulaTerm(functional_term(log, i), term(be.formulation) + term(be.period) + term(be.sequence)) for i in be.vars]
364+
end
365+
elseif estimator == "mm"
366+
if be.logt
367+
return [FormulaTerm(term(i), term(be.formulation) + term(be.period) + term(be.sequence) + functional_term(|, 1, be.subject)) for i in be.vars]
368+
else
369+
return [FormulaTerm(functional_term(log, i), term(be.formulation) + term(be.period) + term(be.sequence) + functional_term(|, 1, be.subject)) for i in be.vars]
370+
end
371+
else
372+
error("Unknown estimator!")
373+
end
374+
end
375+
339376
"""
340377
estimate(be; estimator = "auto", method = "auto", supresswarn = false, alpha = 0.05)
341378
@@ -388,7 +425,7 @@ EMA: [GUIDELINE ON THE INVESTIGATION OF BIOEQUIVALENCE](https://www.ema.europa.e
388425
EMA: [GUIDELINE ON THE INVESTIGATION OF BIOEQUIVALENCE, Annex I](https://www.ema.europa.eu/en/documents/other/31-annex-i-statistical-analysis-methods-compatible-ema-bioequivalence-guideline_en.pdf)
389426
390427
"""
391-
function estimate(be; estimator = "auto", method = "auto", supresswarn = false, alpha = 0.05)
428+
function estimate(be; estimator = "auto", method = "auto", supresswarn = false, alpha = 0.05, rsabe = :none, ntid = false)
392429

393430
length(be.formulations) > 2 && error("More than 2 formulations not supported yet")
394431
design = be.design
@@ -432,17 +469,8 @@ function estimate(be; estimator = "auto", method = "auto", supresswarn = false,
432469
if method != "P" && !supresswarn @warn("Method not P (parallel), for parallel simple GLM model will be used!") end
433470
estimator = "glm"
434471
method = "P"
435-
if be.logt
436-
models = [@eval @formula($i ~ $(be.formulation)) for i in be.vars]
437-
else
438-
models = [begin
439-
rfo = @eval @formula(0 ~ $(be.formulation))
440-
lhs = functional_term(log, i)
441-
FormulaTerm(lhs, rfo.rhs)
442-
end for i in be.vars
443-
#@eval @formula(log(Term($i)) ~ $(be.formulation)) for i in be.vars
444-
]
445-
end
472+
473+
models = models_create(be, estimator, design)
446474

447475
elseif design in ("2X2", "2X2X2")
448476

@@ -458,39 +486,7 @@ function estimate(be; estimator = "auto", method = "auto", supresswarn = false,
458486
end
459487
end
460488

461-
if estimator == "glm"
462-
if be.logt
463-
models = [@eval @formula($i ~ $(be.formulation) + $(be.period) + $(be.sequence) + $(be.subject)) for i in be.vars]
464-
else
465-
models = [begin
466-
rfo = @eval @formula(0 ~ $(be.formulation) + $(be.period) + $(be.sequence) + $(be.subject))
467-
lhs = functional_term(log, i)
468-
FormulaTerm(lhs, rfo.rhs)
469-
end for i in be.vars]
470-
end
471-
elseif estimator == "met"
472-
if be.logt
473-
models = [@eval @formula($i ~ $(be.formulation) + $(be.period) + $(be.sequence)) for i in be.vars]
474-
else
475-
models = [begin
476-
rfo = @eval @formula(0 ~ $(be.formulation) + $(be.period) + $(be.sequence))
477-
lhs = functional_term(log, i)
478-
FormulaTerm(lhs, rfo.rhs)
479-
end for i in be.vars]
480-
end
481-
elseif estimator == "mm"
482-
if be.logt
483-
models = [@eval @formula($i ~ $(be.formulation) + $(be.period) + $(be.sequence) + (1| $(be.subject) )) for i in be.vars]
484-
else
485-
models = [begin
486-
rfo = @eval @formula(0 ~ $(be.formulation) + $(be.period) + $(be.sequence) + (1| $(be.subject) ))
487-
lhs = functional_term(log, i)
488-
FormulaTerm(lhs, rfo.rhs)
489-
end for i in be.vars]
490-
end
491-
else
492-
error("Unknown estimator!")
493-
end
489+
models = models_create(be, estimator)
494490

495491
else
496492
if !(method in ("A", "B", "C"))
@@ -513,47 +509,16 @@ function estimate(be; estimator = "auto", method = "auto", supresswarn = false,
513509
estimator = "met"
514510
end
515511

516-
if estimator == "glm"
517-
if be.logt
518-
models = [@eval @formula($i ~ $(be.formulation) + $(be.period) + $(be.sequence) + $(be.subject)) for i in be.vars]
519-
else
520-
models = [begin
521-
rfo = @eval @formula(0 ~ $(be.formulation) + $(be.period) + $(be.sequence) + $(be.subject))
522-
lhs = functional_term(log, i)
523-
FormulaTerm(lhs, rfo.rhs)
524-
end for i in be.vars]
525-
end
526-
elseif estimator == "met"
527-
if be.logt
528-
models = [@eval @formula($i ~ $(be.formulation) + $(be.period) + $(be.sequence)) for i in be.vars]
529-
else
530-
models = [begin
531-
rfo = @eval @formula(0 ~ $(be.formulation) + $(be.period) + $(be.sequence))
532-
lhs = functional_term(log, i)
533-
FormulaTerm(lhs, rfo.rhs)
534-
end for i in be.vars]
535-
end
536-
elseif estimator == "mm"
537-
if be.logt
538-
models = [@eval @formula($i ~ $(be.formulation) + $(be.period) + $(be.sequence) + (1| $(be.subject) )) for i in be.vars]
539-
else
540-
models = [begin
541-
rfo = @eval @formula(0 ~ $(be.formulation) + $(be.period) + $(be.sequence) + (1| $(be.subject) ))
542-
lhs = functional_term(log, i)
543-
FormulaTerm(lhs, rfo.rhs)
544-
end for i in be.vars]
545-
546-
end
547-
else
548-
error("Unknown estimator!")
549-
end
512+
models = models_create(be, estimator)
550513
end
551514

515+
if rsabe != :none && design in ("parallel", "2X2", "2X2X2") && !supresswarn @warn("rsabe option used with unsupported design, nothing will be done!") end
516+
552517
####################################
553518
# ESTIMATION (fitting)
554519
####################################
555520
df = DataFrame(Parameter = String[], Metric = String[], PE = Float64[], SE = Float64[], DF = Float64[], lnLCI = Float64[], lnUCI = Float64[], GMR = Float64[], LCI = Float64[], UCI = Float64[], level = Float64[])
556-
dfdict = Dict(:result => df)
521+
dfdict = Dict{Symbol, DataFrame}(:result => df)
557522
# If GLM used
558523
if estimator == "glm"
559524

@@ -585,9 +550,6 @@ function estimate(be; estimator = "auto", method = "auto", supresswarn = false,
585550
σ²,
586551
cvfromvar(σ²) * 100
587552
))
588-
589-
590-
591553
end
592554

593555
# If Metida Used
@@ -598,14 +560,14 @@ function estimate(be; estimator = "auto", method = "auto", supresswarn = false,
598560
if method == "B"
599561

600562
results = [fit!(LMM(m, be.data;
601-
random = Metida.VarEffect(@eval(Metida.@covstr(1| $(be.subject))), Metida.SI),
563+
random = Metida.VarEffect(functional_term(|, 1, be.subject), Metida.SI),
602564
contrasts = Dict(be.formulation => DummyCoding(base = be.reference)))) for m in models]
603565

604566
elseif method == "C"
605-
567+
606568
results = [fit!(LMM(m, be.data;
607-
random = Metida.VarEffect(@eval(Metida.@covstr($(be.formulation)|$(be.subject))), Metida.CSH),
608-
repeated = Metida.VarEffect(@eval(Metida.@covstr($(be.formulation)|$(be.subject))), Metida.DIAG),
569+
random = Metida.VarEffect(functional_term(|, be.formulation, be.subject), Metida.CSH),
570+
repeated = Metida.VarEffect(functional_term(|, be.formulation, be.subject), Metida.DIAG),
609571
contrasts = Dict(be.formulation => DummyCoding(base = be.reference)))) for m in models]
610572

611573
else
@@ -686,6 +648,60 @@ function estimate(be; estimator = "auto", method = "auto", supresswarn = false,
686648

687649
end
688650

651+
if rsabe != :none && !(design in ("parallel", "2X2", "2X2X2"))
652+
if !(rsabe in (:ema, :fda, :eeu)) rsabe = :unknown end
653+
cvdf = dfdict[:scaled_ci] = DataFrame(Metric = String[], Reference = String[], σ²= Float64[], CV = Float64[], Method = String[], LCIL = Float64[], UCIL = Float64[] )
654+
655+
if be.logt
656+
cv_models = [FormulaTerm(term(i), term(be.period) + term(be.sequence) + term(be.subject)) for i in be.vars]
657+
else
658+
cv_models = [begin
659+
rfo = FormulaTerm(0, term(be.period) + term(be.sequence) + term(be.subject))
660+
lhs = functional_term(log, i)
661+
FormulaTerm(lhs, rfo.rhs)
662+
end for i in be.vars]
663+
end
664+
cv_data = filter(be.formulation => x -> x == be.reference, be.data)
665+
666+
cv_results = [fit(LinearModel, m, cv_data; dropcollinear = true) for m in cv_models]
667+
668+
for cvr in cv_results
669+
σ² = GLM.dispersion(cvr.model, true)
670+
s_wr = sqrt(σ²)
671+
cv_w = cvfromvar(σ²) * 100
672+
673+
if rsabe in (:ema, :eeu)
674+
lcil = cv_w > 30 ? max(69.84, round(exp(-0.760 * s_wr) * 100, digits = 2)) : 80.0
675+
ucil = cv_w > 30 ? min(143.19, round(exp(0.760 * s_wr) * 100, digits = 2)) : 125.00
676+
elseif rsabe == :fda
677+
if ntid
678+
lcil = s_wr > 0.1 ? round(exp(-0.893 * s_wr) * 100, digits = 2) : 90.0
679+
ucil = s_wr > 0.1 ? round(exp(0.893 * s_wr) * 100, digits = 2) : 111.11
680+
else
681+
lcil = cv_w >= 30 ? round(exp(-0.893 * s_wr) * 100, digits = 2) : 80.0
682+
ucil = cv_w >= 30 ? round(exp(0.893 * s_wr) * 100, digits = 2) : 125.00
683+
end
684+
else
685+
rsabe = :unknown
686+
lcil = NaN
687+
ucil = NaN
688+
end
689+
690+
push!(cvdf, (
691+
coefnames(cvr.mf.f.lhs),
692+
string(be.reference),
693+
σ²,
694+
cv_w,
695+
string(rsabe),
696+
lcil,
697+
ucil,
698+
))
699+
end
700+
701+
end
702+
703+
704+
689705
BEResults(be, results, dfdict, estimator, method)
690706
end
691707

@@ -696,4 +712,13 @@ Returns dataframe with bioequivalence results.
696712
"""
697713
function result(beres::BEResults)
698714
beres.df[:result]
715+
end
716+
717+
"""
718+
variance(beres::BEResults)
719+
720+
Returns dataframe with variance.
721+
"""
722+
function variance(beres::BEResults)
723+
beres.df[:var]
699724
end

‎test/tests.jl‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -57,6 +57,8 @@ rds1missing[36, :logVar] = missing
5757
@test berAglm.method == "A"
5858
@test berAglm.estimator== "glm"
5959

60+
berAglmrsabe = MetidaBioeq.estimate(be; estimator = "met", method = "C", rsabe = :ema)
61+
6062
@test_nowarn MetidaBioeq.result(berCmet)
6163
# Inappropriate
6264
# Try "glm" && "B" "C"

0 commit comments

Comments
 (0)