Merge pull request #41 from PharmCat/master-test

fix Int, minor changes

src/fit.jl

@@ -132,7 +132,7 @@ function fit!(lmm::LMM{T}; kwargs...) where T
             error("init length $(length(init)) != θ length $(length(θ))")
         end
     else
-        initθ = sqrt(initvar(lmm.data.yv, lmm.mm.m)[1])/(length(lmm.covstr.random)+1)
+        initθ = sqrt(initvar(lmm.data.yv, lmm.data.xv)[1])/(length(lmm.covstr.random)+1)
         for i = 1:length(θ)
             if lmm.covstr.ct[i] == :var
                 θ[i] = initθ

src/lmm.jl

@@ -26,7 +26,7 @@ struct LMM{T<:AbstractFloat} <: MetidaModel
     model::FormulaTerm
     mf::ModelFrame
     mm::ModelMatrix
-    covstr::CovStructure{T}
+    covstr::CovStructure
     data::LMMData{T}
     dv::LMMDataViews{T}
     nfixed::Int
@@ -38,15 +38,15 @@ struct LMM{T<:AbstractFloat} <: MetidaModel
     function LMM(model::FormulaTerm,
         mf::ModelFrame,
         mm::ModelMatrix,
-        covstr::CovStructure{T},
+        covstr::CovStructure,
         data::LMMData{T},
         dv::LMMDataViews{T},
         nfixed::Int,
         rankx::Int,
         result::ModelResult,
         maxvcbl::Int,
         log::Vector{LMMLogMsg}) where T
-        new{eltype(mm.m)}(model, mf, mm, covstr, data, dv, nfixed, rankx, result, maxvcbl, log)
+        new{T}(model, mf, mm, covstr, data, dv, nfixed, rankx, result, maxvcbl, log)
     end
     function LMM(model, data; contrasts=Dict{Symbol,Any}(),  random::Union{Nothing, VarEffect, Vector{VarEffect}} = nothing, repeated::Union{Nothing, VarEffect} = nothing)
         #need check responce - Float
@@ -86,11 +86,13 @@ struct LMM{T<:AbstractFloat} <: MetidaModel
         end
         lmmdata = LMMData(modelmatrix(mf), response(mf))
         covstr = CovStructure(random, repeated, data)
+        coefn = size(lmmdata.xv, 2)
         rankx =  rank(lmmdata.xv)
-        if rankx != size(lmmdata.xv, 2)
+        if rankx != coefn
+            @warn "Fixed-effect matrix not full-rank!"
             lmmlog!(lmmlog, 1, LMMLogMsg(:WARN, "Fixed-effect matrix not full-rank!"))
         end
-        mres = ModelResult(false, nothing, fill(NaN, covstr.tl), NaN, fill(NaN, rankx), nothing, fill(NaN, rankx, rankx), fill(NaN, rankx), nothing, false)
+        mres = ModelResult(false, nothing, fill(NaN, covstr.tl), NaN, fill(NaN, coefn), nothing, fill(NaN, coefn, coefn), fill(NaN, coefn), nothing, false)
         LMM(model, mf, mm, covstr, lmmdata, LMMDataViews(lmmdata.xv, lmmdata.yv, covstr.vcovblock), nfixed, rankx, mres, findmax(length, covstr.vcovblock)[1], lmmlog)
     end
     function LMM(f::LMMformula, data; contrasts=Dict{Symbol,Any}(), kwargs...)

src/lmmdata.jl

@@ -5,9 +5,12 @@ struct LMMData{T}
     xv::Matrix{Float64}
     # Responce vector
     yv::Vector{T}
-    function LMMData(xa::Matrix{Float64}, ya::Vector{T}) where T
+    function LMMData(xa::AbstractMatrix{Float64}, ya::AbstractVector{T}) where T
         new{T}(xa, ya)
     end
+    function LMMData(xa::AbstractMatrix{Float64}, ya::AbstractVector{Int})
+        LMMData(xa, float.(ya))
+    end
 end
 
 struct LMMDataViews{T} <: AbstractLMMDataBlocks

src/reml.jl

@@ -64,7 +64,7 @@ function reml_sweep_β(lmm, data, θ::Vector{T}; maxthreads::Int = 16) where T #
             @inbounds for j ∈ 1:d + (t ≤ r)
                 i    =  offset + j
                 q    = length(lmm.covstr.vcovblock[i])
-                qswm = q + lmm.rankx
+                qswm = q + p
                 Vp   = zeros(T, qswm, qswm)
                 V    = view(Vp, 1:q, 1:q)
                 Vx   = view(Vp, 1:q, q+1:qswm)
@@ -103,19 +103,43 @@ function reml_sweep_β(lmm, data, θ::Vector{T}; maxthreads::Int = 16) where T #
                 θ₃ += mulθ₃(data.yv[i], data.xv[i], β, V⁻¹[i])
             end
         else
