Merge pull request #35 from PharmCat/dev

Remove ModelFrame, add weights, Repeated effect vector first step

.github/workflows/CompatHelper.yml

@@ -1,7 +1,7 @@
 name: CompatHelper
 on:
   schedule:
-    - cron: '0 0 * * 1'
+    - cron: '0 0 1 * *'
   workflow_dispatch:
 jobs:
   CompatHelper:

.github/workflows/Tier1.yml

@@ -31,8 +31,8 @@ jobs:
       matrix:
         version:
           - '1.6'
-          - '1.7'
           - '1.8'
+          - '1'
         os:
           - ubuntu-latest
           - macOS-latest

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.14.9"
+version = "0.15.0"
 
 [deps]
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
@@ -20,13 +20,13 @@ StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
 
 [compat]
 DiffResults = "1"
-Distributions = "0.20, 0.21, 0.22, 0.23, 0.24, 0.25"
+Distributions = "0.21, 0.22, 0.23, 0.24, 0.25"
 ForwardDiff = "0.10"
 LineSearches = "7"
 MetidaBase = "0.11, 0.12"
 Optim = "1"
 ProgressMeter = "1"
-StatsBase = "0.29, 0.30, 0.31, 0.32, 0.33, 0.34"
+StatsBase = "0.30, 0.31, 0.32, 0.33, 0.34"
 StatsModels = "0.7"
 julia = "1"
 

docs/src/custom.md

@@ -91,7 +91,7 @@ Example:
 using Metida, DataFrames, CSV, CategoricalArrays
 
 ftdf = CSV.File(joinpath(dirname(pathof(Metida)), "..", "test", "csv",  "1fptime.csv"); types = [String, String, Float64, Float64]) |> DataFrame
-df0 = CSV.File(joinpath(dirname(pathof(Metida)), "..", "test", "csv",  "df0.csv"); types = [String, String, String, String,Float64, Float64]) |> DataFrame
+df0 = CSV.File(joinpath(dirname(pathof(Metida)), "..", "test", "csv",  "df0.csv"); types = [String, String, String, String,Float64, Float64, Float64]) |> DataFrame
 
 struct CustomCovarianceStructure <: Metida.AbstractCovarianceType end
 function Metida.covstrparam(ct::CustomCovarianceStructure, t::Int)::Tuple{Int, Int}

docs/src/details.md

