From aa2ea542c50670e143988e733bc98b4b14e9291c Mon Sep 17 00:00:00 2001 From: PharmCat Date: Thu, 7 Jan 2021 04:12:22 +0300 Subject: [PATCH] update --- Project.toml | 2 +- change.log | 4 ++++ src/fit.jl | 7 ++++--- src/lmm.jl | 12 ++++++++---- src/modelresult.jl | 21 ++++++++++++--------- src/varstruct.jl | 40 +++++++++++++++++++++++++++------------- test/norand.jl | 14 ++++++++++++++ test/test.jl | 4 ++++ 8 files changed, 74 insertions(+), 30 deletions(-) create mode 100644 test/norand.jl diff --git a/Project.toml b/Project.toml index f281f912..edbb987d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Metida" uuid = "a1dec852-9fe5-11e9-361f-8d9fde67cfa2" authors = ["Vladimir Arnautov "] -version = "0.1.8" +version = "0.1.9" [deps] CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597" diff --git a/change.log b/change.log index 0559bb3e..d3ed50b6 100644 --- a/change.log +++ b/change.log @@ -1 +1,5 @@ +- 0.1.9 + Experimental part completed + - 0.1.0 + Initial diff --git a/src/fit.jl b/src/fit.jl index 067270e7..c4f664bf 100644 --- a/src/fit.jl +++ b/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,11 +87,13 @@ 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.") @@ -100,7 +101,7 @@ function fit!(lmm::LMM{T}; verbose::Symbol = :auto, varlinkf = :exp, rholinkf = 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 diff --git a/src/lmm.jl b/src/lmm.jl index 42a5661f..0854abc7 100644 --- a/src/lmm.jl +++ b/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]))") diff --git a/src/modelresult.jl b/src/modelresult.jl index bff3d558..4f544b37 100644 --- a/src/modelresult.jl +++ b/src/modelresult.jl @@ -3,12 +3,13 @@ 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, @@ -16,7 +17,8 @@ mutable struct ModelResult beta, h, c, - se) + se, + hsvds) new(true, optim, theta, @@ -24,10 +26,11 @@ mutable struct ModelResult 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 diff --git a/src/varstruct.jl b/src/varstruct.jl index e807abd6..da207527 100644 --- a/src/varstruct.jl +++ b/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,9 +246,9 @@ 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 ################################################################################ @@ -263,7 +256,11 @@ 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 diff --git a/test/norand.jl b/test/norand.jl new file mode 100644 index 00000000..bbe14edb --- /dev/null +++ b/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 diff --git a/test/test.jl b/test/test.jl index fadae369..2e48347d 100644 --- a/test/test.jl +++ b/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")