Merge pull request #51 from PharmCat/dev

Dev

Project.toml

@@ -3,7 +3,7 @@ uuid = "a1dec852-9fe5-11e9-361f-8d9fde67cfa2"
 keywords = ["lenearmodel", "mixedmodel"]
 desc = "Mixed-effects models with flexible covariance structure."
 authors = ["Vladimir Arnautov <[email protected]>"]
-version = "0.16.2"
+version = "0.16.3"
 
 [deps]
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"

docs/src/api.md

@@ -406,6 +406,13 @@ Metida.SpatialPowerD
 Metida.ScaledWeightedCov
 ```
 
+### Metida.ACOV
+
+```@docs
+Metida.ACOV
+```
+
+
 ### Metida.dof_contain
 
 ```@docs

docs/src/examples.md

@@ -1,10 +1,15 @@
-### Example 1 - Continuous and categorical predictors
+## Example 1 - Continuous and categorical predictors
 
 ```@example lmmexample
 using Metida, CSV, DataFrames, CategoricalArrays, MixedModels;
 
 rds = CSV.File(joinpath(dirname(pathof(Metida)), "..", "test", "csv",  "1fptime.csv"); types = [String, String, Float64, Float64]) |> DataFrame
 
+rds2 = CSV.File(joinpath(dirname(pathof(Metida)), "..", "test", "csv",  "ftdf3.csv"); types = [String,  Float64, Float64, String, String, String, String, String, Float64]) |> DataFrame
+
+
+devday = CSV.File(joinpath(dirname(pathof(Metida)), "..", "test", "csv",  "devday.csv"); types = [Float64, String, String, String ]) |> DataFrame
+
 nothing; # hide
 ```
 
@@ -27,7 +32,7 @@ mm = fit(MixedModel, fm, rds, REML=true)
 println(mm) #hide
 ```
 
-### Example 2 - Two random factors (Penicillin data)
+## Example 2 - Two random factors (Penicillin data)
 
 Metida:
 
@@ -51,7 +56,7 @@ mm = fit(MixedModel, fm2, df, REML=true)
 println(mm) #hide
 ```
 
-### Example 3 - Repeated ARMA/AR/ARH
+## Example 3 - Repeated ARMA/AR/ARH
 
 ```@example lmmexample
 rds = CSV.File(joinpath(dirname(pathof(Metida)), "..", "test", "csv",  "1freparma.csv"); types = [String, String, Float64, Float64]) |> DataFrame
@@ -91,9 +96,9 @@ repeated = VarEffect(@covstr(1|subject&factor), ARH),
 fit!(lmm)
 ```
 
-### Example 4 - SAS relation
+## Example 4 - SAS relation
 
-#### Model 1
+### Model 1
 
 ```
 df0 = CSV.File(joinpath(dirname(pathof(Metida)), "..", "test", "csv", "df0.csv")) |> DataFrame
@@ -116,7 +121,7 @@ REPEATED/GRP=formulation SUB=subject R;
 RUN;
 ```
 
-#### Model 2
+### Model 2
 
 ```
 lmm = LMM(
@@ -138,7 +143,7 @@ REPEATED/GRP=formulation SUB=subject R;
 RUN;
 ```
 
-#### Model 3
+### Model 3
 
 ```
 lmm = LMM(@formula(var ~ sequence + period + formulation), df0;
@@ -157,7 +162,7 @@ RANDOM  subject/TYPE=VC G V;
 RUN;
 ```
 
