From 5327fc156afe6c3bbc76f6b591625d86d8442a99 Mon Sep 17 00:00:00 2001 From: PharmCat Date: Sat, 29 May 2021 17:38:31 +0300 Subject: [PATCH] v0.9.2 experimental spatial exp, internal changes, no varlink for not- var, rho parameters --- Project.toml | 2 +- change.log | 1 + docs/src/api.md | 5 ++++ src/fit.jl | 6 +++-- src/rmat.jl | 34 ++++++++++++++++++++++++++ src/utils.jl | 10 ++++---- src/varstruct.jl | 62 +++++++++++++++++++++++++++++++++--------------- test/test.jl | 8 +++++++ 8 files changed, 102 insertions(+), 26 deletions(-) diff --git a/Project.toml b/Project.toml index a9236850..54fd234e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Metida" uuid = "a1dec852-9fe5-11e9-361f-8d9fde67cfa2" authors = ["Vladimir Arnautov "] -version = "0.9.1" +version = "0.9.2" [deps] CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597" diff --git a/change.log b/change.log index 0512a873..7055017d 100644 --- a/change.log +++ b/change.log @@ -1,6 +1,7 @@ v0.9.2 * z matrix for random effect * contain ddf + * Experimental Spatial Exponential repeated structure v0.9.1 * update deps diff --git a/docs/src/api.md b/docs/src/api.md index f5862392..825e0d71 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -76,6 +76,11 @@ Metida.HeterogeneousToeplitz Metida.HeterogeneousToeplitzParameterized ``` +### Metida.SpatialExponential +```@docs +Metida.SpatialExponential +``` + ### Metida.anova ```@docs Metida.anova diff --git a/src/fit.jl b/src/fit.jl index 9b3f7657..cbf9fae4 100644 --- a/src/fit.jl +++ b/src/fit.jl @@ -89,9 +89,11 @@ function fit!(lmm::LMM{T}; end else initθ = sqrt(initvar(lmm.data.yv, lmm.mm.m)[1])/(length(lmm.covstr.random)+1) - θ .= initθ + #θ .= initθ for i = 1:length(θ) - if lmm.covstr.ct[i] == :rho + if lmm.covstr.ct[i] == :var + θ[i] = initθ + elseif lmm.covstr.ct[i] == :rho θ[i] = 1e-4 end end diff --git a/src/rmat.jl b/src/rmat.jl index f93c2c2e..20b926b9 100644 --- a/src/rmat.jl +++ b/src/rmat.jl @@ -24,6 +24,8 @@ function rmat_base_inc_b!(mx, θ, zrv, covstr) rmatp_toephp!(mx, θ, zrv, covstr.repeated.covtype.p) elseif covstr.repeated.covtype.s == :FUNC covstr.repeated.covtype.p.xmat!(mx, θ, zrv, covstr.repeated.covtype.p) + elseif covstr.repeated.covtype.s == :SPEXP + rmatp_spexp!(mx, θ, zrv, covstr.repeated.covtype.p) end end ################################################################################ @@ -161,3 +163,35 @@ function rmatp_toephp!(mx, θ, rz, p) end nothing end +################################################################################ +function edistance(i::AbstractVector{T1}, j::AbstractVector{T2}) where T1 where T2 + sum = zero(promote_type(T1, T2)) + for c = 1:length(i) + sum += (i[c]-j[c])^2 + end + return sqrt(sum) +end +function edistance(mx::AbstractMatrix{T}, i::Int, j::Int) where T + sum = zero(T) + for c = 1:size(mx, 2) + sum += (mx[i,c] - mx[j,c])^2 + end + return sqrt(sum) +end +################################################################################ +function rmatp_spexp!(mx, θ, rz, p) + σ² = θ[1]^2 + θ = exp(θ[2]) + rn = size(mx, 1) + @inbounds @simd for i = 1:size(mx, 1) + mx[i, i] += σ² + end + if rn > 1 + for m = 1:rn - 1 + @inbounds @simd for n = m + 1:rn + mx[m, n] += σ² * exp(-edistance(rz, m, n) / θ) + end + end + end + nothing +end diff --git a/src/utils.jl b/src/utils.jl index 1fd70b7a..4f9eef4d 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -70,7 +70,7 @@ function varlinkvecapply!(v, p; varlinkf = :exp, rholinkf = :sigm) for i = 1:length(v) if p[i] == :var v[i] = vlink(v[i]) - else + elseif p[i] == :rho if rholinkf == :sigm v[i] = rholinksigmoid(v[i]) elseif rholinkf == :atan @@ -88,7 +88,7 @@ function varlinkrvecapply!(v, p; varlinkf = :exp, rholinkf = :sigm) for i = 1:length(v) if p[i] == :var v[i] = vlinkr(v[i]) - else + elseif p[i] == :rho if rholinkf == :sigm v[i] = rholinksigmoidr(v[i]) elseif rholinkf == :atan @@ -107,7 +107,7 @@ function varlinkvecapply(v, p; varlinkf = :exp, rholinkf = :sigm) for i = 1:length(v) if p[i] == :var s[i] = vlink(v[i]) - else + elseif p[i] == :rho if rholinkf == :sigm s[i] = rholinksigmoid(v[i]) elseif rholinkf == :atan @@ -117,6 +117,8 @@ function varlinkvecapply(v, p; varlinkf = :exp, rholinkf = :sigm) elseif rholinkf == :psigm s[i] = rholinkpsigmoid(v[i]) end + else + s[i] = v[i] end end s @@ -126,7 +128,7 @@ function m2logreml(lmm) lmm.result.reml end function logreml(lmm) - -m2logreml(lmm)/2. + -m2logreml(lmm)/2 end ################################################################################ diff --git a/src/varstruct.jl b/src/varstruct.jl index 47d7d0ba..fb0fbad6 100644 --- a/src/varstruct.jl +++ b/src/varstruct.jl @@ -235,6 +235,27 @@ function HeterogeneousToeplitzParameterized(p::Int) end const TOEPHP(p) = HeterogeneousToeplitzParameterized(p) +""" + SpatialExponential() + +!!! warning + Experimental + +Spatian Exponential covariance structure. Used only for repeated effect. + +```math +R_{i,j} = \\sigma^{2} * exp(-dist(i,j)/exp(\\theta)) +``` + +where `dist` - Euclidean distance between row-vector of repeated effect matrix for subject `i` and `j`. + +SPEXP = SpatialExponential() + +""" +function SpatialExponential() + CovarianceType(:SPEXP) +end +const SPEXP = SpatialExponential() #Spatial Power ? #= """ @@ -311,37 +332,39 @@ CovarianceType(cm::AbstractCovmatMethod) = CovarianceType(:FUNC, cm) #CovarianceType(cmg::AbstractCovmatMethod, cmr::AbstractCovmatMethod) = FullCovarianceType(CovarianceType(:FUNC, cmg), CovarianceType(:FUNC, cmr)) -function covstrparam(ct::CovarianceType, t::Int, q::Int)::Tuple{Int, Int} +function covstrparam(ct::CovarianceType, t::Int, q::Int) if ct.s == :SI - return (1, 0) + return (1, 0, 0) elseif ct.s == :DIAG - return (t, 0) - elseif ct.s == :VC - return (q, 0) + return (t, 0, 0) + #elseif ct.s == :VC + # return (q, 0, q) elseif ct.s == :AR - return (1, 1) + return (1, 1, 0) elseif ct.s == :ARH - return (t, 1) + return (t, 1, 0) elseif ct.s == :ARMA - return (1, 2) + return (1, 2, 0) elseif ct.s == :CS - return (1, 1) + return (1, 1, 0) elseif ct.s == :CSH - return (t, 1) + return (t, 1, 0) elseif ct.s == :TOEP - return (1, t - 1) + return (1, t - 1, 0) elseif ct.s == :TOEPH - return (t, t - 1) + return (t, t - 1, 0) elseif ct.s == :TOEPP - return (1, ct.p - 1) + return (1, ct.p - 1, 0) elseif ct.s == :TOEPHP - return (t, ct.p - 1) + return (t, ct.p - 1, 0) elseif ct.s == :UN - return (t, t * (t + 1) / 2 - t) + return (t, t * (t + 1) / 2 - t, 0) elseif ct.s == :ZERO - return (0, 0) + return (0, 0, 0) elseif ct.s == :FUNC return ct.p.nparamf(t, q) + elseif ct.s == :SPEXP + return (1, 0, 1) else error("Unknown covariance type!") end @@ -467,7 +490,7 @@ struct CovStructure{T} <: AbstractCovarianceStructure ztemp = modelcols(MatrixTerm(schema[i]), data) q[i] = size(ztemp, 2) csp = covstrparam(random[i].covtype, q[i], random[i].p) - t[i] = csp[1] + csp[2] + t[i] = sum(csp) z = hcat(z, ztemp) fillur!(zrndur, i, q) fillur!(tr, i, t) @@ -488,7 +511,7 @@ struct CovStructure{T} <: AbstractCovarianceStructure sn[end] = size(subjz[end], 2) q[end] = size(rz, 2) csp = covstrparam(repeated.covtype, q[end], repeated.p) - t[end] = csp[1] + csp[2] + t[end] = sum(csp) tr[end] = UnitRange(sum(t[1:end-1]) + 1, sum(t[1:end-1]) + t[end]) updatenametype!(ct, rcnames, csp, schema[end], repeated.covtype.s) #Theta length @@ -544,7 +567,8 @@ end function updatenametype!(ct, rcnames, csp, schema, s) append!(ct, fill!(Vector{Symbol}(undef, csp[1]), :var)) append!(ct, fill!(Vector{Symbol}(undef, csp[2]), :rho)) - append!(rcnames, rcoefnames(schema, csp[1]+csp[2], s)) + if length(csp) == 3 append!(ct, fill!(Vector{Symbol}(undef, csp[3]), :theta)) end + append!(rcnames, rcoefnames(schema, sum(csp), s)) end ################################################################################ function rcoefnames(s, t, ve) diff --git a/test/test.jl b/test/test.jl index 6f2b3212..4eb7f117 100644 --- a/test/test.jl +++ b/test/test.jl @@ -429,3 +429,11 @@ end @test iA ≈ iAss atol=1E-6 @test iAs ≈ iAb atol=1E-6 end + +@testset " Experimental " begin + lmm = Metida.LMM(@formula(response ~ 1), ftdf; + repeated = Metida.VarEffect(Metida.@covstr(response+time|subject), Metida.SPEXP), + ) + Metida.fit!(lmm) + @test Metida.m2logreml(lmm) ≈ 1528.7150702624508 atol=1E-6 +end