From 97268930d0fbd62031d621de3910acf7b38a9575 Mon Sep 17 00:00:00 2001 From: PharmCat Date: Mon, 4 Dec 2023 21:32:46 +0300 Subject: [PATCH] fix Int, minor changes --- src/fit.jl | 2 +- src/lmm.jl | 12 +++++++----- src/lmmdata.jl | 5 ++++- src/reml.jl | 37 +++++++++++++++++++++++++++++++------ src/utils.jl | 2 +- test/test.jl | 7 +++++++ 6 files changed, 51 insertions(+), 14 deletions(-) diff --git a/src/fit.jl b/src/fit.jl index 83e020a0..baa1118b 100644 --- a/src/fit.jl +++ b/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θ diff --git a/src/lmm.jl b/src/lmm.jl index af444c3f..2181f89a 100644 --- a/src/lmm.jl +++ b/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,7 +38,7 @@ 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, @@ -46,7 +46,7 @@ struct LMM{T<:AbstractFloat} <: MetidaModel 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...) diff --git a/src/lmmdata.jl b/src/lmmdata.jl index 7f346834..b83aca34 100644 --- a/src/lmmdata.jl +++ b/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 diff --git a/src/reml.jl b/src/reml.jl index fc7e74c3..445add8d 100644 --- a/src/reml.jl +++ b/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,9 +103,32 @@ 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 @@ -113,9 +136,10 @@ function reml_sweep_β_nlopt(lmm, data, θ::Vector{T}; maxthreads::Int = 16) whe 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) diff --git a/src/utils.jl b/src/utils.jl index 3f66d54f..c76694a6 100644 --- a/src/utils.jl +++ b/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 diff --git a/test/test.jl b/test/test.jl index 303ad020..a609fdbf 100644 --- a/test/test.jl +++ b/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