update

Project.toml

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

change.log

@@ -1 +1,5 @@
+- 0.1.9
+  Experimental part completed
+
 - 0.1.0
+  Initial

src/fit.jl

@@ -37,7 +37,6 @@ function fit!(lmm::LMM{T}; verbose::Symbol = :auto, varlinkf = :exp, rholinkf =
     θ  = zeros(T, lmm.covstr.tl)
     #θ                      .= initθ
     θ                      .= initθ ./2
-    #θ[lmm.covstr.tr[end]]  .= 0.01
     #θ .= initθ / (length(lmm.covstr.random) + 1)
     for i = 1:length(θ)
         if lmm.covstr.ct[i] == :rho θ[i] = 0.0 end
@@ -88,19 +87,21 @@ function fit!(lmm::LMM{T}; verbose::Symbol = :auto, varlinkf = :exp, rholinkf =
                 hsvd.S[i] = 0
             end
         end
+        lmm.result.hsvds = copy(hsvd.S)
         rhsvd = hsvd.U * Diagonal(hsvd.S) * hsvd.Vt
+        theta = copy(lmm.result.theta)
         for i = 1:length(lmm.result.theta)
             if rhsvd[i,i] < 1E-10
                 if lmm.covstr.ct[i] == :var
-                    lmm.result.theta[i] = 0
+                    theta[i] = 0
                     push!(lmm.warn, "Variation parameter ($(i)) set to zero.")
                 elseif lmm.covstr.ct[i] == :rho
                     push!(lmm.warn, "Rho SVD value ($(i)) is less than 1e-10.")
                 end
             end
         end
         #-2 LogREML, β, iC
-        lmm.result.reml, lmm.result.beta, iC = optfunc(lmm, lmm.result.theta)
+        lmm.result.reml, lmm.result.beta, iC = optfunc(lmm, theta)
         #Variance-vovariance matrix of β
         lmm.result.c            = pinv(iC)
         #SE

src/lmm.jl

@@ -41,13 +41,13 @@ struct LMM{T} <: MetidaModel
         warn = Vector{String}(undef, 0)
         mf   = ModelFrame(model, data; contrasts = contrasts)
         mm   = ModelMatrix(mf)
-        if random === nothing
-            random = VarEffect()
-            #random = RZero()
-        end
         if repeated === nothing
             repeated = VarEffect()
         end
+        if random === nothing
+            random = VarEffect(Metida.@covstr(0), Metida.RZero(), subj = repeated.subj)
+            #random = RZero()
+        end
         if !isa(random, Vector) random = [random] end
         #blocks
         intsub, eq = intersectsubj(random, repeated)
@@ -79,6 +79,10 @@ function Base.show(io::IO, lmm::LMM)
     println(io, "Linear Mixed Model: ", lmm.model)
     for i = 1:length(lmm.covstr.random)
         println(io, "Random $i: ")
+        if lmm.covstr.random[i].covtype.s == :ZERO
+            println(io, "   No")
+            continue
+        end
         println(io, "   Model: $(lmm.covstr.random[i].model === nothing ? "nothing" : lmm.covstr.random[i].model)")
         println(io, "   Type: $(lmm.covstr.random[i].covtype.s) ($(lmm.covstr.t[i]))")
         #println(io, "   Coefnames: $(coefnames(lmm.covstr.schema[i]))")

src/modelresult.jl

@@ -3,31 +3,34 @@
 mutable struct ModelResult
     fit::Bool
     optim
-    theta
-    reml
-    beta
-    h
-    c
-    se
+    theta::Union{Vector, Nothing}
+    reml::Union{Float64, Nothing}
+    beta::Union{Vector, Nothing}
+    h::Union{Matrix, Nothing}
+    c::Union{Matrix, Nothing}
+    se::Union{Vector, Nothing}
+    hsvds::Union{Vector, Nothing}
     function ModelResult(
     optim,
     theta,
     reml,
     beta,
     h,
     c,
-    se)
+    se,
+    hsvds)
         new(true,
         optim,
         theta,
         reml,
         beta,
         h,
         c,
-        se)
+        se,
+        hsvds)
     end
     function ModelResult()
-        new(false, nothing, nothing, nothing, nothing, nothing, nothing, nothing)
+        new(false, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing)
     end
 end
 

src/varstruct.jl

@@ -139,7 +139,7 @@ struct CovStructure{T} <: AbstractCovarianceStructure
     block::Vector{Vector{Vector{UInt32}}}
     # Z matrix
     z::Matrix{T}
