lmmformula

change.log

@@ -1,4 +1,5 @@
-v0.14
+uv0.14
+  * lmmformula
   * Unstructured Covariance (UN)
   * initial functions for analytical gradient vector
   * docs

docs/src/api.md

@@ -16,6 +16,11 @@
 
 [`fit!`](@ref)
 
+### @lmmformula
+```@docs
+Metida.@lmmformula
+```
+
 ### Metida.CovarianceType
 ```@docs
 Metida.CovarianceType

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.

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")

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)

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

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