replace fitlog with a column table, always populate it (#850)

* replace fitlog with a column table, always populate it

* NEWS entry

* address new and old fitlog formats in serialization

* docs

* remove fitlog from bootstrap

* no thin

Author: palday

NEWS.md

@@ -3,6 +3,8 @@ MixedModels v5.0.0 Release Notes
 - Optimization is now performed _without constraints_. In a post-fitting step, the Cholesky factor is canonicalized to have non-negative diagonal elements. [#840]
 - The default optimizer has changed to NLopt's implementation of NEWUOA where possible. NLopt's implementation fails on 1-dimensional problems, so in the case of a single, scalar random effect, BOBYQA is used instead. In the future, the default optimizer backend will likely change to PRIMA and NLopt support will be moved to an extension. Blocking this change in backend is an issue with PRIMA.jl when running in VSCode's built-in REPL on Linux. [#840]
 - [BREAKING] Support for constrained optimization has been completely removed, i.e. the field `lowerbd` has been removed from `OptSummary`.
+- [BREAKING] A fitlog is always kept -- the deprecated keyword argument `thin` has been removed as has the `fitlog` keyword argument. [#850]
+- The fitlog is now stored as Tables.jl-compatible column table. [#850]
 - Internal code around the default optimizer has been restructured. In particular, the NLopt backend has been moved to a submodule, which will make it easier to move it to an extension if we promote another backend to the default. [#853]
 - Internal code around optimization in profiling has been restructuring so that fitting done during calls to `profile` respect the `backend` and `optimizer` settings. [#853]
 - The `prfit!` convenience function has been removed. [#853]

docs/src/optimization.md

@@ -190,10 +190,8 @@ DisplayAs.Text(ans) # hide
 ```
 
 More detailed information about the intermediate steps of the nonlinear optimizer can be obtained the `fitlog` field.
-By default, `fitlog` is not populated, but passing the keyword argument `fitlog=true` to `fit!` and `refit!` will result in it being populated with the values obtained at each step of optimization:
 
 ```@example Main
-refit!(fm2; fitlog=true)
 first(fm2.optsum.fitlog, 5)
 DisplayAs.Text(ans) # hide
 ```
@@ -250,7 +248,7 @@ Suppose, for example, that the user wishes to try a [Nelder-Mead](https://en.wik
 ```@example Main
 fm2nm = LinearMixedModel(@formula(reaction ~ 1+days+(1+days|subj)), sleepstudy);
 fm2nm.optsum.optimizer = :LN_NELDERMEAD;
-fit!(fm2nm; fitlog=true)
+fit!(fm2nm)
 fm2nm.optsum
 DisplayAs.Text(ans) # hide
 ```

ext/MixedModelsPRIMAExt.jl

@@ -26,10 +26,10 @@ _optimizer!(::Val{:lincoa}, args...; kwargs...) = PRIMA.lincoa!(args...; kwargs.
 _optimizer!(::Val{:newuoa}, args...; kwargs...) = PRIMA.newuoa!(args...; kwargs...)
 
 function MixedModels.optimize!(m::LinearMixedModel, ::PRIMABackend;
-    progress::Bool=true, fitlog::Bool=false, kwargs...)
+    progress::Bool=true, kwargs...)
     optsum = m.optsum
     prog = ProgressUnknown(; desc="Minimizing", showspeed=true)
-    fitlog && empty!(optsum.fitlog)
+    empty!(optsum.fitlog)
 
     function obj(x)
         val = if x == optsum.initial
@@ -49,15 +49,14 @@ function MixedModels.optimize!(m::LinearMixedModel, ::PRIMABackend;
             end
         end
         progress && ProgressMeter.next!(prog; showvalues=[(:objective, val)])
-        fitlog && push!(optsum.fitlog, (copy(x), val))
+        push!(optsum.fitlog, (; θ=copy(x), objective=val))
         return val
     end
 
     maxfun = optsum.maxfeval > 0 ? optsum.maxfeval : 500 * length(optsum.initial)
     info = _optimizer!(Val(optsum.optimizer), obj, optsum.final;
         maxfun,
         optsum.rhoend, optsum.rhobeg)
-    ProgressMeter.finish!(prog)
     optsum.feval = info.nf
     optsum.fmin = info.fx
     optsum.returnvalue = Symbol(info.status)
@@ -66,12 +65,12 @@ function MixedModels.optimize!(m::LinearMixedModel, ::PRIMABackend;
 end
 
 function MixedModels.optimize!(m::GeneralizedLinearMixedModel, ::PRIMABackend;
-    progress::Bool=true, fitlog::Bool=false,
+    progress::Bool=true,
     fast::Bool=false, verbose::Bool=false, nAGQ=1,
     kwargs...)
     optsum = m.optsum
     prog = ProgressUnknown(; desc="Minimizing", showspeed=true)
-    fitlog && empty!(opstum.fitlog)
+    empty!(optsum.fitlog)
 
     function obj(x)
         val = try
@@ -83,7 +82,7 @@ function MixedModels.optimize!(m::GeneralizedLinearMixedModel, ::PRIMABackend;
             x == optsum.initial && rethrow()
             m.optsum.finitial
         end
-        fitlog && push!(optsum.fitlog, (copy(x), val))
+        push!(optsum.fitlog, (; θ=copy(x), objective=val))
         verbose && println(round(val; digits=5), " ", x)
         progress && ProgressMeter.next!(prog; showvalues=[(:objective, val)])
         return val

src/MixedModelsNLoptExt.jl

@@ -19,11 +19,11 @@ end
 const NLoptBackend = Val{:nlopt}
 
 function MixedModels.optimize!(m::LinearMixedModel, ::NLoptBackend;
-    progress::Bool=true, fitlog::Bool=false,
+    progress::Bool=true,
     kwargs...)
     optsum = m.optsum
     prog = ProgressUnknown(; desc="Minimizing", showspeed=true)
-    fitlog && empty!(optsum.fitlog)
+    empty!(optsum.fitlog)
 
     function obj(x, g)
         isempty(g) || throw(ArgumentError("g should be empty for this objective"))
@@ -44,7 +44,7 @@ function MixedModels.optimize!(m::LinearMixedModel, ::NLoptBackend;
             end
         end
         progress && ProgressMeter.next!(prog; showvalues=[(:objective, val)])
-        fitlog && push!(optsum.fitlog, (copy(x), val))
+        push!(optsum.fitlog, (; θ=copy(x), objective=val))
         return val
     end
 
@@ -59,12 +59,12 @@ function MixedModels.optimize!(m::LinearMixedModel, ::NLoptBackend;
 end
 
 function MixedModels.optimize!(m::GeneralizedLinearMixedModel, ::NLoptBackend;
-    progress::Bool=true, fitlog::Bool=false,
+    progress::Bool=true,
     fast::Bool=false, verbose::Bool=false, nAGQ=1,
     kwargs...)
     optsum = m.optsum
     prog = ProgressUnknown(; desc="Minimizing", showspeed=true)
-    fitlog && empty!(optsum.fitlog)
+    empty!(optsum.fitlog)
 
     function obj(x, g)
         isempty(g) || throw(ArgumentError("g should be empty for this objective"))
@@ -77,7 +77,7 @@ function MixedModels.optimize!(m::GeneralizedLinearMixedModel, ::NLoptBackend;
             x == optsum.initial && rethrow()
             optsum.finitial
         end
-        fitlog && push!(optsum.fitlog, (copy(x), val))
+        push!(optsum.fitlog, (; θ=copy(x), objective=val))
         verbose && println(round(val; digits=5), " ", x)
         progress && ProgressMeter.next!(prog; showvalues=[(:objective, val)])
         return val

src/bootstrap.jl

@@ -254,7 +254,7 @@ function parametricbootstrap(
     )
     samp = replicate(n; progress) do
         simulate!(rng, m; β, σ, θ)
-        refit!(m; progress=false, fitlog=false)
+        refit!(m; progress=false)
         (
             objective=ftype.(m.objective),
             σ=ismissing(m.σ) ? missing : ftype(m.σ),

src/generalizedlinearmixedmodel.jl

@@ -213,7 +213,6 @@ end
 """
     fit!(m::GeneralizedLinearMixedModel; fast=false, nAGQ=1,
                                          verbose=false, progress=true,
-                                         thin::Int=1,
                                          init_from_lmm=Set())
 
 Optimize the objective function for `m`.
@@ -228,9 +227,6 @@ 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.
 
-The `thin` argument is ignored: it had no impact on the final model fit and the logic around
-thinning the `fitlog` was needlessly complicated for a trivial performance gain.
-
 By default, the starting values for model fitting are taken from a (non mixed,
 i.e. marginal ) GLM fit. Experience with larger datasets (many thousands of
 observations and/or hundreds of levels of the grouping variables) has suggested
@@ -256,8 +252,6 @@ function StatsAPI.fit!(
     fast::Bool=false,
     nAGQ::Integer=1,
     progress::Bool=true,
-    thin::Int=typemax(Int),
-    fitlog::Bool=false,
     init_from_lmm=Set(),
     backend::Symbol=m.optsum.backend,
     optimizer::Symbol=m.optsum.optimizer,
@@ -294,7 +288,7 @@ function StatsAPI.fit!(
     optsum.backend = backend
     optsum.optimizer = optimizer
 
-    xmin, fmin = optimize!(m; progress, fitlog, fast, verbose, nAGQ)
+    xmin, fmin = optimize!(m; progress, fast, verbose, nAGQ)
 
     θopt = length(xmin) == length(θ) ? xmin : view(xmin, (length(β) + 1):lastindex(xmin))
     rectify!(m.LMM)                  # flip signs of columns of m.λ elements with negative diagonal els
@@ -312,7 +306,7 @@ function StatsAPI.fit!(
             (fmin + optsum.ftol_zero_abs)
             fmin = zeroobj
             copyto!(xmin, xmin_)
-            fitlog && push!(optsum.fitlog, (copy(xmin), fmin))
+            push!(optsum.fitlog, (; θ=copy(xmin), objective=fmin))
         end
     end
 

src/linearmixedmodel.jl

@@ -443,24 +443,17 @@ end
 
 """
     fit!(m::LinearMixedModel; progress::Bool=true, REML::Bool=m.optsum.REML,
-                              σ::Union{Real, Nothing}=m.optsum.sigma,
-                              thin::Int=typemax(Int),
-                              fitlog::Bool=true)
+                              σ::Union{Real, Nothing}=m.optsum.sigma)
 
 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.
-
-The `thin` argument is ignored: it had no impact on the final model fit and the logic around
-thinning the `fitlog` was needlessly complicated for a trivial performance gain.
 """
 function StatsAPI.fit!(
     m::LinearMixedModel{T};
     progress::Bool=true,
     REML::Bool=m.optsum.REML,
     σ::Union{Real,Nothing}=m.optsum.sigma,
-    thin::Int=typemax(Int),
-    fitlog::Bool=false,
     backend::Symbol=m.optsum.backend,
     optimizer::Symbol=m.optsum.optimizer,
 ) where {T}
@@ -495,7 +488,7 @@ function StatsAPI.fit!(
         optsum.finitial = objective!(m, optsum.initial)
     end
 
-    xmin, fmin = optimize!(m; progress, fitlog)
+    xmin, fmin = optimize!(m; progress)
 
     setθ!(m, xmin)                   # ensure that the parameters saved in m are xmin
     rectify!(m)                      # flip signs of columns of m.λ elements with negative diagonal els
@@ -512,7 +505,7 @@ function StatsAPI.fit!(
         if (zeroobj = objective!(m, xmin_)) ≤ (fmin + optsum.ftol_zero_abs)
             fmin = zeroobj
             copyto!(xmin, xmin_)
-            fitlog && push!(optsum.fitlog, (copy(xmin), fmin))
+            push!(optsum.fitlog, (; θ=copy(xmin), objective=fmin))
         end
     end
 

src/optsummary.jl

@@ -47,7 +47,7 @@ Similarly, the list of applicable optimization parameters can be inspected with
 * `rhoend`: as in PRIMA, not used in NLopt
 
 ## MixedModels-specific fields, unrelated to the optimizer
-* `fitlog`: A vector of tuples of parameter and objectives values from steps in the optimization
+* `fitlog`: A Tables.jl-compatible table with columns `θ` and `objective`. The precise type is an implementation detail
 * `nAGQ`: number of adaptive Gauss-Hermite quadrature points in deviance evaluation for GLMMs
 * `REML`: use the REML criterion for LMM fits
 * `sigma`: a priori value for the residual standard deviation for LMM
@@ -90,7 +90,9 @@ Base.@kwdef mutable struct OptSummary{T<:AbstractFloat}
     rhoend::T = rhobeg / 1_000_000
 
     # not SVector because we would need to parameterize on size (which breaks GLMM)
-    fitlog::Vector{Tuple{Vector{T},T}} = Vector{Tuple{Vector{T},T}}()
+    # we could parameterize more strictly, but this is already a mutable struct
+    # and not isbits, so I don't think we're going to see any real advantage
+    fitlog::Table = Table(; θ=Vector{Vector{T}}(), objective=T[])
     nAGQ::Int = 1
     REML::Bool = false
     sigma::Union{T,Nothing} = nothing
@@ -121,7 +123,7 @@ and `par` is a vector of parameter numbers.
 """
 function Tables.columntable(s::OptSummary; stack::Bool=false)
     fitlog = s.fitlog
-    val = (; iter=axes(fitlog, 1), objective=last.(fitlog), θ=first.(fitlog))
+    val = (; iter=axes(fitlog, 1), objective=fitlog.objective, θ=fitlog.θ)
     stack || return val
     θ1 = first(val.θ)
     k = length(θ1)

src/profile/fixefpr.jl

@@ -54,7 +54,7 @@ function betaprofile!(
     pr::FeProfile{T}, tc::TableColumns{T}, βⱼ::T, j::Integer, obj::T, neg::Bool
 ) where {T}
     prm = pr.m
-    refit!(prm, mul!(copyto!(prm.y, pr.y₀), pr.xⱼ, βⱼ, -1, 1); progress=false, fitlog=false)
+    refit!(prm, mul!(copyto!(prm.y, pr.y₀), pr.xⱼ, βⱼ, -1, 1); progress=false)
     (; positions, v) = tc
     v[1] = (-1)^neg * sqrt(prm.objective - obj)
     getθ!(view(v, positions[:θ]), prm)

src/profile/profile.jl

@@ -33,7 +33,8 @@ Profiling starts at the parameter estimate and continues until reaching a parame
 value of ζ exceeds `threshold`.
 """
 function profile(m::LinearMixedModel; threshold=4)
-    isfitted(m) || refit!(m; progress=false, fitlog=false)
+    isfitted(m) || refit!(m; progress=false)
+    fitlog = copy(m.optsum.fitlog)
     final = copy(m.optsum.final)
     profile = try
         tc = TableColumns(m)
@@ -57,6 +58,7 @@ function profile(m::LinearMixedModel; threshold=4)
         copyto!(m.optsum.final, final)
         m.optsum.fmin = objective(m)
         m.optsum.sigma = nothing
+        m.optsum.fitlog = fitlog
     end
     return profile
 end