Merge pull request #33 from PharmCat/dev14

Dev14 (StatsModels 0.7)

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 <[email protected]>"]
-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]

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"

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)
+```

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,

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
+=#

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

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)

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
 ################################################################################

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

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,20 +453,23 @@ 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
     lmm = Metida.LMM(@formula(response ~ 1 + factor), ftdf3;
     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
     lmm = Metida.LMM(@formula(response ~ 1 + factor), ftdf3;
     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