+            #=
+            qrd  = qr(θs₂)
+            vec  = collect(1:length(βm))
+            dr   = diag(qrd.R)
+            tval = mean(dr)*sqrt(eps())
+            inds = Int[]
+            for i = 1:length(βm)
+                if dr[i] > tval 
+                    push!(inds, i) 
+                else
+                    θ₂[:, i] .= zero(T)
+                    θ₂[i, :] .= zero(T)
+                    βm[i]     = zero(T)
+                end
+            end
+            rθs₂ = θs₂[inds, inds]
+            mul!(view(β, inds), inv(rθs₂), view(βm, inds))
+            logdetθ₂ = logdet(rθs₂)
+            =#
             β       .= NaN
             return   Inf, β, Inf, Inf, false
         end
+        # θ₃
+        #@inbounds @simd for i = 1:n
+        #    θ₃ += mulθ₃(data.yv[i], data.xv[i], β, V⁻¹[i])
+        #end
     return   θ₁ + logdetθ₂ + θ₃ + c, β, θs₂, θ₃, noerror #REML, β, iC, θ₃, errors
 end
 # Using BLAS, LAPACK - non ForwardDiff, used by MetidaNLopt
 function reml_sweep_β_nlopt(lmm, data, θ::Vector{T}; maxthreads::Int = 16) where T
     n             = length(lmm.covstr.vcovblock)
     N             = length(lmm.data.yv)
     c             = (N - lmm.rankx)*log(2π)
+    p             = size(lmm.data.xv, 2)
     #---------------------------------------------------------------------------
     θ₁            = zero(T)
-    θ₂            = zeros(T, lmm.rankx, lmm.rankx)
+    θ₂            = zeros(T, p, p)
     θ₃            = zero(T)
     A             = Vector{Matrix{T}}(undef, n)
     logdetθ₂      = zero(T)
@@ -130,8 +154,8 @@ function reml_sweep_β_nlopt(lmm, data, θ::Vector{T}; maxthreads::Int = 16) whe
         d, r = divrem(n, ncore)
         Base.Threads.@threads for t = 1:ncore
             offset   = min(t-1, r) + (t-1)*d
-            accθ₂[t] = zeros(T, lmm.rankx, lmm.rankx)
-            accβm[t] = zeros(T, lmm.rankx)
+            accθ₂[t] = zeros(T, p, p)
+            accβm[t] = zeros(T, p)
             @inbounds for j ∈ 1:d+(t ≤ r)
                 i =  offset + j
                 q    = length(lmm.covstr.vcovblock[i])
@@ -188,6 +212,7 @@ end
 ################################################################################
 function core_sweep_β(lmm, data, θ::Vector{T}, β, n; maxthreads::Int = 16) where T
     ncore     = min(num_cores(), n, maxthreads)
+    p         = size(lmm.data.xv, 2)
     accθ₁     = zeros(T, ncore)
     accθ₂     = Vector{Matrix{T}}(undef, ncore)
     accθ₃     = zeros(T, ncore)
@@ -197,11 +222,11 @@ function core_sweep_β(lmm, data, θ::Vector{T}, β, n; maxthreads::Int = 16) wh
     d, r = divrem(n, ncore)
     Base.Threads.@threads for t = 1:ncore
         offset = min(t-1, r) + (t-1)*d
-        accθ₂[t] = zeros(T, lmm.rankx, lmm.rankx)
+        accθ₂[t] = zeros(T, p, p)
         @inbounds for j ∈ 1:d+(t ≤ r)
             i =  offset + j
             q    = length(lmm.covstr.vcovblock[i])
-            qswm = q + lmm.rankx
+            qswm = q + p
             Vp   = zeros(T, qswm, qswm)
             V    = view(Vp, 1:q, 1:q)
             Vx   = view(Vp, 1:q, q+1:qswm)

src/utils.jl

@@ -51,7 +51,7 @@ L-contrast matrix for `i` fixed effect.
 """
 function lcontrast(lmm::LMM, i::Int)
     n = length(lmm.mf.f.rhs.terms)
-    p = size(lmm.mm.m, 2)
+    p = size(lmm.data.xv, 2)
     if i > n || n < 1 error("Factor number out of range 1-$(n)") end
     inds = findall(x -> x==i, lmm.mm.assign)
     if typeof(lmm.mf.f.rhs.terms[i]) <: CategoricalTerm

test/test.jl

@@ -169,6 +169,13 @@ include("testdata.jl")
     @test tt.ndf[2] ≈ 3.0 atol=1E-5
     @test tt.df[2] ≈ 3.39086 atol=1E-5
     @test tt.pval[2] ≈ 0.900636 atol=1E-5
+
+    # Int dependent variable, function Term in random part
+    df0.varint = Int.(ceil.(df0.var2))
+    lmmint = Metida.fit(Metida.LMM, Metida.@lmmformula(varint~formulation,
+    random = 1+var^2|subject:Metida.SI), df0)
+    Metida.fit!(lmmint)
+    @test Metida.m2logreml(lmmint) ≈ 84.23373276096902 atol=1E-6
 end
 ################################################################################
 #                                  df0