matrix wts, docs, swc cov type, chnges in rmat!

PharmCat · Aug 21, 2024 · 5f05ec7 · 5f05ec7
1 parent 0989e76
commit 5f05ec7
Show file tree

Hide file tree

Showing 10 changed files with 183 additions and 39 deletions.
diff --git a/Project.toml b/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.15.1"
+version = "0.16.0"
 
 [deps]
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"

diff --git a/docs/src/custom.md b/docs/src/custom.md
@@ -26,10 +26,10 @@ function Metida.gmat!(mx, θ, ::YourCovarianceStruct)
 end
 ```
 
-Function `rmat!` have 4 arguments and add repeated effect to V': V = V' + R (so V = Z * G * Z' + R), `mx` - V' matrix, `θ` - theta vector for this effect, `rz` - subject effect matrix, `ct` - your covariance type object. For example, `rmat!` for Heterogeneous Toeplitz Parameterized structure is specified bellow (`TOEPHP_  <: AbstractCovarianceType`).
+Function `rmat!` have 5 arguments and add repeated effect to V': V = V' + R (so V = Z * G * Z' + R), `mx` - V' matrix, `θ` - theta vector for this effect, `rz` - subject effect matrix, `ct` - your covariance type object, `sb` = block number. For example, `rmat!` for Heterogeneous Toeplitz Parameterized structure is specified bellow (`TOEPHP_  <: AbstractCovarianceType`).
 
 ```
-function Metida.rmat!(mx, θ, rz, ct::TOEPHP_)
+function Metida.rmat!(mx, θ, rz, ct::TOEPHP_, ::Int)
     l     = size(rz, 2)
     vec   = rz * (θ[1:l])
     s   = size(mx, 1)
@@ -123,7 +123,7 @@ Metida.fit!(lmm)
 
 # for R matrix
 
-function Metida.rmat!(mx, θ, rz, ::CustomCovarianceStructure)
+function Metida.rmat!(mx, θ, rz, ::CustomCovarianceStructure, ::Int)
     vec = Metida.tmul_unsafe(rz, θ)
     rn    = size(mx, 1)
     if rn > 1

diff --git a/docs/src/details.md b/docs/src/details.md
@@ -13,6 +13,33 @@ In matrix notation a mixed effect model can be represented as:
 y = X\beta + Zu + \epsilon
 ```
 
+where:
+
+
+```math
+\epsilon \sim N(0, R)
+
+\\
+
+u \sim N(0, G)
+
+\\
+
+y \sim N(X\beta, V)
+
+```
+
+where V depends on covariance sructure and parameters ``\theta``:
+
+```math
+V = CovStruct(\theta)
+
+```
+
+The unknown parameters include the regression parameters in ``\beta`` and covariance parameters in ``\theta``.
+
+Estimation of these model parameters relies on the use of a Newton-Ralphson (by default) algorithm. When we use either algorithm for finding REML solutions, we need to compute ``V^{-1}`` and its derivatives with respect to ``\theta``, which are computationally difficult for large ``n``, therefor SWEEP (see https://github.com/joshday/SweepOperator.jl) algorithm used to meke oprtimization less computationaly expensive.
+
 #### V
 
 ```math
