Merge pull request #7 from PharmCat/dev

Dev

.github/workflows/Documenter.yml

@@ -23,7 +23,7 @@ jobs:
     runs-on: ubuntu-latest
     timeout-minutes: 30
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v3
       - uses: julia-actions/setup-julia@v1
         with:
           version: "1.8"

Project.toml

@@ -3,7 +3,7 @@ uuid = "00e1d80e-504f-4d7b-946f-0aea2e7abfe3"
 keywords = ["bioequivalence"]
 desc = "Bioequivalence"
 authors = ["Vladimir Arnautov <mail@pharmcat.net>"]
-version = "0.2.2"
+version = "0.2.3"
 
 [deps]
 CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597"
@@ -19,10 +19,10 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 CategoricalArrays = "0.10"
 DataFrames = "1"
 Distributions = "0.21, 0.22, 0.23, 0.24, 0.25"
-GLM = "1.8.2"
-Metida = "0.14.6, 0.14.7"
-MetidaNCA = "0.5.7"
-MixedModels = "4.11, 4.12, 4.13, 4.14"
+GLM = "1.8"
+Metida = "0.15"
+MetidaNCA = "0.5"
+MixedModels = "4"
 StatsBase = "0.30, 0.31, 0.32, 0.33, 0.34"
 julia = "1"
 

docs/src/api.md

@@ -6,12 +6,17 @@
 MetidaBioeq.bioequivalence
 ```
 
-### MetidaBioeq.makeseq
+### MetidaBioeq.estimate
 ```@docs
-MetidaBioeq.makeseq
+MetidaBioeq.estimate
 ```
 
 ### MetidaBioeq.result
 ```@docs
 MetidaBioeq.result
