diff --git a/Project.toml b/Project.toml index 86d8102..3f7ad97 100644 --- 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 "] -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 index 7789d09..6ad8de7 100644 --- 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 index 4fceee5..a6959f5 100644 --- 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 index cdeb5c7..9edbcc2 100644 --- 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 index bf07fa7..98f6805 100644 --- 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 \ No newline at end of file diff --git a/src/rmat.jl b/src/rmat.jl index 958c3ac..5cac108 100644 --- 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 index 3a8efad..f2de44c 100644 --- 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 ##################################################################### diff --git a/src/varstruct.jl b/src/varstruct.jl index 49a426a..1fbb96b 100644 --- a/src/varstruct.jl +++ b/src/varstruct.jl @@ -394,6 +394,11 @@ struct CovStructure{T, T2} <: AbstractCovarianceStructure end esb = EffectSubjectBlock(sblock, nblock) ####################################################################### + # Modify repeated effect covariance type for some types + for r in repeated + applycovschema!(r.covtype.s, blocks) + end + ####################################################################### new{eltype(z), T2}(random, repeated, schema, rcnames, blocks, rn, rtn, rpn, z, esb, zrndur, rz, q, t, tr, tl, ct, emap, sn, maxn) end end diff --git a/src/vartypes.jl b/src/vartypes.jl index 1b4707b..d4abda2 100644 --- a/src/vartypes.jl +++ b/src/vartypes.jl @@ -5,6 +5,10 @@ ################################################################################ struct SI_ <: AbstractCovarianceType end +mutable struct SWC_{W<:AbstractMatrix, B<:Vector{<:AbstractMatrix}} <: AbstractCovarianceType + wtsm::W + wtsb::B +end struct DIAG_ <: AbstractCovarianceType end struct AR_ <: AbstractCovarianceType end struct ARH_ <: AbstractCovarianceType end @@ -64,6 +68,15 @@ function ScaledIdentity() CovarianceType(SI_()) end const SI = ScaledIdentity() + +# docs need +# Experimental +function ScaledWeightedCov(wtsm::AbstractMatrix{T}) where T + wtsb = Matrix{T}[] + CovarianceType(SWC_(wtsm, wtsb)) +end +const SWC = ScaledWeightedCov + """ Diag() @@ -362,6 +375,9 @@ end function covstrparam(ct::SI_, ::Int)::Tuple{Int, Int} return (1, 0) end +function covstrparam(ct::SWC_, ::Int)::Tuple{Int, Int} + return (1, 0) +end function covstrparam(ct::DIAG_, t::Int, )::Tuple{Int, Int} return (t, 0) end @@ -412,6 +428,9 @@ end function rcoefnames(s, t, ct::SI_) return ["σ² "] end +function rcoefnames(s, t, ct::SWC_) + return ["σ² "] +end function rcoefnames(s, t, ct::DIAG_) if isa(coefnames(s), AbstractArray{T,1} where T) l = length(coefnames(s)) else l = 1 end return fill!(Vector{String}(undef, l), "σ² ") .* string.(coefnames(s)) @@ -510,6 +529,21 @@ function rcoefnames(s, t, ct::AbstractCovarianceType) v .= "Val " return v end +################################################################################ +# APPLY COV SCHEMA +################################################################################ +function applycovschema!(::AbstractCovarianceType, ::Any) + nothing +end + +function applycovschema!(ct::SWC_{<:AbstractMatrix{T}}, vcovblock) where T + if length(ct.wtsb) == 0 + for i in eachindex(vcovblock) + push!(ct.wtsb, ct.wtsm[vcovblock[i], vcovblock[i]]) + end + end + ct +end ################################################################################ # SHOW @@ -524,6 +558,9 @@ end function Base.show(io::IO, ct::SI_) print(io, "SI") end +function Base.show(io::IO, ct::SWC_) + print(io, "SWC") +end function Base.show(io::IO, ct::DIAG_) print(io, "DIAG") end diff --git a/test/test.jl b/test/test.jl index afcf3be..0913efd 100644 --- a/test/test.jl +++ b/test/test.jl @@ -209,6 +209,18 @@ include("testdata.jl") fit!(lmm) @test Metida.m2logreml(lmm) ≈ 17.823729 atol=1E-6 # TEST WITH SPSS 28 + # Matrix wts + matwts = Symmetric(rand(size(df0, 1), size(df0, 1))) + lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0; + random = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.DIAG), + wts = matwts) + @test_nowarn fit!(lmm) + + # experimental weighted covariance + lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0; + repeated = Metida.VarEffect(Metida.@covstr(1|subject), Metida.SWC(matwts))) + @test_nowarn fit!(lmm) + # Repeated vector lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0; @@ -654,7 +666,7 @@ end reml = Metida.m2logreml(lmm) @test reml_c ≈ reml - 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