-    subjz::Vector{BitArray{2}}
+    #subjz::Vector{BitArray{2}}
     # Blocks for each blocking subject, each effect, each effect subject
     sblock::Vector{Vector{Vector{Vector{UInt32}}}}
     #unit range z column range for each random effect
@@ -169,18 +169,11 @@ struct CovStructure{T} <: AbstractCovarianceStructure
         z       = Matrix{Float64}(undef, size(data, 1), 0)
         subjz   = Vector{BitArray{2}}(undef, alleffl)
         sblock  = Vector{Vector{Vector{Vector{UInt32}}}}(undef, length(blocks))
-        zrndur  = Vector{UnitRange{Int}}(undef, alleffl)
+        zrndur  = Vector{UnitRange{Int}}(undef, alleffl - 1)
         rz      = Matrix{Float64}(undef, size(data, 1), 0)
         #
         # RANDOM EFFECTS
         for i = 1:length(random)
-            if random[i].covtype.s == :ZERO
-                q[i] = 0
-                t[i] = 0
-                tr[i] = 0:0
-                zrndur[i] = 0:0
-                continue
-            end
             if length(random[i].coding) == 0 && random[i].fulldummy
                 fill_coding_dict!(random[i].model, random[i].coding, data)
             end
@@ -243,7 +236,7 @@ struct CovStructure{T} <: AbstractCovarianceStructure
             end
         end
         view(rcnames, tr[end]) .= rcoefnames(schema[end], t[end], Val{repeated.covtype.s}())
-        if random[1].covtype.s != :ZERO
+
             for i = 1:length(blocks)
                 sblock[i] = Vector{Vector{Vector{UInt32}}}(undef, alleffl)
                 for s = 1:alleffl
@@ -253,17 +246,21 @@ struct CovStructure{T} <: AbstractCovarianceStructure
                     end
                 end
             end
-        end
+
         #
-        new{eltype(z)}(random, repeated, schema, rcnames, block, z, subjz, sblock, zrndur, rz, q, t, tr, tl, ct)
+        new{eltype(z)}(random, repeated, schema, rcnames, block, z, sblock, zrndur, rz, q, t, tr, tl, ct)
     end
 end
 ################################################################################
 function fillur!(ur, i, v)
     if i > 1
         ur[i]   = UnitRange(sum(v[1:i-1]) + 1, sum(v[1:i-1]) + v[i])
     else
-        ur[1]   = UnitRange(1, v[1])
+        if v[1] > 0
+            ur[1]   = UnitRange(1, v[1])
+        else
+            ur[1]   = UnitRange(0, 0)
+        end
     end
 end
 ################################################################################
@@ -352,3 +349,20 @@ function Base.show(io::IO, e::VarEffect)
     println(io, "FullDummy: ", e.fulldummy)
     println(io, "Subject:", e.subj)
 end
+
+
+function Base.show(io::IO, cs::CovStructure)
+    println(io, "Covariance Structure:")
+    for i = 1:length(cs.random)
+        println(io, "Random $(i):", cs.random[i])
+    end
+    println(io, "Repeated: ", cs.repeated)
+    println(io, "Random effect range in complex Z: ", cs.zrndur)
+    println(io, "Size of Z: ", cs.q)
+    println(io, "Parameter number for each effect: ", cs.t)
+    println(io, "Theta length:", cs.tl)
+end
+
+function Base.show(io::IO, ct::CovarianceType)
+    println(io, "Covariance Type: $(ct.s)")
+end

test/norand.jl

@@ -0,0 +1,14 @@
+@testset "  No random                            " begin
+    lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
+    repeated = Metida.VarEffect(Metida.@covstr(formulation), Metida.DIAG),
+    subject = :subject
+    )
+    Metida.fit!(lmm)
+
+    lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
+    repeated = Metida.VarEffect(Metida.@covstr(formulation), Metida.DIAG, subj = :subject)
+    )
+    Metida.fit!(lmm)
+
+    @test lmm.result.reml ≈ 25.000777869122338 atol=1E-8
+end

test/test.jl

@@ -156,3 +156,7 @@ end
     Metida.fit!(lmm)
     @test lmm.result.reml ≈ 25.00077786912234 atol=1E-4 #need check
 end
+
+include("ar.jl")
+include("lme4.jl")
+include("norand.jl")