+```
+
+### MetidaBioeq.makeseq
+```@docs
+MetidaBioeq.makeseq
 ```

docs/src/example.md

@@ -4,8 +4,8 @@
 using MetidaBioeq, CSV, DataFrames, CategoricalArrays;
 
 bedf2x2 = CSV.File(joinpath(dirname(pathof(MetidaBioeq)), "..", "test", "csv",  "2x2rds1.csv")) |> DataFrame
-transform!(bedf2x2, :Subject => categorical, renamecols = false)
-transform!(bedf2x2, :Period => categorical, renamecols = false)
+transform!(bedf2x2, :Subj => categorical, renamecols = false)
+transform!(bedf2x2, :Per => categorical, renamecols = false)
 bedf2x2.logVar = log.(bedf2x2.Var)
 
 bedf2x2x4 = CSV.File(joinpath(dirname(pathof(MetidaBioeq)), "..", "test", "csv",  "2x2x4rds1.csv")) |> DataFrame
@@ -28,10 +28,10 @@ Simple 2x2 study. Data (`var`) not log-transgormed.
 ```@example beexample
     be1 = MetidaBioeq.bioequivalence(bedf2x2, 
     vars = :Var, 
-    subject = :Subject, 
-    formulation = :Formulation, 
-    period = :Period,
-    sequence = :Sequence, 
+    subject = :Subj, 
+    formulation = :Trt, 
+    period = :Per,
+    sequence = :Seq, 
     reference = "R",
     design = "2x2",
     autoseq = true,

docs/src/index.md

@@ -13,7 +13,103 @@ The United States Food and Drug Administration (FDA) has defined bioequivalence
 * Julia ≥1.8+
 * MixedModels ≥4.11
 * GLM ≥1.8.2
-* Metida ≥0.14.6
+* Metida ≥0.15
+
+### Models 
+
+Basic models for bioequivalence.
+
+#### A
+
+Standard fixed-effect model. For parallel design Welch's correction not used.
+
+```julia
+# Parallel design GLM
+fit(LinearModel, @formula(var ~ formulation), data; 
+contrasts = Dict(formulation => DummyCoding(base = reference)),
+)
+
+
+# Cross-over designs GLM
+fit(LinearModel, @formula(var ~ formulation + period + sequence + subject), data; 
+contrasts = Dict(formulation => DummyCoding(base = reference)),
+)
+```
+
+Validated with:
+
+* Schütz, H., Labes, D., & Fuglsang, A. (2014). Reference datasets for 2-treatment, 2-sequence, 2-period bioequivalence studies. The AAPS journal, 16(6), 1292–1297. https://doi.org/10.1208/s12248-014-9661-0
+* Fuglsang, A., Schütz, H., & Labes, D. (2015). Reference datasets for bioequivalence trials in a two-group parallel design. The AAPS journal, 17(2), 400–404. https://doi.org/10.1208/s12248-014-9704-6
+
+
+#### Type B
+
+Random-effect model, no covariance.
+
+```julia
+# MixedModels
+
+fit(MixedModel, @formula(var ~ formulation + period + sequence + (1|subject)), data; 
+    contrasts = Dict(formulation => DummyCoding(base = reference)),
+    REML=true
+)
+
+# Metida
+lmm = LMM(@formula(var~sequence+period+formulation), data;
+random = VarEffect(@covstr(1|subject), SI),
+contrasts = Dict(formulation => DummyCoding(base = reference)),
+)
+
+```
+
+Validated with: 
+
+* Schütz, H., Labes, D., Tomashevskiy, M., la Parra, M. G., Shitova, A., & Fuglsang, A. (2020). Reference Datasets for Studies in a Replicate Design Intended for Average Bioequivalence with Expanding Limits. The AAPS journal, 22(2), 44. https://doi.org/10.1208/s12248-020-0427-6
+
+
+#### Type C
+
+Random-effect model, with covariance (FDA reference code).
+
+```julia
+# Metida
+lmm =LMM(@formula(var~sequence+period+formulation), data;
+random = VarEffect(@covstr(formulation|subject), CSH),
+repeated = VarEffect(@covstr(formulation|subject), DIAG),
+contrasts = Dict(formulation => DummyCoding(base = reference)),
+)
+```
+
+Validated against SPSS model:
+
+```
+MIXED var BY period formulation sequence subject
+  /CRITERIA=CIN(90) MXITER(200) MXSTEP(20) SCORING(2) SINGULAR(0.000000000001) HCONVERGE(0,
+    RELATIVE) LCONVERGE(0.0000000000001, RELATIVE) PCONVERGE(0, RELATIVE)
+  /FIXED=period formulation sequence | SSTYPE(3)
+  /METHOD=REML
+  /RANDOM=formulation | SUBJECT(subject) COVTYPE(CSH)
+  /REPEATED=formulation | SUBJECT(subject*period) COVTYPE(DIAG)
+  /EMMEANS=TABLES(formulation) COMPARE REFCAT(FIRST) ADJ(LSD).
+```
+
+What complies with the FDA code ([FDA Guidance for Industry: Statistical Approaches to Establishing Bioequivalence](https://www.fda.gov/media/70958/download), APPENDIX F):
+
+```
+PROC MIXED;
+CLASSES SEQ SUBJ PER TRT;
+MODEL Y = SEQ PER TRT/ DDFM=SATTERTH;
+RANDOM TRT/TYPE=FA0(2) SUB=SUBJ G;
+REPEATED/GRP=TRT SUB=SUBJ;
+ESTIMATE 'T vs. R' TRT 1 -1/CL ALPHA=0.1;
+```
+
+**In the Random statement, TYPE=FA0(2) could possibly be replaced by TYPE=CSH. This
+guidance recommends that TYPE=UN not be used, as it could result in an invalid (i.e., not non-
+negative definite) estimated covariance matrix**
+
+[Exported SPSS results](https://github.com/PharmCat/Metida.jl/blob/f76cd21fe7aceb01e5d0fcf75137c9ff6f3a8377/validation/RDS-OUTPUT.docx).
+
 
 ### Reference:
 

src/bioequivalence.jl

@@ -20,19 +20,23 @@ function functional_term(f, arg_expr...)
 end
 
 """
-    bioquivalence(data;
-    vars = nothing,
-    subject = :subject,
-    period = :period,
-    formulation = :formulation,
-    sequence = :sequence,
-    stage = nothing,
-    reference = nothing,
-    design = nothing,
-    io::IO = stdout,
-    seqcheck = true,
-    dropcheck = true,
-    logt = true)
+    bioequivalence(data;
+        vars = nothing,
+        subject::Union{String, Symbol},
+        period::Union{String, Symbol, Nothing} = nothing,
+        formulation::Union{String, Symbol},
+        sequence::Union{String, Symbol, Nothing} = nothing,
+        stage::Union{String, Symbol, Nothing} = nothing,
+        reference::Union{String, Symbol, Nothing} = nothing,
+        design::Union{String, Symbol, Nothing} = nothing,
+        io::IO = stdout,
+        seqcheck::Bool = true,
+        designcheck::Bool = true,
+        dropcheck::Bool = true,
+        info::Bool = true,
+        warns::Bool = true,
+        autoseq::Bool = false,
+        logt::Bool = true)
 
 * `vars` - variabel's column(s);
 * `subject` - subject's column;
@@ -42,8 +46,12 @@ end
 * `stage` - stage's column;
 * `reference` - reference value for `formulation` column;
 * `design` - design: "parallel", "2X2", "2X2X2", "2X2X4", ets. (formulations X sequences X periods);
-* `seqcheck` - check sequencs;     
+* `seqcheck` - check sequencs; 
+* `designcheck` - check design correctness;  
 * `dropcheck` - dropuot check;
+* `info` - show information;
+* `warns` - show warnings;
+* `autoseq` - try to make sequence collumn;
 * `logt` - if `true` (default) data is already log-transformed, else `log()` will be used.
 """
 function bioequivalence(data;
@@ -57,6 +65,7 @@ function bioequivalence(data;
     design::Union{String, Symbol, Nothing} = nothing,
     io::IO = stdout,
     seqcheck::Bool = true,
+    designcheck::Bool = true,
     dropcheck::Bool = true,
     info::Bool = true,
     warns::Bool = true,
@@ -189,7 +198,7 @@ function bioequivalence(data;
         if seqcheck && !isnothing(sequence)
             for i = 1:obsnum
                 if getcol(data, sequence)[i] != subjdict[getcol(data, subject)[i]]
-                    error("Sequence error or data is incomplete! \n Subject: $(getcol(data, subject)[i]), Sequence: $(getcol(data, sequence)[i]), auto: $(subjdict[getcol(data, subject)[i]])")
+                    error("Sequence error or data is incomplete! \n Subject: $(getcol(data, subject)[i]), Sequence: $(getcol(data, sequence)[i]), auto: $(subjdict[getcol(data, subject)[i]]), use `seqcheck = false` keyword to disable sequence check.")
                 end
             end
             if length(unique(length.(sequences))) > 1
@@ -202,13 +211,13 @@ function bioequivalence(data;
             info && @info "Trying to find out the design..."
             design = "$(length(formulations))X$(length(sequences))X$(length(periods))"
             info && @info  "Seems design type is: $design"
-        else
+        elseif designcheck
             if design == "2X2" design = "2X2X2" end
             spldes = split(design, "X")
             if length(spldes) != 3 error("Unknown design type. Use fXsXp format or \"2Х2\".") end
             if length(formulations) != parse(Int, spldes[1]) error("Design error: formulations count wrong!") end
             if length(sequences) != parse(Int, spldes[2]) error("Design error: sequences count wrong!") end
-            if length(periods) != parse(Int, spldes[3]) error("Design error: periods count wrong!") end
+            if length(periods) != parse(Int, spldes[3]) error("Design error: periods count wrong! length(periods) = $(length(periods)), desigh = $design , use `designcheck = false` keyword to disable design check.") end
             info && @info "Design type seems fine..."
         end
     else
@@ -366,8 +375,9 @@ function estimate(be; estimator = "auto", method = "auto", supresswarn = false,
 
         end
     end
-
+    ##############################
     # MODEL SELECTION
+    ##############################
     if design == "parallel"
 
         if estimator != "glm" && !supresswarn @warn("Design is parallel, but estimator not GLM, GLM will be used!") end 
@@ -377,7 +387,13 @@ function estimate(be; estimator = "auto", method = "auto", supresswarn = false,
         if be.logt
             models = [@eval @formula($i ~ $(be.formulation)) for i in be.vars]
         else
-            models = [@eval @formula(log(Term($i)) ~ $(be.formulation)) for i in be.vars]
+            models = [begin
+                rfo =  @eval  @formula(0 ~ $(be.formulation))
+                lhs = functional_term(log, i)
+                FormulaTerm(lhs, rfo.rhs) 
+            end for i in be.vars
+            #@eval @formula(log(Term($i)) ~ $(be.formulation)) for i in be.vars
+            ]
         end
 
     elseif design in ("2X2", "2X2X2") 
@@ -457,7 +473,7 @@ function estimate(be; estimator = "auto", method = "auto", supresswarn = false,
                 rfo = @eval @formula(0 ~ $(be.formulation) + $(be.period) + $(be.sequence) + $(be.subject)) 
                 lhs = functional_term(log, i)
                 FormulaTerm(lhs, rfo.rhs)
-            end for i in be.vars]
+                end for i in be.vars]
             end
         elseif estimator == "met"
             if be.logt
@@ -476,7 +492,6 @@ function estimate(be; estimator = "auto", method = "auto", supresswarn = false,
                 models = [begin 
                 rfo = @eval @formula(0 ~ $(be.formulation) + $(be.period) + $(be.sequence) + (1| $(be.subject) ))
                 lhs = functional_term(log, i)
-                #lhs = term(i)
                 FormulaTerm(lhs, rfo.rhs) 
                 end for i in be.vars]
 
@@ -486,6 +501,9 @@ function estimate(be; estimator = "auto", method = "auto", supresswarn = false,
         end
     end
 
+    ####################################
+    # ESTIMATION (fitting)
+    ####################################
     # If GLM used 
     if estimator == "glm"
 
@@ -539,7 +557,7 @@ function estimate(be; estimator = "auto", method = "auto", supresswarn = false,
             lnLCI = PE - SE * quantile(dist, 1-alpha)
             lnUCI = PE + SE * quantile(dist, 1-alpha)
             push!(df, (string(coefnames(i)[2], " - ", be.reference),
-                coefnames(i.mf.f.lhs),
+                responsename(i),
                 PE,
                 SE,
                 DF,
@@ -583,9 +601,8 @@ function estimate(be; estimator = "auto", method = "auto", supresswarn = false,
                 ))
         end
 
-
     end
-    BEResults(results, df, estimator, method)
+    BEResults(results, Dict(:result => df), estimator, method)
 end
 
 """
@@ -594,5 +611,5 @@ end
 Returns dataframe with bioequivalence results.
 """
 function result(beres::BEResults)
-    beres.df
+    beres.df[:result]
 end

src/types.jl

@@ -18,7 +18,7 @@ end
 
 struct BEResults
     models::Vector
-    df::DataFrame
+    df::Dict
     estimator::String
     method::String
 end
@@ -48,6 +48,7 @@ function Base.show(io::IO, obj::Bioequivalence)
 end
 
 function Base.show(io::IO, obj::BEResults)
-    println(io, obj.df)
+    println(io, "  Result:")
+    println(io, obj.df[:result])
     print(io, "Estimator: $(obj.estimator) ($(obj.method))")
 end

test/csv/2x2rds1.csv

@@ -1,4 +1,4 @@
-Subject	Sequence	Period	Formulation	Var
+Subj	Seq	Per	Trt	Var
 3	TR	1	T	225.95
 1	RT	1	R	181.09
 2	RT	1	R	114.48