@@ -33,7 +60,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,16 +78,36 @@ 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
+#### [Weights](@id weights_header)
 
 If weights defined:
 
 ```math
+
+W_{i} = diag(wts_{i})
+
+\\
+
 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.
+
+
+If wts is matrix then:
+
+
+```math
+
+W_{i} = wts_{i}
+
+\\
+
+V_{i} = Z_{i} G Z_i'+  R_{i} \circ W_{i}
+```
+
+where ``\circ`` - element-wise multiplication.
 
-where ```W``` - diagonal matrix of weights.
 
 
 #### Multiple random and repeated effects

diff --git a/src/lmm.jl b/src/lmm.jl
@@ -11,7 +11,7 @@ struct ModelStructure
 end
 
 """
-    LMM(model, data; contrasts=Dict{Symbol,Any}(),  random::Union{Nothing, VarEffect, Vector{VarEffect}} = nothing, repeated::Union{Nothing, VarEffect} = nothing, wts::Union{Nothing, AbstractVector, AbstractString, Symbol} = nothing)
+    LMM(model, data; contrasts=Dict{Symbol,Any}(),  random::Union{Nothing, VarEffect, Vector{VarEffect}} = nothing, repeated::Union{Nothing, VarEffect} = nothing, wts::Union{Nothing, AbstractVector, AbstractMatrix, AbstractString, Symbol} = nothing)
 
 Make Linear-Mixed Model object.
 
@@ -27,6 +27,10 @@ Make Linear-Mixed Model object.
 
 `wts`: regression weights (residuals).
 
+Weigts can be set as `Symbol` or `String`, in this case weights taken from tabular data.
+If weights is vector then this vector applyed to R-side part of covariance matrix (see [Weights details](@ref weights_header)).
+If weights is matrix then R-side part of covariance matrix multiplied by corresponding part of weight-matrix.
+
 See also: [`@lmmformula`](@ref)
 """
 struct LMM{T <: AbstractFloat, W <: Union{LMMWts, Nothing}} <: MetidaModel
@@ -57,7 +61,11 @@ struct LMM{T <: AbstractFloat, W <: Union{LMMWts, Nothing}} <: MetidaModel
         log::Vector{LMMLogMsg}) where T where W <: Union{LMMWts, Nothing}
         new{T, W}(model, f, modstr, covstr, data, dv, nfixed, rankx, result, maxvcbl, wts, log)
     end
-    function LMM(model, data; contrasts=Dict{Symbol,Any}(),  random::Union{Nothing, VarEffect, Vector{VarEffect}} = nothing, repeated::Union{Nothing, VarEffect, Vector{VarEffect}} = nothing, wts::Union{Nothing, AbstractVector, AbstractString, Symbol} = nothing)
+    function LMM(model, data; 
+        contrasts=Dict{Symbol,Any}(),  
+        random::Union{Nothing, VarEffect, Vector{VarEffect}} = nothing, 
+        repeated::Union{Nothing, VarEffect, Vector{VarEffect}} = nothing, 
+        wts::Union{Nothing, AbstractVector, AbstractMatrix, AbstractString, Symbol} = nothing)
         #need check responce - Float
         if !Tables.istable(data) error("Data not a table!") end
         if repeated === nothing && random === nothing
@@ -122,12 +130,22 @@ struct LMM{T <: AbstractFloat, W <: Union{LMMWts, Nothing}} <: MetidaModel
             if wts isa Symbol
                 wts = Tables.getcolumn(data, wts)
             end
-            if length(lmmdata.yv) == length(wts)
-                if any(x -> x <= zero(x), wts) error("Only cases with positive weights allowed!") end
-                lmmwts = LMMWts(wts, covstr.vcovblock)
-            else
-                @warn "wts count not equal observations count! wts not used."
-                lmmwts = nothing
+            if wts isa AbstractVector
+                if length(lmmdata.yv) == length(wts)
+                    if any(x -> x <= zero(x), wts) error("Only cases with positive weights allowed!") end
+                    lmmwts = LMMWts(wts, covstr.vcovblock)
+                else
+                    @warn "wts count not equal observations count! wts not used."
+                    lmmwts = nothing
+                end
+            elseif wts isa AbstractMatrix
+                if length(lmmdata.yv) == LinearAlgebra.checksquare(wts)
+                    if any(x -> x <= zero(x), wts) error("Only positive values allowed!") end
+                    lmmwts = LMMWts(wts, covstr.vcovblock)
+                else
+                    @warn "wts count not equal observations count! wts not used."
+                    lmmwts = nothing
+                end
             end
         end
 

diff --git a/src/lmmdata.jl b/src/lmmdata.jl
@@ -35,16 +35,23 @@ struct LMMDataViews{T<:AbstractFloat} <: AbstractLMMDataBlocks
     end
 end
 
-struct LMMWts{T<:AbstractFloat} 
-    sqrtwts::Vector{Vector{T}}
-    function LMMWts(sqrtwts::Vector{Vector{T}}) where T
+struct LMMWts{T} 
+    sqrtwts::Vector{T}
+    function LMMWts(sqrtwts::Vector{T}) where T
         new{T}(sqrtwts)
     end
-    function LMMWts(wts::Vector{T}, vcovblock) where T
+    function LMMWts(wts::AbstractVector{T}, vcovblock) where T
         sqrtwts = Vector{Vector{T}}(undef, length(vcovblock))
         for i in eachindex(vcovblock)
             sqrtwts[i] = @. inv(sqrt($(view(wts, vcovblock[i]))))
         end
         LMMWts(sqrtwts)
     end
+    function LMMWts(wts::AbstractMatrix{T}, vcovblock) where T
+        sqrtwts = Vector{Matrix{T}}(undef, length(vcovblock))
+        for i in eachindex(vcovblock)
+            sqrtwts[i] = wts[vcovblock[i], vcovblock[i]]
+        end
+        LMMWts(sqrtwts)
+    end
 end
diff --git a/src/rmat.jl b/src/rmat.jl
@@ -10,25 +10,40 @@
         zblock    = view(covstr.rz[j], block, :)
         @simd for i = 1:subjn(covstr, en, bi)
             sb = getsubj(covstr, en, bi, i)
-            rmat!(view(mx, sb, sb), view(θ, rθ[j]), view(zblock, sb, :), covstr.repeated[j].covtype.s)
+            rmat!(view(mx, sb, sb), view(θ, rθ[j]), view(zblock, sb, :), covstr.repeated[j].covtype.s, bi)
         end
     end
     mx
 end
 ################################################################################
-function rmat!(::Any, ::Any, ::Any,  ::AbstractCovarianceType)
+function rmat!(::Any, ::Any, ::Any,  ::AbstractCovarianceType, ::Int)
     error("No rmat! method defined for thit structure!")
 end
 #SI
-Base.@propagate_inbounds function rmat!(mx, θ, ::AbstractMatrix, ::SI_)
+Base.@propagate_inbounds function rmat!(mx, θ, ::AbstractMatrix, ::SI_, ::Int)
     val = θ[1]^2
     @inbounds @simd for i ∈ axes(mx, 1)
             mx[i, i] += val
     end
     mx
 end
+#SWC
+function rmat!(mx, θ, ::AbstractMatrix, ct::SWC_, sbj::Int)
+    s  = size(mx, 1)
+    de  = θ[1] ^ 2
+    if s > 1
+        for n = 1:s
+            @inbounds @simd for m = 1:n
+                mx[m, n] += de * ct.wtsb[sbj][m, n]
+            end
+        end
+    else
+        @inbounds mx[1, 1] += de * ct.wtsb[sbj][1, 1]
+    end
+    mx
+end
 #DIAG
-function rmat!(mx, θ, rz, ::DIAG_)
+function rmat!(mx, θ, rz, ::DIAG_, ::Int)
     for i ∈ axes(mx, 1)
         @inbounds @simd for c ∈ axes(θ, 1)
             mx[i, i] += rz[i, c] * θ[c] ^ 2
@@ -37,7 +52,7 @@ function rmat!(mx, θ, rz, ::DIAG_)
     mx
 end
 #AR
-function rmat!(mx, θ, ::AbstractMatrix, ::AR_)
+function rmat!(mx, θ, ::AbstractMatrix, ::AR_, ::Int)
     s  = size(mx, 1)
     de  = θ[1] ^ 2
     @inbounds @simd for m = 1:s
@@ -53,7 +68,7 @@ function rmat!(mx, θ, ::AbstractMatrix, ::AR_)
     mx
 end
 #ARH
-function rmat!(mx, θ, rz, ::ARH_)
+function rmat!(mx, θ, rz, ::ARH_, ::Int)
     vec = tmul_unsafe(rz, θ)
     s    = size(mx, 1)
     if s > 1
@@ -69,7 +84,7 @@ function rmat!(mx, θ, rz, ::ARH_)
     mx
 end
 #CS
-function rmat!(mx, θ, ::AbstractMatrix,  ::CS_)
+function rmat!(mx, θ, ::AbstractMatrix,  ::CS_, ::Int)
     s    = size(mx, 1)
     θsq   =  θ[1] * θ[1]
     θsqp  =  θsq  * θ[2]
@@ -86,7 +101,7 @@ function rmat!(mx, θ, ::AbstractMatrix,  ::CS_)
     mx
 end
 #CSH
-function rmat!(mx, θ, rz, ::CSH_)
+function rmat!(mx, θ, rz, ::CSH_, ::Int)
     vec = tmul_unsafe(rz, θ)
     s    = size(mx, 1)
     if s > 1
@@ -105,7 +120,7 @@ function rmat!(mx, θ, rz, ::CSH_)
 end
 ################################################################################
 #ARMA
-function rmat!(mx, θ, ::AbstractMatrix,  ::ARMA_)
+function rmat!(mx, θ, ::AbstractMatrix,  ::ARMA_, ::Int)
     s  = size(mx, 1)
     de  = θ[1] ^ 2
     @inbounds @simd for m = 1:s
@@ -123,7 +138,7 @@ function rmat!(mx, θ, ::AbstractMatrix,  ::ARMA_)
 end
 ################################################################################
 #TOEPP
-function rmat!(mx, θ, ::AbstractMatrix, ct::TOEPP_)
+function rmat!(mx, θ, ::AbstractMatrix, ct::TOEPP_, ::Int)
     de  = θ[1] ^ 2    #diagonal element
     s   = size(mx, 1) #size
     @inbounds @simd for i = 1:s
@@ -140,7 +155,7 @@ function rmat!(mx, θ, ::AbstractMatrix, ct::TOEPP_)
 end
 ################################################################################
 #TOEPHP
-function rmat!(mx, θ, rz, ct::TOEPHP_)
+function rmat!(mx, θ, rz, ct::TOEPHP_, ::Int)
     l     = size(rz, 2)
     vec   = rz * (θ[1:l])
     s   = size(mx, 1) #size
@@ -181,7 +196,7 @@ function edistance(mx::AbstractMatrix{T}, i::Int, j::Int) where T
 end
 ################################################################################
 #SPEXP
-function rmat!(mx, θ, rz,  ::SPEXP_)
+function rmat!(mx, θ, rz,  ::SPEXP_, ::Int)
     σ²    = θ[1]^2
     #θe    = exp(θ[2])
     θe    = θ[2]
@@ -202,7 +217,7 @@ function rmat!(mx, θ, rz,  ::SPEXP_)
 end
 ################################################################################
 #SPPOW
-function rmat!(mx, θ, rz,  ::SPPOW_)
+function rmat!(mx, θ, rz,  ::SPPOW_, ::Int)
     σ²    = θ[1]^2
     ρ     = θ[2]
     s    = size(mx, 1)
@@ -222,7 +237,7 @@ function rmat!(mx, θ, rz,  ::SPPOW_)
 end
 
 #SPGAU
-function rmat!(mx, θ, rz,  ::SPGAU_)
+function rmat!(mx, θ, rz,  ::SPGAU_, ::Int)
     σ²    = θ[1]^2
     #θe    = exp(θ[2])
     θe    = θ[2]
@@ -244,7 +259,7 @@ function rmat!(mx, θ, rz,  ::SPGAU_)
 end
 ################################################################################
 #SPEXPD cos(pidij)
-function rmat!(mx, θ, rz,  ::SPEXPD_)
+function rmat!(mx, θ, rz,  ::SPEXPD_, ::Int)
     σ²    = θ[2]^2
     σ²d   = θ[1]^2 + σ²
     θe    = θ[3]
@@ -263,7 +278,7 @@ function rmat!(mx, θ, rz,  ::SPEXPD_)
     mx
 end
 #SPPOWD
-function rmat!(mx, θ, rz,  ::SPPOWD_)
+function rmat!(mx, θ, rz,  ::SPPOWD_, ::Int)
     σ²    = θ[2]^2
     σ²d   = θ[1]^2 + σ²
     ρ     = θ[3]
@@ -281,7 +296,7 @@ function rmat!(mx, θ, rz,  ::SPPOWD_)
     mx
 end
 #SPGAUD
-function rmat!(mx, θ, rz,  ::SPGAUD_)
+function rmat!(mx, θ, rz,  ::SPGAUD_, ::Int)
     σ²    = θ[2]^2
     σ²d   = θ[1]^2 + σ²
     θe    = θ[3]
@@ -320,7 +335,7 @@ function unrmat(θ::AbstractVector{T}, rz) where T
     end
     Symmetric(mx)
 end
-function rmat!(mx, θ, rz::AbstractMatrix, ::UN_)
+function rmat!(mx, θ, rz::AbstractMatrix, ::UN_, ::Int)
     vec    = tmul_unsafe(rz, θ)
     rm     = size(mx, 1)
     rcov  = unrmat(θ, rz)

diff --git a/src/utils.jl b/src/utils.jl
@@ -272,10 +272,13 @@ function applywts!(::Any, ::Int, ::Nothing)
     nothing
 end
 
-function applywts!(V::AbstractMatrix, i::Int, wts::LMMWts)
+function applywts!(V::AbstractMatrix, i::Int, wts::LMMWts{<:Vector})
     mulβdαβd!(V, wts.sqrtwts[i])
 end
 
+function applywts!(V::AbstractMatrix, i::Int, wts::LMMWts{<:Matrix})
+    V .*= wts.sqrtwts[i]
+end
 
 #####################################################################