Add a progress meter for LMM fits

.github/workflows/Tier1.yml

@@ -19,7 +19,7 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        julia-version: [1.4, 1.6]
+        julia-version: [1.6]
         julia-arch: [x64]
         os: [ubuntu-18.04, macos-10.15, windows-2019]
     steps:

NEWS.md

@@ -24,6 +24,10 @@ MixedModels v4.0.0 Release Notes
 * `saveoptsum` and `restoreoptsum!` provide for saving and restoring the `optsum`
   field of a `LinearMixedModel` as a JSON file, allowing for recreating a model fit
   that may take a long time for parameter optimization. [#506]
+* Verbose output now uses `ProgressMeter`, which gives useful information about the timing
+  of each iteration and does not flood stdio. The `verbose` argument has been renamed `progress`
+  and the default changed to `true`. [#539]
+* Support for Julia < 1.6 has been dropped. [#539]
 
 Run-time formula syntax
 -----------------------
@@ -257,3 +261,4 @@ Package dependencies
 [#524]: https://github.com/JuliaStats/MixedModels.jl/issues/524
 [#536]: https://github.com/JuliaStats/MixedModels.jl/issues/536
 [#537]: https://github.com/JuliaStats/MixedModels.jl/issues/537
+[#539]: https://github.com/JuliaStats/MixedModels.jl/issues/539

Project.toml

@@ -34,14 +34,14 @@ JSON3 = "1"
 LazyArtifacts = "1"
 NLopt = "0.5, 0.6"
 PooledArrays = "0.5, 1"
-ProgressMeter = "1.5"
+ProgressMeter = "1.7"
 StaticArrays = "0.11, 0.12, 1"
 StatsBase = "0.31, 0.32, 0.33"
 StatsFuns = "0.8, 0.9"
 StatsModels = "0.6.23"
 StructTypes = "1"
 Tables = "1"
-julia = "1.4"
+julia = "1.6"
 
 [extras]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"

README.md

@@ -49,11 +49,11 @@ Typical distribution forms are _Bernoulli_ for binary data or _Poisson_ for coun
 
 |OS|OS Version|Arch|Julia|Tier|
 |:-:|:-:|:-:|:-:|:-:|
-|Linux|Ubuntu 18.04|x64|v1.4, 1.5, 1.6|1|
-|macOS|Catalina 10.15|x64|v1.4, 1.5, 1.6|1|
-|Windows|Server 2019|x64|v1.4, 1.5, 1.6|1|
+|Linux|Ubuntu 18.04|x64|v1.6|1|
+|macOS|Catalina 10.15|x64v1.6|1|
+|Windows|Server 2019|x64|v1.6|1|
 
-Upon release of the next Julia LTS, Tier 1 will become Tier 2 and Julia LTS will become Tier 1.
+Note that previous releases still support older Julia versions.
 
 ## Version 3.0.0
 

src/generalizedlinearmixedmodel.jl

@@ -171,6 +171,7 @@ fit(
     verbose::Bool = false,
     fast::Bool = false,
     nAGQ::Integer = 1,
+    progress::Bool = true,
 ) = fit(
     GeneralizedLinearMixedModel,
     f,
@@ -183,6 +184,7 @@ fit(
     verbose = verbose,
     fast = fast,
     nAGQ = nAGQ,
+    progress = progress,
 )
 
 fit(
@@ -197,6 +199,7 @@ fit(
     verbose::Bool = false,
     fast::Bool = false,
     nAGQ::Integer = 1,
+    progress::Bool = true,
 ) = fit!(
     GeneralizedLinearMixedModel(
         f,
@@ -210,6 +213,7 @@ fit(
     verbose = verbose,
     fast = fast,
     nAGQ = nAGQ,
+    progress = progress,
 )
 
 
@@ -226,6 +230,7 @@ fit(
     REML::Bool = false,
     fast::Bool = false,
     nAGQ::Integer = 1,
+    progress::Bool = true,
 ) = fit(
     GeneralizedLinearMixedModel,
     f,
@@ -238,22 +243,30 @@ fit(
     verbose = verbose,
     fast = fast,
     nAGQ = nAGQ,
+    progress = progress,
 )
 
 """
-    fit!(m::GeneralizedLinearMixedModel[, verbose = false, fast = false, nAGQ = 1])
+    fit!(m::GeneralizedLinearMixedModel[, verbose=false, fast=false, nAGQ=1, progress=true])
 
 Optimize the objective function for `m`.
 
 When `fast` is `true` a potentially much faster but slightly less accurate algorithm, in
 which `pirls!` optimizes both the random effects and the fixed-effects parameters,
 is used.
+
+If `progress` is `true`, the default, a `ProgressMeter.ProgressUnknown` counter is displayed.
+during the iterations to minimize the deviance.  There is a delay before this display is initialized
+and it may not be shown at all for models that are optimized quickly.
+
+If `verbose` is `true`, then both the intermediate results of both the nonlinear optimization and PIRLS are also displayed on standard output.
 """
 function fit!(
     m::GeneralizedLinearMixedModel{T};
     verbose::Bool = false,
     fast::Bool = false,
     nAGQ::Integer = 1,
+    progress::Bool = true,
 ) where {T}
     β = m.β
     lm = m.LMM
@@ -269,16 +282,19 @@ function fit!(
         optsum.final = copy(optsum.initial)
     end
     setpar! = fast ? setθ! : setβθ!
+    prog = ProgressUnknown("Minimizing"; showspeed=true)
     function obj(x, g)
         isempty(g) || throw(ArgumentError("g should be empty for this objective"))
         val = deviance(pirls!(setpar!(m, x), fast, verbose), nAGQ)
         verbose && println(round(val, digits = 5), " ", x)
+        progress && ProgressMeter.next!(prog; showvalues = [(:objective, val),])
         val
     end
     opt = Opt(optsum)
     NLopt.min_objective!(opt, obj)
     optsum.finitial = obj(optsum.initial, T[])
     fmin, xmin, ret = NLopt.optimize(opt, copyto!(optsum.final, optsum.initial))
+    ProgressMeter.finish!(prog)
     ## check if very small parameter values bounded below by zero can be set to zero
     xmin_ = copy(xmin)
     for i in eachindex(xmin_)
@@ -561,33 +577,30 @@ LinearAlgebra.rank(m::GeneralizedLinearMixedModel) = m.LMM.feterm.rank
 
 """
     refit!(m::GeneralizedLinearMixedModel[, y::Vector];
-          fast::Bool = (length(m.θ) == length(m.optsum.final)),
-          nAGQ::Integer = m.optsum.nAGQ))
+           fast::Bool = (length(m.θ) == length(m.optsum.final)),
+           nAGQ::Integer = m.optsum.nAGQ,
+           kwargs...)
 
 Refit the model `m` after installing response `y`.
 
 If `y` is omitted the current response vector is used.
 
 If not specified, the `fast` and `nAGQ` options from the previous fit are used.
-
+`kwargs` are the same as [`fit!`](@ref)
 """
-function refit!(m::GeneralizedLinearMixedModel{T};
+function refit!(m::GeneralizedLinearMixedModel;
                 fast::Bool = (length(m.θ) == length(m.optsum.final)),
-                nAGQ::Integer = m.optsum.nAGQ)  where T
-
-    fit!(unfit!(m); fast=fast, nAGQ=nAGQ)
+                nAGQ::Integer = m.optsum.nAGQ, kwargs...)
+    fit!(unfit!(m); fast=fast, nAGQ=nAGQ, kwargs...)
 end
 
-function refit!(m::GeneralizedLinearMixedModel{T}, y;
-                fast::Bool = (length(m.θ) == length(m.optsum.final)),
-                nAGQ::Integer = m.optsum.nAGQ) where T
+function refit!(m::GeneralizedLinearMixedModel, y; kwargs...)
     m_resp_y = m.resp.y
     length(y) == size(m_resp_y, 1) || throw(DimensionMismatch(""))
     copyto!(m_resp_y, y)
-    refit!(m)
+    refit!(m; kwargs...)
 end
 
-
 """
     setβθ!(m::GeneralizedLinearMixedModel, v)
 

src/linearmixedmodel.jl

@@ -166,15 +166,15 @@ fit(
     tbl;
     wts = [],
     contrasts = Dict{Symbol,Any}(),
-    verbose::Bool = false,
+    progress::Bool = true,
     REML::Bool = false,
 ) = fit(
     LinearMixedModel,
     f,
     Tables.columntable(tbl),
     wts = wts,
     contrasts = contrasts,
-    verbose = verbose,
+    progress = progress,
     REML = REML,
 )
 
@@ -184,11 +184,11 @@ fit(
     tbl::Tables.ColumnTable;
     wts = wts,
     contrasts = contrasts,
-    verbose = verbose,
+    progress = progress,
     REML = REML,
 ) = fit!(
     LinearMixedModel(f, tbl, contrasts = contrasts, wts = wts),
-    verbose = verbose,
+    progress = progress,
     REML = REML,
 )
 
@@ -201,7 +201,7 @@ fit(
     tbl;
     wts = [],
     contrasts = Dict{Symbol,Any}(),
-    verbose::Bool = false,
+    progress::Bool = true,
     REML::Bool = false,
     offset = [],
 ) = !isempty(offset) ? _offseterr() : fit(
@@ -210,7 +210,7 @@ fit(
     tbl,
     wts = wts,
     contrasts = contrasts,
-    verbose = verbose,
+    progress = progress,
     REML = REML,
 )
 
@@ -222,7 +222,7 @@ fit(
     l::IdentityLink;
     wts = [],
     contrasts = Dict{Symbol,Any}(),
-    verbose::Bool = false,
+    progress::Bool = true,
     REML::Bool = false,
     offset = [],
     fast::Bool = false,
@@ -233,7 +233,7 @@ fit(
     tbl,
     wts = wts,
     contrasts = contrasts,
-    verbose = verbose,
+    progress = progress,
     REML = REML,
 )
 
@@ -384,28 +384,31 @@ function feL(m::LinearMixedModel)
 end
 
 """
-    fit!(m::LinearMixedModel[; verbose::Bool=false, REML::Bool=false])
+    fit!(m::LinearMixedModel[; progress::Bool=true, REML::Bool=false])
 
-Optimize the objective of a `LinearMixedModel`.  When `verbose` is `true` the values of the
-objective and the parameters are printed on stdout at each function evaluation.
+Optimize the objective of a `LinearMixedModel`.  When `progress` is `true` a
+`ProgressMeter.ProgressUnknown` display is shown during the optimization of the
+objective, if the optimization takes more than one second or so.
 """
-function fit!(m::LinearMixedModel{T}; verbose::Bool = false, REML::Bool = false) where {T}
+function fit!(m::LinearMixedModel{T}; progress::Bool=true, REML::Bool=false) where {T}
     optsum = m.optsum
     # this doesn't matter for LMM, but it does for GLMM, so let's be consistent
     if optsum.feval > 0
         throw(ArgumentError("This model has already been fitted. Use refit!() instead."))
     end
     opt = Opt(optsum)
     optsum.REML = REML
+    prog = ProgressUnknown("Minimizing"; showspeed=true)
     function obj(x, g)
         isempty(g) || throw(ArgumentError("g should be empty for this objective"))
         val = objective(updateL!(setθ!(m, x)))
-        verbose && println(round(val, digits = 5), " ", x)
+        progress && ProgressMeter.next!(prog; showvalues = [(:objective, val),])
         val
     end
     NLopt.min_objective!(opt, obj)
     optsum.finitial = obj(optsum.initial, T[])
     fmin, xmin, ret = NLopt.optimize!(opt, copyto!(optsum.final, optsum.initial))
+    ProgressMeter.finish!(prog)
     ## check if small non-negative parameter values can be set to zero
     xmin_ = copy(xmin)
     lb = optsum.lowerbd
@@ -772,21 +775,22 @@ function reevaluateAend!(m::LinearMixedModel)
 end
 
 """
-    refit!(m::LinearMixedModel[, y::Vector]; REML=m.optsum.REML)
+    refit!(m::LinearMixedModel[, y::Vector]; REML=m.optsum.REML, kwargs...)
 
 Refit the model `m` after installing response `y`.
 
 If `y` is omitted the current response vector is used.
+`kwargs` are the same as [`fit!`](@ref).
 """
-function refit!(m::LinearMixedModel; REML=m.optsum.REML)
-    fit!(unfit!(m); REML=REML)
+function refit!(m::LinearMixedModel; REML=m.optsum.REML, kwargs...)
+    fit!(unfit!(m); REML=REML, kwargs...)
 end
 
-function refit!(m::LinearMixedModel, y; REML=m.optsum.REML)
+function refit!(m::LinearMixedModel, y; kwargs...)
     resp = m.y
     length(y) == length(resp) || throw(DimensionMismatch(""))
     copyto!(resp, y)
-    refit!(m; REML=REML)
+    refit!(m; kwargs...)
 end
 
 """

test/FactorReTerm.jl

@@ -117,8 +117,8 @@ end
         # equality of RE terms not defined so check that they generate same modelcols
         @test modelcols(ff1.rhs[end], slp) == modelcols(ff2.rhs[end], slp)
 
-        m1 = fit(LMM, f1, slp)
-        m2 = fit(LMM, f2, slp)
+        m1 = fit(LMM, f1, slp; progress=false)
+        m2 = fit(LMM, f2, slp; progress=false)
         @test all(m1.λ .== m2.λ)
 
         @test StatsModels.terms(f2.rhs[end]) == [one, d, s]
@@ -140,8 +140,8 @@ end
         # test that zerocorr actually worked
         @test mc1.inds == mc2.inds == [1, 4]
 
-        m1 = fit(LMM, f1, slp)
-        m2 = fit(LMM, f2, slp)
+        m1 = fit(LMM, f1, slp; progress=false)
+        m2 = fit(LMM, f2, slp; progress=false)
         @test all(m1.λ .== m2.λ)
 
         @test StatsModels.terms(f2.rhs[end]) == [one, d, s]
@@ -263,10 +263,10 @@ end
 
         # fitted model to test amalgamate and fnames, and equivalence with other formulations
         psts = dataset("pastes")
-        m = fit(MixedModel, @formula(strength ~ 1 + (1|batch/cask)), psts)
-        m2 = fit(MixedModel, @formula(strength ~ 1 + (1|batch) + (1|batch&cask)), psts)
-        mr = fit(MixedModel, term(:strength) ~ term(1) + (term(1)|term(:batch)/term(:cask)), psts)
-        m2r = fit(MixedModel, term(:strength) ~ term(1) + (term(1)|term(:batch)) + (term(1)|term(:batch)&term(:cask)), psts)
+        m = fit(MixedModel, @formula(strength ~ 1 + (1|batch/cask)), psts; progress=false)
+        m2 = fit(MixedModel, @formula(strength ~ 1 + (1|batch) + (1|batch&cask)), psts; progress=false)
+        mr = fit(MixedModel, term(:strength) ~ term(1) + (term(1)|term(:batch)/term(:cask)), psts; progress=false)
+        m2r = fit(MixedModel, term(:strength) ~ term(1) + (term(1)|term(:batch)) + (term(1)|term(:batch)&term(:cask)), psts; progress=false)
 
         @test fnames(m) == fnames(m2) == fnames(mr) == fnames(m2r) == (Symbol("batch & cask"), :batch)
         @test m.λ == m2.λ == mr.λ == m2r.λ

test/bootstrap.jl

@@ -39,14 +39,14 @@ include("modelcache.jl")
         center(v::AbstractVector) = v .- (sum(v) / length(v))
         grouseticks = DataFrame(dataset(:grouseticks))
         grouseticks.ch = center(grouseticks.height)
-        gm4 = fit(MixedModel, only(gfms[:grouseticks]), grouseticks, Poisson(), fast=true)  # fails in pirls! with fast=false
+        gm4 = fit(MixedModel, only(gfms[:grouseticks]), grouseticks, Poisson(), fast=true, progress=false)  # fails in pirls! with fast=false
         gm4sim = refit!(simulate!(StableRNG(42), deepcopy(gm4)))
         @test isapprox(gm4.β, gm4sim.β; atol=norm(stderror(gm4)))
     end
 
     @testset "Binomial" begin
         cbpp = dataset(:cbpp)
-        gm2 = fit(MixedModel, first(gfms[:cbpp]), cbpp, Binomial(), wts=float(cbpp.hsz))
+        gm2 = fit(MixedModel, first(gfms[:cbpp]), cbpp, Binomial(), wts=float(cbpp.hsz), progress=false)
         gm2sim = refit!(simulate!(StableRNG(42), deepcopy(gm2)), fast=true)
         @test isapprox(gm2.β, gm2sim.β; atol=norm(stderror(gm2)))
     end
@@ -102,7 +102,7 @@ end
         contra = dataset(:contra)
         # need a model with fast=false to test that we only
         # copy the optimizer constraints for θ and not β
-        gm0 = fit(MixedModel, only(gfms[:contra]), contra, Bernoulli(), fast=false)
+        gm0 = fit(MixedModel, only(gfms[:contra]), contra, Bernoulli(), fast=false, progress=false)
         bs = parametricbootstrap(StableRNG(42), 100, gm0)
         # make sure we're not copying
         @test length(bs.lowerbd) == length(gm0.θ)