diff --git a/change.log b/change.log index 4d410f66..110588a4 100644 --- a/change.log +++ b/change.log @@ -1,4 +1,5 @@ -v0.14 +uv0.14 + * lmmformula * Unstructured Covariance (UN) * initial functions for analytical gradient vector * docs diff --git a/docs/src/api.md b/docs/src/api.md index 1293b779..1758e080 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -16,6 +16,11 @@ [`fit!`](@ref) +### @lmmformula +```@docs +Metida.@lmmformula +``` + ### Metida.CovarianceType ```@docs Metida.CovarianceType diff --git a/docs/src/instanduse.md b/docs/src/instanduse.md index ec6e6b3a..e8b03e2c 100644 --- a/docs/src/instanduse.md +++ b/docs/src/instanduse.md @@ -60,7 +60,7 @@ nothing # hide Make model with `@formula` macro from `StatsModels`. Define `random` and `repreated` effects with [`Metida.VarEffect`](@ref) using [`Metida.@covstr`](@ref) macros. Left side of `@covstr` is model of effect and -right side is a effect itself. [`Metida.HeterogeneousCompoundSymmetry`](@ref) and [`Metida.Diagonal`](@ref) in example bellow is a model of variance-covariance structure. +right side is a effect itself. [`Metida.HeterogeneousCompoundSymmetry`](@ref) and [`Metida.Diagonal`](@ref) in example bellow is a model of variance-covariance structure. See also [`@lmmformula`](@ref) macro. !!! note In some cases levels of repeated effect should not be equal inside each level of subject or model will not have any sense. For example, it is assumed that usually CSH or UN (Unstructured) using with levels of repeated effect is different inside each level of subject. diff --git a/src/Metida.jl b/src/Metida.jl index fd5f5920..630c2025 100644 --- a/src/Metida.jl +++ b/src/Metida.jl @@ -63,6 +63,7 @@ include("linearalgebra.jl") include("options.jl") include("modelresult.jl") include("lmmdata.jl") +include("lmmformula.jl") include("lmm.jl") include("reml.jl") include("ml.jl") diff --git a/src/lmm.jl b/src/lmm.jl index af4981d1..70982895 100644 --- a/src/lmm.jl +++ b/src/lmm.jl @@ -20,6 +20,8 @@ Make Linear-Mixed Model object. `random`: vector of random effects or single random effect `repeated`: is a repeated effect (only one) + +See also: [`@lmmformula`](@ref) """ struct LMM{T<:AbstractFloat} <: MetidaModel model::FormulaTerm @@ -75,6 +77,7 @@ struct LMM{T<:AbstractFloat} <: MetidaModel if repeated === nothing repeated = NOREPEAT end + if random === nothing random = VarEffect(Metida.@covstr(0|1), Metida.RZero()) end @@ -96,6 +99,10 @@ struct LMM{T<:AbstractFloat} <: MetidaModel end end +function LMM(f::LMMformula, data; contrasts=Dict{Symbol,Any}()) + LMM(f.formula, data; contrasts=contrasts, random = f.random, repeated = f.repeated) +end + ################################################################################ """ thetalength(lmm::LMM) diff --git a/src/lmmformula.jl b/src/lmmformula.jl new file mode 100644 index 00000000..9d80beb7 --- /dev/null +++ b/src/lmmformula.jl @@ -0,0 +1,123 @@ +struct LMMformula + formula + random + repeated +end + +""" + @lmmformula(formula, args...) + +Macro for made formula with variance-covariance structure representation. `@lmmformula` +could be used for shorter `LMM` construction. + +Example: + +``` +lmm = Metida.LMM(@lmmformula(var~sequence+period+formulation, +random = formulation|subject:CSH, +repeated = formulation|subject:DIAG), +df0) +``` + +equal to: + +``` +lmm = LMM(@formula(var~sequence+period+formulation), df0; +random = Metida.VarEffect(@covstr(formulation|subject), CSH), +repeated = Metida.VarEffect(@covstr(formulation|subject), DIAG), +) +``` + +`@lmmformula` have 3 components - 1'st is a formula for fixed effect, it defined +like in `StstsModels` (1st argument just provided to `@formula` macro). Other arguments +should be defined like keywords. `repeated` keyword define repeated effect part, +`random` define random effect part. You can use several random factors as in example bellow: + +``` +lmm = Metida.LMM(Metida.@lmmformula(var~sequence+period+formulation, +random = formulation|subject:CSH, +random = 1|subject:DIAG, +repeated = formulation|subject:DIAG), +df0) +``` + +`random` or `repeated` structure made by template: + +`effect formula` | `blocking factor` [/ `nested factor`] [: `covariance structure`] + +`|` - devide effect formula form blocking factor definition (necessarily), +`/` and `:` modificator are optional. `/` work like in MixedModels or in RegressionFormulae - +expand factor `f|a/b` to `f|a` + `f|a&b`. It can't be used in repeated effect defenition. + +`:` - covariance structure defined right after `:` (SI, DIAG, CS, CSH, ets...), +if `:` not used then SI used for this effect. + +Terms like `a+b` or `a*b` shuould not be used as a blocking factors. + +""" +macro lmmformula(formula, args...) + f = eval(:(@formula($formula))) + ranfac = Metida.VarEffect[] + repeff = nothing + for ex in args + if ex.head != :(=) error("Error = ") end + a, e = ex.args + ear = e.args + # Check if structure defined, else SI used + if ear[1] == :(:) + ct = eval(ear[3]) + e = ear[2] + else + ct = SI + end + ear = e.args + if a == :random + if ear[1] == :(|) + # For simple random factor + if isa(ear[3], Symbol) + mcs = eval(:(@covstr($e))) + push!(ranfac, VarEffect(mcs, ct)) + # For nested factors defined with "/" + elseif length(ear[3].args) == 3 && ear[3].args[1] == :(/) + o = ear[3].args[2] + mcs = eval(:(@covstr($(ear[2]) | $o))) + push!(ranfac, VarEffect(mcs, ct)) + u = ear[3].args[3] + t = Expr(:call, :&, o, u) + mcs = eval(:(@covstr($(ear[2]) | $t))) + push!(ranfac, VarEffect(mcs, ct)) + # Other cases + else + mcs = eval(:(@covstr($e))) + push!(ranfac, VarEffect(mcs, ct)) + end + else + error("Random term should include `|` operator in syntax") + end + + elseif a == :repeated + if ear[1] == :(|) + if isa(ear[3], Symbol) + + mcs = eval(:(@covstr($e))) + repeff = VarEffect(mcs, ct) + + elseif length(ear[3].args) == 3 && ear[3].args[1] == :(/) + + error("Repeated effect can't include `/`") + + else + mcs = eval(:(@covstr($e))) + repeff = VarEffect(mcs, ct) + end + else + error("Random term should include `|` operator in syntax") + end + + else + error("Unknown type") + end + end + if length(ranfac) == 0 ranfac = nothing end + return LMMformula(f, ranfac, repeff) +end diff --git a/test/test.jl b/test/test.jl index dc131ba7..629dbb86 100644 --- a/test/test.jl +++ b/test/test.jl @@ -236,7 +236,7 @@ end @test std[1] ≈ 0.3345433916523553 atol=1E-8 @test cn[1] ≈ 1.577492862311838 atol=1E-8 end -@testset " Model: CSH/DIAG (rholinkf = :psigm) " begin +@testset " Model: CSH/DIAG (rholinkf = :psigm) & lmmformula " begin lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0; random = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.CSH), repeated = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.DIAG), @@ -248,6 +248,15 @@ end @test std[1] ≈ 0.3345433910321999 atol=1E-8 @test cn[1] ≈ 1.5774928621922844 atol=1E-8 + std = stderror(lmm) + + lmm = Metida.LMM(Metida.@lmmformula(var~sequence+period+formulation, + random = formulation|subject:Metida.CSH, + repeated = formulation|subject:Metida.DIAG), + df0) + Metida.fit!(lmm; rholinkf = :psigm) + @test Metida.m2logreml(lmm) ≈ 10.065239006121315 atol=1E-6 + end ################################################################################ # ftdf / 1fptime.csv