-### Example 5 - Working with Effects.jl
+## Example 5 - Working with Effects.jl
 
 ```
 using Effects, StatsModels
@@ -172,4 +177,73 @@ table_model = StatsModels.TableRegressionModel(lmm, lmm.mf, lmm.mm)
 emmeans(tm)
 
 effects(Dict(:period => ["1", "2", "3", "4"]), tm)
+```
+
+
+## Unstructured covariance 
+
+Unstructured covariance example.
+
+
+Metida result:
+
+```@example lmmexample
+lmm = Metida.LMM(@formula(response~factor), rds2;
+    random = Metida.VarEffect(Metida.@covstr(r1|subject), UN),
+    )
+Metida.fit!(lmm)
+```
+
+MixedModels result:
+
+```@example lmmexample
+mm = fit(MixedModel, @formula(response ~ factor+ (0+r1|subject)), rds2, REML = true)
+println(mm) #hide
+```
+
+## Aumented covariance (Experimental)
+
+Covariance modificator `ACOV()` can be used as second repeated effect. In this case covariance calculated with existed matrix, 
+that was build at previous step. For example addition `ACOV(AR)` to `DIAG` structure is the same as `ARH` if same blocking factor used.
+
+```@example lmmexample
+   lmm1 = Metida.LMM(@formula(response ~ 1), rds2;
+    repeated = [Metida.VarEffect(Metida.@covstr(r1|subject), Metida.DIAG), Metida.VarEffect(Metida.@covstr(1|subject), Metida.ACOV(Metida.AR))]
+    )
+    Metida.fit!(lmm1)
+```
+
+```@example lmmexample
+    lmm2 = Metida.LMM(@formula(response ~ 1), rds2;
+    repeated = [Metida.VarEffect(Metida.@covstr(r1|subject), Metida.ARH)]
+    )
+    Metida.fit!(lmm2)
+```
+
+R-part of variance-covariance matrix:
+
+```@example lmmexample
+Metida.rmatrix(lmm1, 1)
+```
+
+If nested blocking factor used - covariance modification applyed only within that blocks:
+
+```@example lmmexample
+   lmm = Metida.LMM(@formula(response ~ 1), rds2;
+    repeated = [Metida.VarEffect(Metida.@covstr(r1|subject), Metida.DIAG), Metida.VarEffect(Metida.@covstr(1|subject), Metida.ACOV(Metida.AR))]
+    )
+Metida.fit!(lmm)
+Metida.rmatrix(lmm, 1)
+```
+
+For nested models covariance structure can be expanded as follows:
+* the first layer describes unstructured the device-device covariance;
+* the second layer adds the time covariance for each device
+
+```@example lmmexample
+lmm = Metida.LMM(@formula(resp ~ 0 + device), devday;
+    repeated = [Metida.VarEffect(Metida.@covstr(device|subj&day), Metida.UN), 
+    Metida.VarEffect(Metida.@covstr(1|subj&device), Metida.ACOV(Metida.AR))]
+    )
+    Metida.fit!(lmm)
 ```

src/Metida.jl

@@ -36,6 +36,7 @@ SPPOW, SpatialPower,
 SPGAU, SpatialGaussian, RawCoding,
 UN, Unstructured,
 CovarianceType,
+ACOV,
 fit, fit!, LMM, VarEffect, theta, logreml, m2logreml, m2logml, logml, thetalength, dof_satter, dof_contain, rankx, caic, lcontrast, typeiii, estimate, contrast,
 gmatrix, rmatrix, vmatrix!, responsename, nblocks, raneff,
 AbstractCovarianceType, AbstractCovmatMethod, MetidaModel,

src/rmat.jl

@@ -342,6 +342,76 @@ function rmat!(mx, θ, rz::AbstractMatrix, ::UN_, ::Int)
     mulαβαtinc!(mx, rz, rcov)
     return mx
 end
+
+# ACOV (AR)
+function rmat!(mx::AbstractMatrix{T}, θ, ::AbstractMatrix, ct::ACOV_{AR_}, ::Int) where T
+    s  = size(mx, 1)
+    ρ  = θ[1]
+    if s > 1
+        for n = 2:s
+            mxnn = mx[n, n]
+            @inbounds @simd for m = 1:n-1
+                mxmm = mx[m, m]
+                mxmn = mx[m, n]
+                cov  = sqrt(mxnn * mxmm) * ρ ^ (n - m)
+                if ct.a != 0
+                    if ct.a == 1
+                        mxmn = zero(T)
+                    elseif ct.a == 2
+                        cov = zero(T)
+                    elseif ct.a == 3
+                        @warn "covariance have existed value!"
+                    elseif ct.a == 4
+                        @warn "covariance have existed value! replaced..."
+                        mxmn = zero(T)
+                    elseif ct.a == 5
+                        @warn "covariance have existed value! Save existed..."
+                        cov = zero(T)
+                    else
+                        error("covariance have existed value!")
+                    end
+                end
+                mx[m, n] = mxmn + cov
+            end
+        end
+    end
+    return mx
+end
+
+# ACOV (CS)
+function rmat!(mx::AbstractMatrix{T}, θ, ::AbstractMatrix, ct::ACOV_{CS_}, ::Int) where T
+    s  = size(mx, 1)
+    ρ  = θ[1]
+    if s > 1
+        for n = 2:s
+            mxnn = mx[n, n]
+            @inbounds @simd for m = 1:n-1
+                mxmm = mx[m, m]
+                mxmn = mx[m, n]
+                cov  = sqrt(mxnn * mxmm) * ρ 
+                if ct.a != 0
+                    if ct.a == 1
+                        mxmn = zero(T)
+                    elseif ct.a == 2
+                        cov = zero(T)
+                    elseif ct.a == 3
+                        @warn "covariance have existed value!"
+                    elseif ct.a == 4
+                        @warn "covariance have existed value! replaced..."
+                        mxmn = zero(T)
+                    elseif ct.a == 5
+                        @warn "covariance have existed value! Save existed..."
+                        cov = zero(T)
+                    else
+                        error("covariance have existed value!")
+                    end
+                end
+                mx[m, n] = mxmn + cov
+            end
+        end
+    end
+    return mx
+end
 ###############################################################################
 ###############################################################################
 ###############################################################################

