diff --git a/Project.toml b/Project.toml index 005d3d62..0597345a 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.13.0" +version = "0.14.0" [deps] DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" diff --git a/change.log b/change.log index fca86da8..4d410f66 100644 --- a/change.log +++ b/change.log @@ -1,5 +1,9 @@ -v0.13 +v0.14 + * Unstructured Covariance (UN) + * initial functions for analytical gradient vector + * docs +v0.13 * Bootstrap * MI * Changes in types diff --git a/docs/src/index.md b/docs/src/index.md index 5a4681d9..c38faa25 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -27,6 +27,7 @@ Implemented covariance structures: * Spatial Exponential (SPEXP) * Spatial Power (SPPOW) * Spatial Gaussian (SPGAU) + * Unstructured (UN) * Custom Covariance Type ## Limitations diff --git a/src/Metida.jl b/src/Metida.jl index f5742b8f..fd5f5920 100644 --- a/src/Metida.jl +++ b/src/Metida.jl @@ -16,7 +16,7 @@ import MetidaBase.PrettyTables: TextFormat, pretty_table, tf_borderless, ft_prin import LinearAlgebra:checksquare, BlasFloat import StatsModels: @formula, termvars, ModelFrame, FunctionTerm, AbstractTerm, CategoricalTerm, AbstractContrasts, ConstantTerm, InterceptTerm, Term, InteractionTerm, FormulaTerm, ModelMatrix, schema, apply_schema, MatrixTerm, modelcols import StatsBase: fit, fit!, coef, coefnames, confint, nobs, dof_residual, dof, loglikelihood, aic, bic, aicc, isfitted, vcov, mean, var, stderror, modelmatrix, response, CoefTable, coeftable -import Base:show, rand, ht_keyindex +import Base:show, rand, ht_keyindex, getproperty import Random: default_rng, AbstractRNG, rand! export @formula, @covstr, diff --git a/src/gmat.jl b/src/gmat.jl index 1258729c..893d0c7b 100644 --- a/src/gmat.jl +++ b/src/gmat.jl @@ -8,7 +8,7 @@ function gmat_switch!(G, θ, covstr, r) G end function gmatvec(θ::Vector{T}, covstr) where T - gt = Tuple(Symmetric(Matrix{T}(undef, covstr.q[r], covstr.q[r])) for r in 1:covstr.rn) + gt = [Symmetric(Matrix{T}(undef, covstr.q[r], covstr.q[r])) for r in 1:covstr.rn] #<: fix to vec for r = 1:covstr.rn fill!(gt[r].data, zero(T)) if covstr.random[r].covtype.z @@ -18,9 +18,8 @@ function gmatvec(θ::Vector{T}, covstr) where T gt end ################################################################################ -function zgz_base_inc!(mx::AbstractArray{T}, θ::AbstractArray{T}, covstr, block, sblock) where T +function zgz_base_inc!(mx::AbstractArray, θ::AbstractArray{T}, covstr, block, sblock) where T if covstr.random[1].covtype.z - #length of random to covstr for r = 1:covstr.rn G = fill!(Symmetric(Matrix{T}(undef, covstr.q[r], covstr.q[r])), zero(T)) gmat_switch!(G.data, θ, covstr, r) @@ -30,24 +29,17 @@ function zgz_base_inc!(mx::AbstractArray{T}, θ::AbstractArray{T}, covstr, block end end end - # fill!(mx, zero(T)) - #end mx end -function zgz_base_inc!(mx::AbstractArray{T}, G, θ::AbstractArray{T}, covstr, block, sblock) where T +function zgz_base_inc!(mx::AbstractArray, G, θ, covstr, block, sblock) if covstr.random[1].covtype.z - #length of random to covstr for r = 1:covstr.rn - #G = fill!(Symmetric(Matrix{T}(undef, covstr.q[r], covstr.q[r])), zero(T)) - #gmat_switch!(G.data, θ, covstr, r) zblock = view(covstr.z, block, covstr.zrndur[r]) @inbounds for i = 1:length(sblock[r]) mulαβαtinc!(view(mx, sblock[r][i], sblock[r][i]), view(zblock, sblock[r][i], :), G[r]) end end end - # fill!(mx, zero(T)) - #end mx end ################################################################################ @@ -64,18 +56,15 @@ Base.@propagate_inbounds function gmat!(mx, θ, ::SI_) @inbounds @simd for i = 1:size(mx, 1) mx[i, i] = val end - nothing + mx end #DIAG function gmat!(mx, θ, ::DIAG_) @inbounds @simd for i = 1:size(mx, 1) mx[i, i] = θ[i] ^ 2 end - nothing + mx end -#function gmat_vc!(mx, θ::Vector{T}, ::Int, ::CovarianceType) where T -# nothing -#end #AR function gmat!(mx, θ, ::AR_) de = θ[1] ^ 2 @@ -90,7 +79,7 @@ function gmat!(mx, θ, ::AR_) end end end - nothing + mx end #ARH function gmat!(mx, θ, ::ARH_) @@ -108,7 +97,7 @@ function gmat!(mx, θ, ::ARH_) @inbounds @simd for m = 1:s mx[m, m] = mx[m, m] * mx[m, m] end - nothing + mx end #CS function gmat!(mx, θ, ::CS_) @@ -121,7 +110,7 @@ function gmat!(mx, θ, ::CS_) end end end - nothing + mx end #CSH function gmat!(mx, θ, ::CSH_) @@ -139,7 +128,7 @@ function gmat!(mx, θ, ::CSH_) @inbounds @simd for m = 1:s mx[m, m] = mx[m, m] * mx[m, m] end - nothing + mx end ################################################################################ #ARMA @@ -156,7 +145,7 @@ function gmat!(mx, θ, ::ARMA_) end end end - nothing + mx end #TOEP function gmat!(mx, θ, ::TOEP_) @@ -172,7 +161,7 @@ function gmat!(mx, θ, ::TOEP_) end end end - nothing + mx end function gmat!(mx, θ, ct::TOEPP_) de = θ[1] ^ 2 #diagonal element @@ -187,7 +176,7 @@ function gmat!(mx, θ, ct::TOEPP_) end end end - nothing + mx end #TOEPH function gmat!(mx, θ, ::TOEPH_) @@ -205,7 +194,7 @@ function gmat!(mx, θ, ::TOEPH_) @inbounds @simd for m = 1:s mx[m, m] = mx[m, m] * mx[m, m] end - nothing + mx end #TOEPHP function gmat!(mx, θ, ct::TOEPHP_) @@ -223,7 +212,7 @@ function gmat!(mx, θ, ct::TOEPHP_) @inbounds @simd for m = 1:s mx[m, m] = mx[m, m] * mx[m, m] end - nothing + mx end #UN function gmat!(mx, θ, ::UN_) @@ -242,7 +231,7 @@ function gmat!(mx, θ, ::UN_) v = mx[m, m] mx[m, m] = v * v end - nothing + mx end function tpnum(m, n, s) @@ -252,3 +241,44 @@ function tpnum(m, n, s) end b -= s - n end + +################################################################################ + +# Grads +#= +function sim_gmat(s, θ::AbstractVector{T}, ct::AbstractCovarianceType) where T + mxt = zeros(T, s, s) + Symmetric(gmat!(mxt, θ, ct)) +end + +function grad_gmat!(s, θ, ct::AbstractCovarianceType) + mxg = zeros(s, s, length(θ)) + gm = ForwardDiff.jacobian!(mxg, x-> sim_gmat(s, x, ct), θ) + mxg +end + +function grad_gmatvec(θ::Vector{T}, covstr) where T + gt = [Matrix{T}(undef, covstr.q[r], covstr.q[r], length(θ)) for r in 1:covstr.rn] #<: fix to vec + for r = 1:covstr.rn + if covstr.random[r].covtype.z + gt[r] = grad_gmat!(covstr.q[r], view(θ, covstr.tr[r]), covstr.random[r].covtype.s) + end + end + gt +end + +function grad_zgz_base_inc!(mx::AbstractArray, G, θ::AbstractArray{T}, covstr, block, sblock) where T + if covstr.random[1].covtype.z + gl = size(mx, 3) + for r = 1:covstr.rn + zblock = view(covstr.z, block, covstr.zrndur[r]) + @inbounds for i = 1:length(sblock[r]) + for gn = 1:gl + mulαβαtinc!(view(mx, sblock[r][i], sblock[r][i], gn), view(zblock, sblock[r][i], :), view(G[r], :, :, gn)) + end + end + end + end + mx +end +=# diff --git a/src/linearalgebra.jl b/src/linearalgebra.jl index e9318cb8..47c876e9 100644 --- a/src/linearalgebra.jl +++ b/src/linearalgebra.jl @@ -19,7 +19,7 @@ end Change θ (only upper triangle). B is symmetric. """ -function mulαβαtinc!(θ::AbstractMatrix{T}, A::AbstractMatrix, B::AbstractMatrix) where T +function mulαβαtinc!(θ::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix) axb = axes(B, 1) sa = size(A, 1) for m ∈ 1:sa @@ -38,7 +38,7 @@ end Change θ (only upper triangle). B is symmetric. """ -function mulαβαtinc!(θ::AbstractMatrix{T}, A::AbstractMatrix, B::AbstractMatrix, alpha) where T +function mulαβαtinc!(θ::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, alpha) if !(size(B, 1) == size(B, 2) == size(A, 2)) || !(size(A, 1) == size(θ, 1) == size(θ, 2)) throw(ArgumentError("Wrong dimentions!")) end axb = axes(B, 1) sa = size(A, 1) @@ -58,7 +58,7 @@ end Change θ (only upper triangle). B is symmetric. """ -function mulαβαtinc!(θ::AbstractVector{T}, A::AbstractMatrix, B::AbstractMatrix, a::AbstractVector, b::AbstractVector, alpha) where T +function mulαβαtinc!(θ::AbstractVector, A::AbstractMatrix, B::AbstractMatrix, a::AbstractVector, b::AbstractVector, alpha) if !(size(B, 2) == length(a) == length(b)) || size(B, 1) != size(A, 2) || size(A, 1) != length(θ) throw(ArgumentError("Wrong dimentions.")) end axb = axes(B, 1) sa = size(A, 1) diff --git a/src/lmm.jl b/src/lmm.jl index 7365f3bc..af4981d1 100644 --- a/src/lmm.jl +++ b/src/lmm.jl @@ -287,3 +287,14 @@ Return fitting log. function getlog(lmm::LMM) lmm.log end + +################################################################################ + +function Base.getproperty(x::LMM, s::Symbol) + if s == :θ + return x.result.theta + elseif s == :β + return x.result.beta + end + getfield(x, s) +end diff --git a/src/reml.jl b/src/reml.jl index d47370e9..732c8c5a 100644 --- a/src/reml.jl +++ b/src/reml.jl @@ -33,12 +33,7 @@ end ################################################################################ # REML without provided β ################################################################################ -#= -function reml_sweep_β(lmm, θ::Vector{T}; syrkblas::Bool = false) where T <: Number - data = LMMDataViews(lmm) - reml_sweep_β(lmm, data, θ; syrkblas = syrkblas) -end -=# + function reml_sweep_β(lmm, data, θ::Vector{T}; syrkblas::Bool = false) where T n = length(lmm.covstr.vcovblock) N = length(lmm.data.yv) @@ -110,8 +105,9 @@ function reml_sweep_β(lmm, data, θ::Vector{T}; syrkblas::Bool = false) where T if issuccess(cθs₂) # β calculation mul!(β, inv(cθs₂), βm) - # final θ₂ + # final θ₂ logdetθ₂ = logdet(cθs₂) + # θ₃ @inbounds @simd for i = 1:n θ₃ += mulθ₃(data.yv[i], data.xv[i], β, V⁻¹[i]) end @@ -119,20 +115,13 @@ function reml_sweep_β(lmm, data, θ::Vector{T}; syrkblas::Bool = false) where T β .= NaN return Inf, β, Inf, Inf, false end - # θ₃ - return θ₁ + logdetθ₂ + θ₃ + c, β, θs₂, θ₃, noerror #REML, β, iC, θ₃, errors end ################################################################################ # REML with provided β ################################################################################ -#= -function reml_sweep_β(lmm, θ::Vector{T}, β::Vector) where T <: Number - data = LMMDataViews(lmm) - reml_sweep_β(lmm, data, θ, β) -end -=# + function core_sweep_β(lmm, data, θ::Vector{T}, β, n) where T ncore = min(num_cores(), n, 16) accθ₁ = zeros(T, ncore) diff --git a/src/rmat.jl b/src/rmat.jl index 4ff9561b..5429f2a4 100644 --- a/src/rmat.jl +++ b/src/rmat.jl @@ -25,7 +25,7 @@ Base.@propagate_inbounds function rmat!(mx, θ, ::AbstractMatrix, ::SI_) @inbounds @simd for i ∈ axes(mx, 1) mx[i, i] += θsq end - nothing + mx end #DIAG function rmat!(mx, θ, rz, ::DIAG_) @@ -34,7 +34,7 @@ function rmat!(mx, θ, rz, ::DIAG_) mx[i, i] += rz[i, c] * θ[c] * θ[c] end end - nothing + mx end #AR function rmat!(mx, θ, ::AbstractMatrix, ::AR_) @@ -50,7 +50,7 @@ function rmat!(mx, θ, ::AbstractMatrix, ::AR_) end end end - nothing + mx end #ARH function rmat!(mx, θ, rz, ::ARH_) @@ -66,7 +66,7 @@ function rmat!(mx, θ, rz, ::ARH_) @inbounds for m ∈ axes(mx, 1) mx[m, m] += vec[m] * vec[m] end - nothing + mx end #CS function rmat!(mx, θ, ::AbstractMatrix, ::CS_) @@ -83,19 +83,10 @@ function rmat!(mx, θ, ::AbstractMatrix, ::CS_) end end end - nothing + mx end #CSH function rmat!(mx, θ, rz, ::CSH_) - #vec = rz * (θ[1:end-1]) - #= - vec = zeros(T, size(rz, 1)) - #=@turbo=# @inbounds for r ∈ axes(rz, 1) - for i ∈ axes(rz, 2) - vec[r] += rz[r, i] * θ[i] - end - end - =# vec = tmul_unsafe(rz, θ) rn = size(mx, 1) if rn > 1 @@ -108,7 +99,7 @@ function rmat!(mx, θ, rz, ::CSH_) #=@turbo=# @inbounds for m ∈ axes(mx, 1) mx[m, m] += vec[m] * vec[m] end - nothing + mx end ################################################################################ #ARMA @@ -125,7 +116,7 @@ function rmat!(mx, θ, ::AbstractMatrix, ::ARMA_) end end end - nothing + mx end ################################################################################ #TOEPP @@ -142,7 +133,7 @@ function rmat!(mx, θ, ::AbstractMatrix, ct::TOEPP_) end end end - nothing + mx end ################################################################################ #TOEPHP @@ -160,7 +151,7 @@ function rmat!(mx, θ, rz, ct::TOEPHP_) @inbounds @simd for m = 1:s mx[m, m] += vec[m] * vec[m] end - nothing + mx end ################################################################################ #= @@ -199,7 +190,7 @@ function rmat!(mx, θ, rz, ::SPEXP_) end end end - nothing + mx end ################################################################################ #SPPOW @@ -219,7 +210,7 @@ function rmat!(mx, θ, rz, ::SPPOW_) end end end - nothing + mx end #SPGAU @@ -241,7 +232,7 @@ function rmat!(mx, θ, rz, ::SPGAU_) end end end - nothing + mx end ################################################################################ #SPEXPD cos(pidij) @@ -261,7 +252,7 @@ function rmat!(mx, θ, rz, ::SPEXPD_) end end end - nothing + mx end #SPPOWD function rmat!(mx, θ, rz, ::SPPOWD_) @@ -279,7 +270,7 @@ function rmat!(mx, θ, rz, ::SPPOWD_) end end end - nothing + mx end #SPGAUD function rmat!(mx, θ, rz, ::SPGAUD_) @@ -298,11 +289,10 @@ function rmat!(mx, θ, rz, ::SPGAUD_) end end end - nothing + mx end #UN - function unrmat(θ::AbstractVector{T}, rz) where T rm = size(rz, 2) mx = zeros(T, rm, rm) @@ -326,5 +316,5 @@ function rmat!(mx, θ, rz::AbstractMatrix, ::UN_) rm = size(mx, 1) rcov = unrmat(θ, rz) mulαβαtinc!(mx, rz, rcov) - nothing + mx end diff --git a/src/utils.jl b/src/utils.jl index e6e03e12..cfc72127 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -237,8 +237,8 @@ function vmatrix!(V, G, θ, lmm::LMM, i::Int) zgz_base_inc!(V, G, θ, lmm.covstr, lmm.covstr.vcovblock[i], lmm.covstr.sblock[i]) rmat_base_inc!(V, θ[lmm.covstr.tr[end]], lmm.covstr, lmm.covstr.vcovblock[i], lmm.covstr.sblock[i]) end -function vmatrix(θ, lmm::LMM, i::Int) - V = zeros(length(lmm.covstr.vcovblock[i]), length(lmm.covstr.vcovblock[i])) +function vmatrix(θ::AbstractVector{T}, lmm::LMM, i::Int) where T + V = zeros(T, length(lmm.covstr.vcovblock[i]), length(lmm.covstr.vcovblock[i])) gvec = gmatvec(θ, lmm.covstr) vmatrix!(V, gvec, θ, lmm, i) Symmetric(V) @@ -259,6 +259,18 @@ function vmatrix(θ::Vector, covstr::CovStructure, i::Int) rmat_base_inc!(V, θ[covstr.tr[end]], covstr, covstr.vcovblock[i], covstr.sblock[i]) Symmetric(V) end +#= +function grad_vmatrix(θ::AbstractVector{T}, lmm::LMM, i::Int) + V = zeros(T, length(lmm.covstr.vcovblock[i]), length(lmm.covstr.vcovblock[i])) + gvec = gmatvec(θ, lmm.covstr) + vmatrix!(V, gvec, θ, lmm, i) + Symmetric(V) +end +=# + + + + function nblocks(lmm::LMM) return length(lmm.covstr.vcovblock) end diff --git a/src/vartypes.jl b/src/vartypes.jl index 6aa16233..174d630a 100644 --- a/src/vartypes.jl +++ b/src/vartypes.jl @@ -319,7 +319,13 @@ function SpatialGaussianD() CovarianceType(SPGAUD_()) end const SPGAUD = SpatialGaussianD() -#NOT IMPLEMENTED +""" + Unstructured() + +Unstructured covariance structure with `t*(t+1)/2-t` paremeters where `t` - number of factor levels, `t*(t+1)/2-2t` of them is covariance (ρ) patemeters. + +UN = Unstructured() +""" function Unstructured() CovarianceType(UN_()) end diff --git a/test/test.jl b/test/test.jl index f79b01b5..dc131ba7 100644 --- a/test/test.jl +++ b/test/test.jl @@ -46,6 +46,9 @@ include("testdata.jl") ############################################################################ l = [0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0] @test Metida.logreml(lmm) ≈ -8.120556322253035 atol=1E-6 + @test Metida.theta(lmm) ≈ [0.4473222800422779, 0.3673667558902593, 0.1850675552332174] + @test lmm.θ ≈ [0.4473222800422779, 0.3673667558902593, 0.1850675552332174] + @test lmm.β ≈ coef(lmm) @test isfitted(lmm) == true @test islinear(lmm) == true @test bic(lmm) ≈ 24.558878811225412 atol=1E-6