From 880b296239a651e0fc875d186c4b06cf22a75b12 Mon Sep 17 00:00:00 2001 From: PharmCat Date: Thu, 26 Jan 2023 15:44:01 +0300 Subject: [PATCH 1/9] sm 0.7 --- Project.toml | 4 ++-- src/reml.jl | 2 +- test/devtest.jl | 8 ++++++++ 3 files changed, 11 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index d1ba0a2d..b4e9efb2 100644 --- a/Project.toml +++ b/Project.toml @@ -21,11 +21,11 @@ StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d" DiffResults = "1" ForwardDiff = "0.10" LineSearches = "7" -MetidaBase = "0.10.1, 0.10.2" +MetidaBase = "0.10.1, 0.10.2, 0.11" Optim = "1" ProgressMeter = "1" StatsBase = "0.29, 0.30, 0.31, 0.32, 0.33" -StatsModels = "0.6" +StatsModels = "0.6, 0.7" julia = "1" [extras] 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/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 From 2f71754caafe0c78296395e5f4e8a0dd51440fcf Mon Sep 17 00:00:00 2001 From: PharmCat Date: Thu, 26 Jan 2023 17:04:17 +0300 Subject: [PATCH 2/9] update to StatsModels 0.7 --- Project.toml | 4 ++-- src/varstruct.jl | 41 ++++++++++++++++++++++++++++++++++++++--- test/test.jl | 15 +++------------ 3 files changed, 43 insertions(+), 17 deletions(-) diff --git a/Project.toml b/Project.toml index b4e9efb2..ea8b8209 100644 --- a/Project.toml +++ b/Project.toml @@ -21,11 +21,11 @@ StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d" DiffResults = "1" ForwardDiff = "0.10" LineSearches = "7" -MetidaBase = "0.10.1, 0.10.2, 0.11" +MetidaBase = "0.11" Optim = "1" ProgressMeter = "1" StatsBase = "0.29, 0.30, 0.31, 0.32, 0.33" -StatsModels = "0.6, 0.7" +StatsModels = "0.7" julia = "1" [extras] diff --git a/src/varstruct.jl b/src/varstruct.jl index 69382b3b..e9f77f73 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) @@ -360,7 +363,7 @@ 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 +375,38 @@ 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/test.jl b/test/test.jl index a4024034..ce5ecaf6 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 ############################################################################ @@ -173,7 +168,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 +286,6 @@ end df0) Metida.fit!(lmm; rholinkf = :psigm) @test Metida.m2logreml(lmm) ≈ 10.065239006121315 atol=1E-6 - end ################################################################################ # ftdf / 1fptime.csv @@ -459,6 +452,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-8 @test Metida.m2logreml(lmm) ≈ 713.0655862252027 atol=1E-8 end @testset " Model: &, DIAG/SI " begin @@ -466,6 +460,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 +468,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 +627,6 @@ end ) Metida.fit!(lmm) @test Metida.m2logreml(lmm) ≈ 8.740095378772942 atol=1E-8 - end @testset " Model: Spatial Exponential " begin @@ -662,7 +657,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,7 +766,6 @@ end @test iAs ≈ iAb atol=1E-6 end - @testset " Experimental " begin io = IOBuffer(); lmm = Metida.LMM(@formula(r2 ~ f), spatdf; @@ -855,10 +848,8 @@ 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 - end From 5988cebcb7cea522381b3941984dd9bd1871ce67 Mon Sep 17 00:00:00 2001 From: PharmCat Date: Wed, 15 Mar 2023 00:56:07 +0300 Subject: [PATCH 3/9] update --- Project.toml | 2 +- docs/src/examples.md | 17 +++++++++++++++++ src/lmm.jl | 6 ++++++ 3 files changed, 24 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index ea8b8209..42200b29 100644 --- a/Project.toml +++ b/Project.toml @@ -3,7 +3,7 @@ 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" diff --git a/docs/src/examples.md b/docs/src/examples.md index 5e9c075c..9c0931f3 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -156,3 +156,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/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 From e0ecc602194dabe2af50bb5bf60c3230a770a1ee Mon Sep 17 00:00:00 2001 From: PharmCat Date: Wed, 15 Mar 2023 01:18:32 +0300 Subject: [PATCH 4/9] update --- Project.toml | 2 ++ src/Metida.jl | 15 +++++++-------- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/Project.toml b/Project.toml index 42200b29..013d4ad7 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ 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,6 +20,7 @@ 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.11" 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, From 0f17718ac3c921f71523505f9720a15731dbae7e Mon Sep 17 00:00:00 2001 From: PharmCat Date: Wed, 15 Mar 2023 01:39:04 +0300 Subject: [PATCH 5/9] update --- test/test.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test.jl b/test/test.jl index ce5ecaf6..c7f6fb4b 100644 --- a/test/test.jl +++ b/test/test.jl @@ -452,7 +452,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-8 + @test Metida.theta(lmm) ≈ [2.796694409004289, 2.900485570555582, 3.354913215348968, 2.0436114769223237, 1.8477830405766895, 2.0436115732330955, 1.0131934233937254] atol=1E-7 # atol=1E-8 ! @test Metida.m2logreml(lmm) ≈ 713.0655862252027 atol=1E-8 end @testset " Model: &, DIAG/SI " begin From 8c5e150fa803c8a87cf4d03fea4fe45663fe2fd3 Mon Sep 17 00:00:00 2001 From: PharmCat Date: Wed, 15 Mar 2023 01:48:59 +0300 Subject: [PATCH 6/9] test update --- test/test.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test.jl b/test/test.jl index c7f6fb4b..ced3de9f 100644 --- a/test/test.jl +++ b/test/test.jl @@ -452,7 +452,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-7 # atol=1E-8 ! + @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 From 4333a3f2fc7632a53190045cbc3b79792621bb44 Mon Sep 17 00:00:00 2001 From: PharmCat Date: Wed, 15 Mar 2023 02:21:24 +0300 Subject: [PATCH 7/9] update --- src/varstruct.jl | 4 ++++ test/test.jl | 7 +++++++ 2 files changed, 11 insertions(+) diff --git a/src/varstruct.jl b/src/varstruct.jl index e9f77f73..297c1bb9 100644 --- a/src/varstruct.jl +++ b/src/varstruct.jl @@ -355,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) @@ -363,6 +364,7 @@ function fill_coding_dict!(t::T, d::Dict, data) where T <: InteractionTerm end d end +=# function fill_coding_dict_ct!(t, d, data) for i in t if isa(i, Term) @@ -375,9 +377,11 @@ function fill_coding_dict_ct!(t, d, data) 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 diff --git a/test/test.jl b/test/test.jl index ced3de9f..6c7f424e 100644 --- a/test/test.jl +++ b/test/test.jl @@ -94,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 @@ -767,6 +768,8 @@ end end @testset " Experimental " begin + + io = IOBuffer(); lmm = Metida.LMM(@formula(r2 ~ f), spatdf; repeated = Metida.VarEffect(Metida.@covstr(x+y|1), Metida.SPEXP), @@ -852,4 +855,8 @@ end 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 From 45f1e6cbe4b9a0cf145b890f38afc44b4bfb5c2e Mon Sep 17 00:00:00 2001 From: PharmCat Date: Sat, 18 Mar 2023 04:11:23 +0300 Subject: [PATCH 8/9] MixedModels --- docs/Project.toml | 4 +--- docs/src/examples.md | 6 ++++++ 2 files changed, 7 insertions(+), 3 deletions(-) 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 9c0931f3..567ec3e1 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -3,6 +3,12 @@ ```@example lmmexample using Metida, CSV, DataFrames, MixedModels, 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 nothing; # hide From 4d5007386083e5f8e4a68ba3b591758ee3905e7a Mon Sep 17 00:00:00 2001 From: PharmCat Date: Sat, 18 Mar 2023 04:22:45 +0300 Subject: [PATCH 9/9] fix --- docs/src/examples.md | 2 +- src/random.jl | 10 ++++++++++ 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/docs/src/examples.md b/docs/src/examples.md index 567ec3e1..21087aff 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -1,7 +1,7 @@ ### Example 1 - Continuous and categorical predictors ```@example lmmexample -using Metida, CSV, DataFrames, MixedModels, CategoricalArrays; +using Metida, CSV, DataFrames, CategoricalArrays; import Pkg Pkg.activate("MixedModels") 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