v0.9.2 experimental spatial exp, internal changes, no varlink for not…

…- var, rho parameters

Project.toml

@@ -1,7 +1,7 @@
 name = "Metida"
 uuid = "a1dec852-9fe5-11e9-361f-8d9fde67cfa2"
 authors = ["Vladimir Arnautov <[email protected]>"]
-version = "0.9.1"
+version = "0.9.2"
 
 [deps]
 CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597"

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

docs/src/api.md

@@ -76,6 +76,11 @@ Metida.HeterogeneousToeplitz
 Metida.HeterogeneousToeplitzParameterized
 ```
 
+### Metida.SpatialExponential
+```@docs
+Metida.SpatialExponential
+```
+
 ### Metida.anova
 ```@docs
 Metida.anova

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

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

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
 ################################################################################
 

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)

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