From 9d3bcb0d617a1b3a1bf01394e2c30ac6369c4b22 Mon Sep 17 00:00:00 2001 From: PharmCat Date: Sun, 10 Jan 2021 15:48:16 +0300 Subject: [PATCH] update --- src/gmat.jl | 5 ++++- src/linearalgebra.jl | 18 ++++++++++++++++++ src/reml.jl | 20 ++++++++------------ src/reml_gh_deprecated.jl | 7 ++++++- src/rmat.jl | 4 ++++ src/statsbase.jl | 5 +++-- src/utils.jl | 15 ++++++++------- src/varstruct.jl | 10 ++++++++-- test/ar.jl | 11 +++++++++-- test/berds.jl | 4 ++-- test/lme4.jl | 23 ++++++++++++++++++----- test/norand.jl | 2 +- test/test.jl | 10 +++++++--- 13 files changed, 96 insertions(+), 38 deletions(-) diff --git a/src/gmat.jl b/src/gmat.jl index 3b8801dc..3ff2095d 100644 --- a/src/gmat.jl +++ b/src/gmat.jl @@ -33,6 +33,7 @@ function gmat_switch!(G, θ, covstr, i) G end ################################################################################ +#= function gmat_base_z!(mx, θ::Vector{T}, covstr) where T q = sum(length.(covstr.block[1])) for r = 1:length(covstr.random) @@ -44,8 +45,9 @@ function gmat_base_z!(mx, θ::Vector{T}, covstr) where T end mx end +=# ################################################################################ - +#= function gmat_base_z2!(mx, θ::Vector{T}, covstr, block) where T q = sum(length.(covstr.block[1])) for r = 1:length(covstr.random) @@ -62,6 +64,7 @@ function gmat_base_z2!(mx, θ::Vector{T}, covstr, block) where T end mx end +=# function gmat_base_z2!(mx, θ::Vector{T}, covstr, block, sblock) where T q = sum(length.(covstr.block[1])) for r = 1:length(covstr.random) diff --git a/src/linearalgebra.jl b/src/linearalgebra.jl index 08bb2650..187d9d98 100644 --- a/src/linearalgebra.jl +++ b/src/linearalgebra.jl @@ -2,6 +2,7 @@ """ A * B * A' + C """ +#= function mulαβαtc(A, B, C) q = size(B, 1) p = size(A, 1) @@ -23,9 +24,11 @@ function mulαβαtc(A, B, C) end Symmetric(mx) end +=# """ A * B * A' """ +#= function mulαβαt(A, B) q = size(B, 1) p = size(A, 1) @@ -92,6 +95,7 @@ function mulαβαtc3(A, B, C, X) mx[p+1:end, 1:p] = X' mx end +=# function mulαβαt3(A, B, X) q = size(B, 1) p = size(A, 1) @@ -118,6 +122,7 @@ end A' * B * A -> +θ A' * B * C -> +β """ +#= function mulθβinc!(θ, β, A::AbstractMatrix, B::AbstractMatrix, C::AbstractVector) q = size(B, 1) p = size(A, 2) @@ -138,10 +143,12 @@ function mulθβinc!(θ, β, A::AbstractMatrix, B::AbstractMatrix, C::AbstractVe end θ, β end +=# #------------------------------------------------------------------------------- """ A' * B * A -> + θ """ +#= function mulαtβαinc!(θ, A, B) q = size(B, 1) p = size(A, 2) @@ -156,6 +163,7 @@ function mulαtβαinc!(θ, A, B) end end end +=# """ A * B * A' -> + θ """ @@ -179,6 +187,7 @@ end """ A' * B * A -> θ """ +#= function mulαtβα!(θ, A, B) q = size(B, 1) p = size(A, 2) @@ -195,9 +204,11 @@ function mulαtβα!(θ, A, B) end θ end +=# """ A * B * A -> θ """ +#= function mulαβα!(θ, A, B) q = size(B, 1) p = size(A, 2) @@ -214,9 +225,11 @@ function mulαβα!(θ, A, B) end θ end +=# """ tr(A * B) """ +#= function trmulαβ(A, B) c = 0 @inbounds for n = 1:size(A,1), m = 1:size(B, 1) @@ -224,11 +237,14 @@ function trmulαβ(A, B) end c end +=# """ tr(H * A' * B * A) """ +#= function trmulhαtβα(H, A, B) end +=# """ (y - X * β)' * V * (y - X * β) """ @@ -252,6 +268,7 @@ end """ (y - X * β) """ +#= function mulr!(v::AbstractVector, y::AbstractVector, X::AbstractMatrix, β::AbstractVector) fill!(v, zero(eltype(v))) q = length(y) @@ -276,6 +293,7 @@ function mulr(y::AbstractVector, X::AbstractMatrix, β::AbstractVector) end return v end +=# """ A' * B -> + θ """ diff --git a/src/reml.jl b/src/reml.jl index 2ea6feef..e46c7fb9 100644 --- a/src/reml.jl +++ b/src/reml.jl @@ -3,6 +3,7 @@ """ -2 log Restricted Maximum Likelihood; """ +#= function reml_sweep(lmm, β, θ::Vector{T})::T where T <: Number n = length(lmm.data.yv) N = sum(length.(lmm.data.yv)) @@ -27,14 +28,16 @@ function reml_sweep(lmm, β, θ::Vector{T})::T where T <: Number θ₂ = logdet(θ2m) return -(θ₁ + θ₂ + θ₃ + c) end - +=# """ -2 log Restricted Maximum Likelihood; β calculation inside """ +#= function reml_sweep_β(lmm, f::Function, θ::Vector{T}) where T <: Number f(θ) reml_sweep_β(lmm, θ) end +=# """ -2 log Restricted Maximum Likelihood; β calculation inside """ @@ -69,9 +72,7 @@ function reml_sweep_β(lmm::LMM{T2}, @nospecialize θ::Vector{T}) where T <: Num catch θ₁ += Inf end - #println("logdet(V) $(i) : ", θ₁) - sweep!(Vp, 1:q) V⁻¹[i] = Symmetric(utriaply!(x -> -x, V)) #----------------------------------------------------------------------- @@ -93,7 +94,6 @@ function reml_sweep_β(lmm::LMM{T2}, @nospecialize θ::Vector{T}) where T <: Num @inbounds θ₃ += mulθ₃(view(lmm.data.yv, lmm.data.block[i]), view(lmm.data.xv, lmm.data.block[i],:), β, V⁻¹[i]) #println("θ₃ : ", θ₃) end - #logdetθ₂ = logdet(θ₂) try logdetθ₂ = logdet(θ₂) @@ -102,7 +102,7 @@ function reml_sweep_β(lmm::LMM{T2}, @nospecialize θ::Vector{T}) where T <: Num end return θ₁ + logdetθ₂ + θ₃ + c, β, θ₂, θ₃ end - +#= function reml_sweep_β2(lmm::LMM{T2}, @nospecialize θ::Vector{T}) where T <: Number where T2 <: Number n::Int = length(lmm.data.block) N::Int = length(lmm.data.yv) @@ -157,7 +157,7 @@ function reml_sweep_β2(lmm::LMM{T2}, @nospecialize θ::Vector{T}) where T <: Nu return θ₁ + logdetθ₂ + θ₃ + c, β, θ₂, θ₃ end - +=# function reml_sweep_β3(lmm::LMM{T2}, @nospecialize θ::Vector{T}) where T <: Number where T2 <: Number n::Int = length(lmm.data.block) N::Int = length(lmm.data.yv) @@ -187,11 +187,8 @@ function reml_sweep_β3(lmm::LMM{T2}, @nospecialize θ::Vector{T}) where T <: Nu Vx .= view(lmm.data.xv, lmm.data.block[i],:) gmat_base_z2!(V, θ, lmm.covstr, lmm.data.block[i], lmm.covstr.sblock[i]) rmat_basep_z2!(V, θ[lmm.covstr.tr[end]], lmm.covstr, lmm.data.block[i], lmm.covstr.sblock[i]) - #θ₁ += logdet(V) - θ₁ += logdet(V) - sweep!(Vp, 1:q) V⁻¹[i] = Symmetric(utriaply!(x -> -x, V)) #----------------------------------------------------------------------- @@ -206,9 +203,7 @@ function reml_sweep_β3(lmm::LMM{T2}, @nospecialize θ::Vector{T}) where T <: Nu #θ₃ += r' * V⁻¹[i] * r @inbounds θ₃ += mulθ₃(view(lmm.data.yv, lmm.data.block[i]), view(lmm.data.xv, lmm.data.block[i],:), β, V⁻¹[i]) end - logdetθ₂ = logdet(θ₂) - return θ₁ + logdetθ₂ + θ₃ + c, β, θ₂, θ₃ end #= @@ -258,6 +253,7 @@ function reml_sweep_β_ub(lmm::LMM{T2}, θ::Vector{T})::Tuple{T, Vector{T}, Matr end =# ################################################################################ +#= function reml_inv_β(lmm::LMM{T2}, θ::Vector{T})::Tuple{T, Vector{T}, Matrix{T}} where T <: Number where T2 <: Number n::Int = length(lmm.data.block) N::Int = length(lmm.data.yv) @@ -308,5 +304,5 @@ function reml_inv_β(lmm::LMM{T2}, θ::Vector{T})::Tuple{T, Vector{T}, Matrix{T} end return θ₁ + logdetθ₂ + θ₃ + c, β, θ₂ end - +=# ################################################################################ diff --git a/src/reml_gh_deprecated.jl b/src/reml_gh_deprecated.jl index 1b4991b6..0274ef34 100644 --- a/src/reml_gh_deprecated.jl +++ b/src/reml_gh_deprecated.jl @@ -2,6 +2,7 @@ """ 2 log Restricted Maximum Likelihood gradient vector """ +#= function reml_grad(yv, Zv, p, Xv, θvec, β) n = length(yv) G = gmat(θvec[3:5]) @@ -40,10 +41,11 @@ function reml_grad(yv, Zv, p, Xv, θvec, β) end return - (θ1 .+ θ2 .+ θ3) end - +=# """ 2 log Restricted Maximum Likelihood hessian matrix """ +#= function reml_hessian(yv, Zv, p, Xv, θvec, β) n = length(yv) G = gmat(θvec[3:5]) @@ -91,7 +93,9 @@ function reml_hessian(yv, Zv, p, Xv, θvec, β) end return - (θ1 .+ θ2 .+ θ3) end +=# ################################################################################ +#= function reml_grad2(yv, Zv, p, Xv, θvec, β) n = length(yv) for i = 1:n @@ -102,3 +106,4 @@ function reml_grad3(yv, Zv, p, Xv, θvec, β) n = length(yv) covmat_grad2(vmatvec, Zv, θvec) end +=# diff --git a/src/rmat.jl b/src/rmat.jl index 38a2c358..f2f38e13 100644 --- a/src/rmat.jl +++ b/src/rmat.jl @@ -21,6 +21,7 @@ function rmat_basep!(mx, θ::AbstractVector{T}, zrv, covstr::CovStructure{T2}) w end end ################################################################################ +#= function rmat_basep_z!(mx, θ::AbstractVector{T}, zrv, covstr) where T for i = 1:length(covstr.block[end]) if covstr.repeated.covtype.s == :SI @@ -41,7 +42,9 @@ function rmat_basep_z!(mx, θ::AbstractVector{T}, zrv, covstr) where T end mx end +=# ################################################################################ +#= function rmat_basep_z2!(mx, θ::AbstractVector{T}, covstr, block) where T subjblock = view(covstr.subjz[end], block, :) zblock = view(covstr.rz, block, :) @@ -53,6 +56,7 @@ function rmat_basep_z2!(mx, θ::AbstractVector{T}, covstr, block) where T end mx end +=# function rmat_basep_z2!(mx, θ::AbstractVector{T}, covstr, block, sblock) where T zblock = view(covstr.rz, block, :) for i = 1:length(sblock[end]) diff --git a/src/statsbase.jl b/src/statsbase.jl index 1da5c5fb..f1055abc 100644 --- a/src/statsbase.jl +++ b/src/statsbase.jl @@ -1,6 +1,5 @@ - - +#= StatsBase.coef(model::LMM) = error("coef is not defined for $(typeof(model)).") StatsBase.coefnames(model::LMM) = error("coefnames is not defined for $(typeof(model)).") @@ -122,3 +121,5 @@ StatsBase.predict(model::LMM) = error("predict is not defined for $(typeof(model StatsBase.predict!(model::LMM) = error("predict! is not defined for $(typeof(model)).") StatsBase.dof_residual(model::LMM) = error("dof_residual is not defined for $(typeof(model)).") + +=# diff --git a/src/utils.jl b/src/utils.jl index d1c1cf53..d20338df 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,6 +1,7 @@ """ Make X, Z matrices and vector y for each subject; """ +#= function subjblocks(df, sbj::Symbol, x::Matrix{T}, z::Matrix{T}, y::Vector{T}, rz::Matrix{T}) where T u = unique(df[!, sbj]) xa = Vector{Matrix{T}}(undef, length(u)) @@ -44,6 +45,7 @@ function subjblocks(df, sbj) r[1] = collect(1:size(df, 1)) r end +=# """ Intersect dataframe. """ @@ -73,7 +75,7 @@ function intersectdf(df, s)::Vector end res end - +#= function intersectsubj(covstr) a = Vector{Vector{Symbol}}(undef, length(covstr.random)+1) eq = true @@ -89,6 +91,7 @@ function intersectsubj(covstr) end intersect(a...), eq end +=# function intersectsubj(random, repeated) a = Vector{Vector{Symbol}}(undef, length(random)+1) eq = true @@ -104,12 +107,12 @@ function intersectsubj(random, repeated) end intersect(a...), eq end - +#= function diffsubj!(a, subj) push!(a, subj) symdiff(a...) end - +=# """ Variance estimate via OLS and QR decomposition. """ @@ -208,7 +211,7 @@ function varlinkrvecapply2!(v, p; varlinkf = :exp, rholinkf = :sigm) end ################################################################################ - +#= function vmatr(lmm, i) θ = lmm.result.theta G = gmat_base(θ, lmm.covstr) @@ -225,7 +228,7 @@ function gmatr(lmm, i) θ = lmm.result.theta gmat_base(θ, lmm.covstr) end - +=# ################################################################################ function m2logreml(lmm) @@ -236,8 +239,6 @@ function logreml(lmm) end ################################################################################ - - function optim_callback(os) false end diff --git a/src/varstruct.jl b/src/varstruct.jl index e494526f..b2617b83 100644 --- a/src/varstruct.jl +++ b/src/varstruct.jl @@ -19,13 +19,14 @@ end function ffxpone(x::T)::T where T x + one(T) end +#= function ffxmone(x::T)::T where T x - one(T) end function ff2xmone(x::T)::T where T 2x - one(T) end - +=# ################################################################################ # COVARIANCE TYPE ################################################################################ @@ -263,6 +264,7 @@ function fillur!(ur, i, v) end end ################################################################################ +#= function schemalength(s) if isa(s, Tuple) return length(s) @@ -270,6 +272,7 @@ function schemalength(s) return 1 end end +=# # function subjmatrix!(subj, data, subjz, i) @@ -326,10 +329,13 @@ end """ Return variance-covariance matrix V """ +#= function vmat(G, R, Z)::AbstractMatrix return mulαβαtc(Z, G, R) end +=# ################################################################################ +#= function vmatvec(G, Z, θ) v = Vector{Matrix{eltype(G)}}(undef, length(Z)) for i = 1:length(Z) @@ -338,7 +344,7 @@ function vmatvec(G, Z, θ) #reduce(vcat, v) v end - +=# ################################################################################ function Base.show(io::IO, e::VarEffect) println(io, "Effect") diff --git a/test/ar.jl b/test/ar.jl index 116575ed..f97446b9 100644 --- a/test/ar.jl +++ b/test/ar.jl @@ -9,7 +9,7 @@ sort!(df, :Day) transform!(df, :Time => categorical, renamecols=false) transform!(df, :Day => categorical, renamecols=false) -@testset " RepeatedPulse.csv" begin +@testset " AR RepeatedPulse.csv " begin lmm = Metida.LMM(@formula(Pulse~1), df; random = Metida.VarEffect(Metida.@covstr(Time), Metida.SI), repeated = Metida.VarEffect(Metida.@covstr(Day), Metida.AR), @@ -24,6 +24,13 @@ transform!(df, :Day => categorical, renamecols=false) ) Metida.fit!(lmm) @test lmm.result.reml ≈ 471.85107712169827 atol=1E-6 + + lmm = Metida.LMM(@formula(Pulse~1), df; + random = Metida.VarEffect(Metida.@covstr(Day), Metida.AR), + subject = :Time + ) + Metida.fit!(lmm) + @test lmm.result.reml ≈ 453.3395560121246 atol=1E-6 end #= @@ -37,7 +44,7 @@ df.Time2 = copy(df.Time) transform!(df, :Time2 => categorical, renamecols=false) transform!(df, :Chick => categorical, renamecols=false) transform!(df, :Diet => categorical, renamecols=false) -@testset " ChickWeight.csv" begin +@testset " ARH ChickWeight.csv " begin lmm = Metida.LMM(@formula(weight~1 + Diet & Time), df; random = Metida.VarEffect(Metida.@covstr(1), Metida.SI), repeated = Metida.VarEffect(Metida.@covstr(Diet), Metida.ARH), diff --git a/test/berds.jl b/test/berds.jl index 7bffe070..ac4d33dd 100644 --- a/test/berds.jl +++ b/test/berds.jl @@ -36,7 +36,7 @@ for i = 1:30 df.lnpk = log.(df.PK) - @testset " Test $(i)" begin + @testset " RDS Test $(i) " begin atol=1E-6 lmm = Metida.LMM(@formula(lnpk~sequence+period+treatment), df; @@ -46,7 +46,7 @@ for i = 1:30 Metida.fit!(lmm) @test lmm.result.reml ≈ remlsb[i] atol=atol - + if i == 13 || i == 15 atol = 1E-4 end lmm = Metida.LMM(@formula(lnpk~sequence+period+treatment), df; random = Metida.VarEffect(Metida.@covstr(treatment), Metida.CSH), diff --git a/test/lme4.jl b/test/lme4.jl index 60f3551c..dc554fbb 100644 --- a/test/lme4.jl +++ b/test/lme4.jl @@ -4,12 +4,11 @@ df = CSV.File(path*"/csv/lme4/sleepstudy.csv") |> DataFrame transform!(df, :Subject => categorical, renamecols=false) transform!(df, :Days => categorical, renamecols=false) - #= SPSS REML 1729.492560 =# -@testset " sleepstudy.csv" begin +@testset " sleepstudy.csv " begin lmm = Metida.LMM(@formula(Reaction~Days), df; random = Metida.VarEffect(Metida.SI), subject = :Subject @@ -18,13 +17,27 @@ REML 1729.492560 @test lmm.result.reml ≈ 1729.4925602367025 atol=1E-6 end +@testset " CS sleepstudy.csv " begin + lmm = Metida.LMM(@formula(Reaction~1), df; + random = Metida.VarEffect(Metida.@covstr(Days), Metida.CS), + subject = :Subject + ) + Metida.fit!(lmm) + @test lmm.result.reml ≈ 1904.3265170722132 atol=1E-6 + + lmm = Metida.LMM(@formula(Reaction~1 + Days), df; + repeated = Metida.VarEffect(Metida.@covstr(1), Metida.CS, subj = :Subject) + ) + Metida.fit!(lmm) + @test lmm.result.reml ≈ 1729.4925602367027 atol=1E-6 +end + df = CSV.File(path*"/csv/lme4/Penicillin.csv") |> DataFrame df.diameter = float.(df.diameter) transform!(df, :plate => categorical, renamecols=false) transform!(df, :sample => categorical, renamecols=false) - -@testset " Penicillin.csv" begin +@testset " Penicillin.csv " begin lmm = Metida.LMM(@formula(diameter~1), df; random = [Metida.VarEffect(Metida.SI, subj = :plate), Metida.VarEffect(Metida.SI, subj = :sample)] ) @@ -38,7 +51,7 @@ transform!(df, :batch => categorical, renamecols=false) transform!(df, :sample => categorical, renamecols=false) transform!(df, :cask=> categorical, renamecols=false) -@testset " Pastes.csv" begin +@testset " Pastes.csv " begin lmm = Metida.LMM(@formula(strength~1), df; random = [Metida.VarEffect(Metida.SI, subj = :batch), Metida.VarEffect(Metida.SI, subj = [:batch, :cask])] ) diff --git a/test/norand.jl b/test/norand.jl index bbe14edb..bc0c64a9 100644 --- a/test/norand.jl +++ b/test/norand.jl @@ -1,4 +1,4 @@ -@testset " No random " begin +@testset " No random " begin lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0; repeated = Metida.VarEffect(Metida.@covstr(formulation), Metida.DIAG), subject = :subject diff --git a/test/test.jl b/test/test.jl index 5f00af35..40ee23a6 100644 --- a/test/test.jl +++ b/test/test.jl @@ -5,15 +5,19 @@ using Test, CSV, DataFrames, StatsModels path = dirname(@__FILE__) include("testdata.jl") -@testset " Model: DIAG + nothing " begin +@testset " Basic test, DIAG / SI " begin lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0; random = Metida.VarEffect(Metida.@covstr(formulation), Metida.DIAG), ) Metida.fit!(lmm) io = IOBuffer(); Base.show(io, lmm) - - @test lmm.result.reml ≈ 25.129480634331063 atol=1E-8 #need check + Base.show(io, lmm.data) + Base.show(io, lmm.result) + Base.show(io, lmm.covstr) + @test Metida.logreml(lmm) ≈ -12.564740317165533 atol=1E-6 + @test Metida.m2logreml(lmm) ≈ 25.129480634331067 atol=1E-6 + @test lmm.result.reml ≈ 25.129480634331063 atol=1E-6 #need check end @testset " Model: DIAG/subject + nothing " begin lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;