src/varstruct.jl

@@ -31,7 +31,7 @@ Macros for random/repeated effect model.
 # Example
 
 ```julia
-@covstr(factor|subject)
+@covstr(model|subject)
 ```
 """
 macro covstr(ex)
@@ -64,7 +64,7 @@ Random/repeated effect.
 
 !!! note
 
-    Categorical factors are coded with `FullDummyCoding()` by default, use `coding` for other contrast codeing.
+    Categorical factors are coded with `FullDummyCoding()` by default, use `coding` for other contrast coding.
 
 # Example
 
@@ -352,6 +352,9 @@ struct CovStructure{T, T2} <: AbstractCovarianceStructure
                         subjblockdict = sabjcrossdicts(subjblockdict, dicts[i])
                     end
                 end
+                if isa(repeated[1].covtype.s, ACOV_)
+                    @warn "Using ACOV covariance additional effect at first position is meaningless."
+                end
             else
                 subjblockdict = nothing
             end
@@ -432,7 +435,15 @@ end
 #                            CONTRAST CODING
 ################################################################################
 
-function fill_coding_dict!(t::T, d::Dict, data) where T <: Union{ConstantTerm, InterceptTerm, FunctionTerm}
+function fill_coding_dict!(t::T, d::Dict, data) where T <: Union{ConstantTerm, InterceptTerm}
+    return d
+end
+function fill_coding_dict!(t::T, d::Dict, data) where T <: FunctionTerm
+    if t.f === +
+        for i in t.args
+            fill_coding_dict!(i, d, data)
+        end
+    end
     return d
 end
 function fill_coding_dict!(t::T, d::Dict, data) where T <: Term

src/vartypes.jl

@@ -32,6 +32,11 @@ struct SPGAUD_ <: AbstractCovarianceType end
 struct UN_ <: AbstractCovarianceType end
 struct ZERO <: AbstractCovarianceType end
 
+struct ACOV_{C <: AbstractCovarianceType} <: AbstractCovarianceType
+    c::C
+    a::Int
+end
+
 ################################################################################
 #                          COVARIANCE TYPE
 ################################################################################
@@ -399,10 +404,37 @@ function Unstructured()
 end
 const UN = Unstructured()
 
+
+"""
+    ACOV(c, action = 0)
+
+!!! warning
+    Experimental
+    
+Augmented (adjusted) covariance. Add additional correlations to existed R-part of variance covariance matrix.
+Can be used with `AR` or `CS` types. 
+
+`action` if existed covariance not equal sero:
+
+* 0 - add
+* 1 - replace
+* 2 - do nothing (use existed value)
+* 3 - warning and add
+* 4 - warning and replace
+* 5 - warning and use existed value
+* other - error
+
+"""
+function ACOV(c; action = 0)
+    CovarianceType(ACOV_(c.s, action))
+end
+
 function RZero()
     CovarianceType(ZERO(), false)
 end
 
+#######################################################################################
+
 function covstrparam(ct::SI_, ::Int)::Tuple{Int, Int}
     return (1, 0)
 end
@@ -446,6 +478,10 @@ function covstrparam(ct::SPPOWD_, ::Int)::Tuple{Int, Int}
     return (2, 1)
 end
 
+function covstrparam(ct::ACOV_{<:Union{CS_, AR_, SPPOW_}}, ::Int)::Tuple{Int, Int}
+    return (0, 1)
+end
+
 function covstrparam(ct::ZERO, ::Int)::Tuple{Int, Int}
     return (0, 0)
 end
@@ -555,6 +591,10 @@ function rcoefnames(s, t, ct::UN_)
     return v
 end
 
+function rcoefnames(s, t, ct::ACOV_{<: Union{CS_, AR_}})
+    return ["ρ "]
+end
+
 function rcoefnames(s, t, ct::AbstractCovarianceType)
     v = Vector{String}(undef, t)
     v .= "Val "
@@ -643,6 +683,9 @@ end
 function Base.show(io::IO, ct::UN_)
     print(io, "UN")
 end
+function Base.show(io::IO, ct::ACOV_)
+    print(io, "ACOV(", ct.c, ")")
+end
 function Base.show(io::IO, ct::ZERO)
     print(io, "No effect")
 end