@@ -33,7 +33,7 @@ logREML(\theta,\beta) = -\frac{N-p}{2} - \frac{1}{2}\sum_{i=1}^nlog|V_{\theta, i
 -\frac{1}{2}log|\sum_{i=1}^nX_i'V_{\theta, i}^{-1}X_i|-\frac{1}{2}\sum_{i=1}^n(y_i - X_{i}\beta)'V_{\theta, i}^{-1}(y_i - X_{i}\beta)
 ```
 
-Actually ```L(\theta) = -2logREML = L_1(\theta) + L_2(\theta) + L_3(\theta) + c`` used for optimization, where:
+Actually ```L(\theta) = -2logREML = L_1(\theta) + L_2(\theta) + L_3(\theta) + c``` used for optimization, where:
 
 ```math
 L_1(\theta) = \frac{1}{2}\sum_{i=1}^nlog|V_{i}| \\
@@ -51,6 +51,27 @@ L_3(\theta) = \frac{1}{2}\sum_{i=1}^n(y_i - X_{i}\beta)'V_i^{-1}(y_i - X_{i}\bet
 \mathcal{H}\mathcal{L}(\theta) =  \mathcal{H}L_1(\theta)  + \mathcal{H}L_2(\theta) +  \mathcal{H} L_3(\theta)
 ```
 
+#### Weights
+
+If weights defined:
+
+```math
+V_{i} = Z_{i} G Z_i'+ W^{- \frac{1}{2}}_i R_{i} W^{- \frac{1}{2}}_i
+```
+
+
+where ```W``` - diagonal matrix of weights.
+
+
+#### Multiple random and repeated effects
+
+If model include multiple effects ( with n random and m repeated effects) final V will be:
+
+```math
+V_{i} = Z_{i, 1} G_{1} Z_{i, 1}' + ... + Z_{i, n} G_{1} Z_{i, n}'+ W^{- \frac{1}{2}}_i ( R_{i, 1} + ... + R_{i, m}) W^{- \frac{1}{2}}_i
+```
+
+
 ##### Initial step
 
 Initial (first) step before optimization may be done:

docs/src/instanduse.md

@@ -63,8 +63,8 @@ Define `random` and `repreated` effects with [`Metida.VarEffect`](@ref) using [`
 right side is a effect itself. [`Metida.HeterogeneousCompoundSymmetry`](@ref) and [`Metida.Diag`](@ref) (Diagonal) in example bellow is a model of variance-covariance structure. See also [`Metida.@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.
-  Metida does not check this!
+
+    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. Metida does not check this!
 
 
 ```@example lmmexample
@@ -73,6 +73,15 @@ random = VarEffect(@covstr(formulation|subject), CSH),
 repeated = VarEffect(@covstr(formulation|subject), DIAG));
 ```
 
+Also [`Metida.@lmmformula`](@ref) macro can be used:
+
+```julia
+lmm = LMM(@lmmformula(var~sequence+period+formulation,
+    random = formulation|subject:CSH,
+    repeated = formulation|subject:DIAG),
+    df)
+```
+
 #### Step 3: Fit
 
 Just fit the model.
@@ -87,6 +96,23 @@ fit!(lmm)
 lmm.log
 ```
 
+#### Confidence intervals for coefficients
+
+
+```@example lmmexample
+confint(lmm)
+```
+
+!!! note
+
+    Satterthwaite approximation for the denominator degrees of freedom used by default.
+
+#### StatsBsae API
+
+StatsBsae API implemented: [`Metida.islinear`](@ref), [`Metida.confint`](@ref), [`Metida.coef`](@ref), [`Metida.coefnames`](@ref), [`Metida.dof_residual`](@ref), [`Metida.dof`](@ref), [`Metida.loglikelihood`](@ref), [`Metida.aic`](@ref), [`Metida.bic`](@ref), [`Metida.aicc`](@ref),  [`Metida.isfitted`](@ref), [`Metida.vcov`](@ref), [`Metida.stderror`](@ref), [`Metida.modelmatrix`](@ref), [`Metida.response`](@ref), [`Metida.crossmodelmatrix`](@ref), [`Metida.coeftable`](@ref), [`Metida.responsename`](@ref)
+
+
+
 ##### Type III Tests of Fixed Effects
 
 !!! warning

src/dof_contain.jl

@@ -41,8 +41,8 @@ end
     Minimum returned. If no random effect found  N - rank(XZ) returned.
 """
 function dof_contain(lmm, i)
-    ind  = lmm.mm.assign[i]
-    sym  = StatsModels.termvars(lmm.mf.f.rhs.terms[ind])
+    ind  = lmm.modstr.assign[i]
+    sym  = StatsModels.termvars(lmm.f.rhs.terms[ind])
     rr   = Vector{Int}(undef, 0)
     for r = 1:length(lmm.covstr.random)
         if length(intersect(sym, StatsModels.termvars(lmm.covstr.random[r].model))) > 0
@@ -57,12 +57,12 @@ function dof_contain(lmm, i)
 end
 
 function dof_contain(lmm)
-    dof   = zeros(Int, length(lmm.mm.assign))
+    dof   = zeros(Int, length(lmm.modstr.assign))
     rrt   = zeros(Int, length(lmm.covstr.random))
     rz   = 0
-    for i = 1:length(lmm.mm.assign)
-        ind  = lmm.mm.assign[i]
-        sym  = StatsModels.termvars(lmm.mf.f.rhs.terms[ind])
+    for i = 1:length(lmm.modstr.assign)
+        ind  = lmm.modstr.assign[i]
+        sym  = StatsModels.termvars(lmm.f.rhs.terms[ind])
         rr   = Vector{Int}(undef, 0)
         for r = 1:length(lmm.covstr.random)
             if length(intersect(sym, StatsModels.termvars(lmm.covstr.random[r].model))) > 0
@@ -87,7 +87,7 @@ end
 
 """
 function dof_contain_f(lmm, i)
-    sym  = StatsModels.termvars(lmm.mf.f.rhs.terms[i])
+    sym  = StatsModels.termvars(lmm.f.rhs.terms[i])
     rr   = Vector{Int}(undef, 0)
     for r = 1:length(lmm.covstr.random)
         if length(intersect(sym, StatsModels.termvars(lmm.covstr.random[r].model))) > 0

src/fit.jl

@@ -51,6 +51,7 @@ Fit LMM model.
 * `varlinkf` - :exp / :sq / :identity [ref](@ref varlink_header)
 * `rholinkf` - :sigm / :atan / :sqsigm / :psigm
 * `aifirst` - first iteration with AI-like method - :default / :ai / :score
+* `aifmax` - maximum pre-optimization steps
 * `g_tol` - absolute tolerance in the gradient
 * `x_tol` - absolute tolerance of theta vector
 * `f_tol` - absolute tolerance in changes of the REML

src/gmat.jl

@@ -89,11 +89,12 @@ end
 #CS
 function gmat!(mx, θ, ::CS_)
     s = size(mx, 1)
-    mx .= θ[1]^2
+    θ₁² = θ[1]^2
+    mx .= θ₁²
     if s > 1
-        mxθ2 = θ[1]^2 * θ[2]
+        mxθ2 = θ₁² * θ[2]
         for n = 2:s
-            @inbounds @simd for m = 1:n-1
+            @inbounds @simd for m = 1:n - 1
                 mx[m, n] = mxθ2
             end
         end
@@ -212,8 +213,8 @@ function gmat!(mx, θ, ::UN_)
     end
     if s > 1
         for n = 2:s
-            @inbounds @simd for m = 1:n-1
-                mx[m, n] = mx[m, m] * mx[n, n] * θ[s+tpnum(m, n, s)]
+            @inbounds @simd for m = 1:n - 1
+                mx[m, n] = mx[m, m] * mx[n, n] * θ[s + tpnum(m, n, s)]
             end
         end
     end

src/linearalgebra.jl

@@ -139,6 +139,20 @@ Change θ.
     end
     θ
 end
+# Diagonal(b) * A * Diagonal(b) - chnage only A upper triangle 
+@noinline function mulβdαβd!(A::AbstractMatrix, b::AbstractVector)
+    q = size(A, 1)
+    p = size(A, 2)
+    if !(q == p == length(b)) throw(DimensionMismatch("size(A, 1) and size(A, 2) should be equal length(b)")) end
+    for n in 1:p
+        @simd for m in 1:n
+            @inbounds A[m, n] *= b[m] * b[n]
+        end
+    end
+    A
+end
+
+
 ################################################################################
 @inline function tmul_unsafe(rz, θ::AbstractVector{T}) where T
     vec = zeros(T, size(rz, 1))