diff --git a/Project.toml b/Project.toml index d1ba0a2d..013d4ad7 100644 --- a/Project.toml +++ b/Project.toml @@ -3,10 +3,11 @@ uuid = "a1dec852-9fe5-11e9-361f-8d9fde67cfa2" keywords = ["lenearmodel", "mixedmodel"] desc = "Mixed-effects models with flexible covariance structure." authors = ["Vladimir Arnautov "] -version = "0.14.5" +version = "0.14.6" [deps] DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" @@ -19,13 +20,14 @@ StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d" [compat] DiffResults = "1" +Distributions = "0.20, 0.21, 0.22, 0.23, 0.24, 0.25" ForwardDiff = "0.10" LineSearches = "7" -MetidaBase = "0.10.1, 0.10.2" +MetidaBase = "0.11" Optim = "1" ProgressMeter = "1" StatsBase = "0.29, 0.30, 0.31, 0.32, 0.33" -StatsModels = "0.6" +StatsModels = "0.7" julia = "1" [extras] diff --git a/docs/Project.toml b/docs/Project.toml index 82b90aae..a31cd3ca 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -3,7 +3,6 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" -MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" @@ -17,8 +16,7 @@ Plots = "≥1" StatsPlots = "≥0.14" CSV = "≥0.8" DataFrames = "≥1" -MixedModels = "3.1.5" -PrettyTables = "0.10, 0.11, 1, 2" +PrettyTables = "1, 2" StatsBase = "≥0.33" Distributions = "≥0.25" CategoricalArrays = "≥0.9, 0.10" diff --git a/docs/src/examples.md b/docs/src/examples.md index 5e9c075c..21087aff 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -1,7 +1,13 @@ ### Example 1 - Continuous and categorical predictors ```@example lmmexample -using Metida, CSV, DataFrames, MixedModels, CategoricalArrays; +using Metida, CSV, DataFrames, CategoricalArrays; + +import Pkg +Pkg.activate("MixedModels") +Pkg.add(name="Example", version="3.1.5") +using MixedModels + rds = CSV.File(joinpath(dirname(pathof(Metida)), "..", "test", "csv", "1fptime.csv"); types = [String, String, Float64, Float64]) |> DataFrame @@ -156,3 +162,20 @@ MODEL var = sequence period formulation/ DDFM=SATTERTH s; RANDOM subject/TYPE=VC G V; RUN; ``` + +### Example 5 - Working with Effects.jl + +``` +using Effects, StatsModels + +lmm = LMM(@formula(var ~ sequence + period + formulation), df0; + random = VarEffect(@covstr(subject|1), SI) + ) +fit!(lmm) + +table_model = StatsModels.TableRegressionModel(lmm, lmm.mf, lmm.mm) + +emmeans(tm) + +effects(Dict(:period => ["1", "2", "3", "4"]), tm) +``` \ No newline at end of file diff --git a/src/Metida.jl b/src/Metida.jl index 3e5acc6d..567993d2 100644 --- a/src/Metida.jl +++ b/src/Metida.jl @@ -4,19 +4,18 @@ __precompile__() module Metida -using ProgressMeter, LinearAlgebra, ForwardDiff, DiffResults, Random, Optim, LineSearches, MetidaBase#, SparseArrays#, Polyester#, LoopVectorization - -import MetidaBase: StatsBase, StatsModels, CategoricalArrays, Distributions +using ProgressMeter, LinearAlgebra, ForwardDiff, DiffResults, Random, Optim, LineSearches, MetidaBase#, SparseArrays#, Polyester#, LoopVectorization +import StatsBase, StatsModels, Distributions +import MetidaBase: CategoricalArrays, Tables, MetidaModel, AbstractCovarianceStructure, AbstractCovmatMethod, AbstractCovarianceType, AbstractLMMDataBlocks, MetidaTable, metida_table, PrettyTables, indsdict! import MetidaBase.CategoricalArrays: CategoricalArray, AbstractCategoricalVector -import MetidaBase.Distributions: Normal, TDist, FDist, Chisq, MvNormal, FullNormal, ccdf, cdf, quantile - -import MetidaBase: Tables, MetidaModel, AbstractCovarianceStructure, AbstractCovmatMethod, AbstractCovarianceType, AbstractLMMDataBlocks, MetidaTable, metida_table, PrettyTables, indsdict! import MetidaBase.PrettyTables: TextFormat, pretty_table, tf_borderless, ft_printf -import LinearAlgebra:checksquare, BlasFloat + +import Distributions: Normal, TDist, FDist, Chisq, MvNormal, FullNormal, ccdf, cdf, quantile +import LinearAlgebra: checksquare, BlasFloat import StatsModels: @formula, termvars, ModelFrame, FunctionTerm, AbstractTerm, CategoricalTerm, AbstractContrasts, ConstantTerm, InterceptTerm, Term, InteractionTerm, FormulaTerm, ModelMatrix, schema, apply_schema, MatrixTerm, modelcols import StatsBase: fit, fit!, coef, coefnames, confint, nobs, dof_residual, dof, loglikelihood, aic, bic, aicc, isfitted, vcov, mean, var, stderror, modelmatrix, response, responsename, CoefTable, coeftable, crossmodelmatrix -import Base:show, rand, ht_keyindex, getproperty +import Base: show, rand, ht_keyindex, getproperty import Random: default_rng, AbstractRNG, rand! export @formula, @covstr, @lmmformula, diff --git a/src/lmm.jl b/src/lmm.jl index 9d7bfc37..94b16627 100644 --- a/src/lmm.jl +++ b/src/lmm.jl @@ -288,3 +288,9 @@ function Base.getproperty(x::LMM, s::Symbol) end getfield(x, s) end + +#= +function Base.convert(::Type{StatsModels.TableRegressionModel}, lmm::LMM) + StatsModels.TableRegressionModel(lmm, lmm.mf, lmm.mm) +end +=# \ No newline at end of file diff --git a/src/random.jl b/src/random.jl index 3ccaa68a..baf011de 100644 --- a/src/random.jl +++ b/src/random.jl @@ -12,6 +12,11 @@ function rand!(rng::AbstractRNG, v::AbstractVector, lmm::LMM) if !lmm.result.fit error("Model not fitted!") end rand!(rng, v, lmm, lmm.result.theta, lmm.result.beta) end +""" + rand!(v::AbstractVector, lmm::LMM) = rand!(default_rng(), v, lmm, lmm.result.theta, lmm.result.beta) + +Generate random responce vector for fitted 'lmm' model, store results in `v`. +""" rand!(v::AbstractVector, lmm::LMM) = rand!(default_rng(), v, lmm, lmm.result.theta, lmm.result.beta) """ rand(rng::AbstractRNG, lmm::LMM{T}; theta) where T @@ -22,6 +27,11 @@ function rand(rng::AbstractRNG, lmm::LMM{T}, theta::AbstractVector) where T v = Vector{T}(undef, nobs(lmm)) rand!(rng, v, lmm, theta) end +""" + rand!(rng::AbstractRNG, lmm::LMM{T}; theta) where T + +Generate random responce vector 'lmm' model, theta covariance vector, and zero means, store results in `v`. +""" function rand!(rng::AbstractRNG, v::AbstractVector, lmm::LMM{T}, theta::AbstractVector) where T n = length(lmm.covstr.vcovblock) if length(v) != nobs(lmm) error("Wrong v length!") end diff --git a/src/reml.jl b/src/reml.jl index fd833470..fc7e74c3 100644 --- a/src/reml.jl +++ b/src/reml.jl @@ -176,7 +176,7 @@ function reml_sweep_β_nlopt(lmm, data, θ::Vector{T}; maxthreads::Int = 16) whe # θ₂ calculation logdetθ₂ = logdet(Cholesky(ldθ₂, 'U', 0)) # θ₃ calculation - @inbounds @simd for i = 1:n + @inbounds for i = 1:n r = mul!(copy(data.yv[i]), data.xv[i], βtc, -1, 1) vr = LinearAlgebra.LAPACK.potrs!('U', A[i], copy(r)) θ₃ += dot(r, vr) diff --git a/src/varstruct.jl b/src/varstruct.jl index 69382b3b..297c1bb9 100644 --- a/src/varstruct.jl +++ b/src/varstruct.jl @@ -1,6 +1,9 @@ +const CType = Union{FunctionTerm{typeof(+), Vector{Term}}, FunctionTerm{typeof(*), Vector{Term}}, FunctionTerm{typeof(&), Vector{Term}}} + ################################################################################ # @covstr macro ################################################################################ + """ @covstr(ex) @@ -16,8 +19,8 @@ macro covstr(ex) return :(@formula(nothing ~ $ex).rhs) end function modelparse(term::FunctionTerm{typeof(|)}) - eff, subj = term.args_parsed - if !isa(subj, AbstractTerm) throw(FormulaException("Subject term type not <: AbstractTerm. Use `term` or `interaction term` only. Maybe you are using something like this: `@covstr(factor|term1*term2)` or `@covstr(factor|(term1+term2))`. Use only `@covstr(factor|term)` or `@covstr(factor|term1&term2)`.")) end + eff, subj = term.args + if !isa(subj, AbstractTerm) || isa(subj, FunctionTerm{typeof(*), Vector{Term}}) throw(FormulaException("Subject term type not <: AbstractTerm. Use `term` or `interaction term` only. Maybe you are using something like this: `@covstr(factor|term1*term2)` or `@covstr(factor|(term1+term2))`. Use only `@covstr(factor|term)` or `@covstr(factor|term1&term2)`.")) end eff, subj end function modelparse(term) @@ -352,6 +355,7 @@ function fill_coding_dict!(t::T, d::Dict, data) where T <: Term end d end +#= function fill_coding_dict!(t::T, d::Dict, data) where T <: InteractionTerm for i in t.terms if typeof(Tables.getcolumn(data, i.sym)) <: AbstractCategoricalVector || !(typeof(Tables.getcolumn(data, i.sym)) <: AbstractVector{V} where V <: Real) @@ -360,7 +364,8 @@ function fill_coding_dict!(t::T, d::Dict, data) where T <: InteractionTerm end d end -function fill_coding_dict!(t::T, d::Dict, data) where T <: Tuple{Vararg{AbstractTerm}} +=# +function fill_coding_dict_ct!(t, d, data) for i in t if isa(i, Term) if typeof(Tables.getcolumn(data, i.sym)) <: AbstractCategoricalVector || !(typeof(Tables.getcolumn(data, i.sym)) <: AbstractVector{V} where V <: Real) @@ -372,6 +377,40 @@ function fill_coding_dict!(t::T, d::Dict, data) where T <: Tuple{Vararg{Abstract end d end +#= +function fill_coding_dict!(t::T, d::Dict, data) where T <: Tuple{Vararg{AbstractTerm}} + fill_coding_dict_ct!(t, d, data) +end +=# +function fill_coding_dict!(t::T, d::Dict, data) where T <: CType + fill_coding_dict_ct!(t.args, d, data) +end +#= +function fill_coding_dict!(t::T, d::Dict, data) where T <: FunctionTerm{typeof(&), Vector{Term}} + for i in t.args + if isa(i, Term) + if typeof(Tables.getcolumn(data, i.sym)) <: AbstractCategoricalVector || !(typeof(Tables.getcolumn(data, i.sym)) <: AbstractVector{V} where V <: Real) + d[i.sym] = StatsModels.FullDummyCoding() + end + else + fill_coding_dict!(i, d, data) + end + end + d +end +function fill_coding_dict!(t::T, d::Dict, data) where T <: FunctionTerm{typeof(+), Vector{Term}} + for i in t.args + if isa(i, Term) + if typeof(Tables.getcolumn(data, i.sym)) <: AbstractCategoricalVector || !(typeof(Tables.getcolumn(data, i.sym)) <: AbstractVector{V} where V <: Real) + d[i.sym] = StatsModels.FullDummyCoding() + end + else + fill_coding_dict!(i, d, data) + end + end + d +end +=# ################################################################################ # SHOW ################################################################################ diff --git a/test/devtest.jl b/test/devtest.jl index 97b59f1a..d0650e47 100644 --- a/test/devtest.jl +++ b/test/devtest.jl @@ -27,6 +27,14 @@ random = Metida.VarEffect(Metida.@covstr(1 + time|subject&factor), Metida.CSH), ) b11 = @benchmark Metida.fit!($lmm, hes = false; maxthreads = 16) seconds = 15 +#= +lmm = Metida.LMM(@formula(response ~1 + factor2), ftdf3; +repeated = Metida.VarEffect(Metida.@covstr(p|subject), Metida.CSH), +) +@benchmark Metida.fit!($lmm, hes = false; maxthreads = 16) seconds = 15 +=# +#:LN_BOBYQA :LN_NEWUOA +#@benchmark Metida.fit!($lmm, hes = false; maxthreads = 16, solver = :nlopt, optmethod = :LN_NEWUOA) seconds = 15 #@benchmark Metida.fit!($lmm, optmethod = Metida.LBFGS_OM, hes = false; maxthreads = 16) seconds = 15 #@benchmark Metida.fit!($lmm, optmethod = Metida.BFGS_OM, hes = false; maxthreads = 16) seconds = 15 #@benchmark Metida.fit!($lmm, optmethod = Metida.CG_OM, hes = false; maxthreads = 16) seconds = 15 diff --git a/test/test.jl b/test/test.jl index a4024034..6c7f424e 100644 --- a/test/test.jl +++ b/test/test.jl @@ -29,8 +29,6 @@ include("testdata.jl") Metida.fit!(lmm) @test Metida.m2logreml(lmm) ≈ 25.00077786912235 atol=1E-6 - - #Missing lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0m; random = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.DIAG), @@ -39,7 +37,6 @@ include("testdata.jl") @test Metida.m2logreml(lmm) ≈ 16.636012616466203 atol=1E-6 #milmm = Metida.MILMM(lmm, df0m) - #Basic, Subject block lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0; random = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.DIAG), @@ -56,13 +53,11 @@ include("testdata.jl") random = formulation|subject:Metida.DIAG), df0) @test Metida.m2logreml(lmm) ≈ 16.241112644506067 atol=1E-6 - t3table = Metida.typeiii(lmm; ddf = :contain) # NOT VALIDATED t3table = Metida.typeiii(lmm; ddf = :residual) t3table = Metida.typeiii(lmm) ############################################################################ - ############################################################################ # API test ############################################################################ @@ -99,6 +94,7 @@ include("testdata.jl") @test sum(Metida.hessian(lmm)) ≈ 1118.160713481362 atol=1E-2 @test Metida.nblocks(lmm) == 5 @test length(coefnames(lmm)) == 6 + @test Metida.gmatrixipd(lmm) @test Metida.confint(lmm)[end][1] ≈ -0.7630380758015894 atol=1E-4 @test Metida.confint(lmm, 6)[1] ≈ -0.7630380758015894 atol=1E-4 @test Metida.confint(lmm; ddf = :residual)[end][1] ≈ -0.6740837049617738 atol=1E-4 @@ -173,7 +169,6 @@ 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 - end ################################################################################ # df0 @@ -292,7 +287,6 @@ end df0) Metida.fit!(lmm; rholinkf = :psigm) @test Metida.m2logreml(lmm) ≈ 10.065239006121315 atol=1E-6 - end ################################################################################ # ftdf / 1fptime.csv @@ -459,6 +453,7 @@ end random = Metida.VarEffect(Metida.@covstr(1 + r2 * r1|subject), Metida.DIAG; coding=Dict(:r1 => DummyCoding(), :r2 => DummyCoding())) ) Metida.fit!(lmm) + @test Metida.theta(lmm) ≈ [2.796694409004289, 2.900485570555582, 3.354913215348968, 2.0436114769223237, 1.8477830405766895, 2.0436115732330955, 1.0131934233937254] atol=1E-6 # atol=1E-8 ! @test Metida.m2logreml(lmm) ≈ 713.0655862252027 atol=1E-8 end @testset " Model: &, DIAG/SI " begin @@ -466,6 +461,7 @@ end random = Metida.VarEffect(Metida.@covstr(r1&r2|subject), Metida.DIAG), ) Metida.fit!(lmm) + @test Metida.theta(lmm) ≈ [3.0325005960015985, 3.343826588448401, 1.8477830405766895, 1.8477830405766895, 1.8477830405766895, 4.462942536844632, 1.0082345219318216] atol=1E-8 @test Metida.m2logreml(lmm) ≈ 719.9413776641368 atol=1E-8 end @testset " Model: INT, +, TOEPHP(3)/SI " begin @@ -473,6 +469,7 @@ end random = Metida.VarEffect(Metida.@covstr(1 + r1 + r2|subject), Metida.TOEPHP(3); coding = Dict(:r1 => DummyCoding(), :r2 => DummyCoding())), ) Metida.fit!(lmm) + @test Metida.theta(lmm) ≈ [2.843269324925114, 3.3598654954863423, 7.582560427911907e-10, 4.133572859333964, -0.24881591201506625, 0.46067672264107506, 1.0091887333170306] atol=1E-8 @test Metida.m2logreml(lmm) ≈ 705.9946274598822 atol=1E-8 end @testset " Model: TOEP/SI " begin @@ -631,7 +628,6 @@ end ) Metida.fit!(lmm) @test Metida.m2logreml(lmm) ≈ 8.740095378772942 atol=1E-8 - end @testset " Model: Spatial Exponential " begin @@ -662,7 +658,6 @@ end Metida.rand!(v, lmm) Metida.rand!(v, lmm, [4.54797, 2.82342, 1.05771, 0.576979]) Metida.rand!(v, lmm, [4.54797, 2.82342, 1.05771, 0.576979], [44.3, 5.3, 0.5, 0.29]) - end @testset " Show functions " begin @@ -772,8 +767,9 @@ end @test iAs ≈ iAb atol=1E-6 end - @testset " Experimental " begin + + io = IOBuffer(); lmm = Metida.LMM(@formula(r2 ~ f), spatdf; repeated = Metida.VarEffect(Metida.@covstr(x+y|1), Metida.SPEXP), @@ -855,10 +851,12 @@ end @test_throws ErrorException Metida.milmm(lmm; n = 10, verbose = false, rng = StableRNG(1234)) - if !(VERSION < v"1.7") mb = Metida.miboot(mi; n = 10, bootn = 10, double = true, verbose = false, rng = StableRNG(1234)) Base.show(io, mb) end + # Other + @test Metida.varlinkvecapply([0.1, 0.1], [:var, :rho]; varlinkf = :exp, rholinkf = :sigm) ≈ [1.1051709180756477, 0.004999958333749888] atol=1E-6 + end