add fitlog

NEWS.md

@@ -1,11 +1,13 @@
 MixedModels v4.1.0 Release Notes
 ========================
 * Add support for specifying a fixed value of `σ`, the residual standard deviation,
-  in `LinearMixedModel`. `fit` takes a keyword-argument `σ`. For models constructed
-  separately, `σ` can be specified by setting `optsum.sigma`. [#551]
-  struct. [#551]
-
-
+  in `LinearMixedModel`. `fit` takes a keyword-argument `σ`. `fit!` does not expose `σ`,
+  but `σ` can be changed after model construction by setting `optsum.sigma`. [#551]
+* Add support for logging the non-linear optimizer's steps via a `thin`
+  keyword-argument for `fit` and `fit!`. The default behavior is 'maximal' thinning,
+  such that only the initial and final values are stored. `OptSummary` has a new field
+  `fitlog` that contains the aforementioned log as a  vector of tuples of parameter and
+  objective values.[#552]
 
 MixedModels v4.0.0 Release Notes
 ========================
@@ -275,3 +277,4 @@ Package dependencies
 [#537]: https://github.com/JuliaStats/MixedModels.jl/issues/537
 [#539]: https://github.com/JuliaStats/MixedModels.jl/issues/539
 [#551]: https://github.com/JuliaStats/MixedModels.jl/issues/551
+[#552]: https://github.com/JuliaStats/MixedModels.jl/issues/552

Project.toml

@@ -1,7 +1,7 @@
 name = "MixedModels"
 uuid = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
 author = ["Phillip Alday <[email protected]>", "Douglas Bates <[email protected]>", "Jose Bayoan Santiago Calderon <[email protected]>"]
-version = "4.0.0"
+version = "4.1.0"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"

docs/src/optimization.md

@@ -189,6 +189,15 @@ fm2.optsum
 DisplayAs.Text(ans) # hide
 ```
 
+More detailed information about the intermediate steps of the nonlinear optimizer can be obtained the `fitlog` field.
+By default, `fitlog` contains entries for only the initial and final steps, but additional information about every nth step can be obtained with the `thin` keyword-argument to `fit`, `fit!` and `refit!`:
+
+```@example Main
+refit!(fm2; thin=1)
+fm2.optsum.fitlog[1:10]
+DisplayAs.Text(ans) # hide
+```
+
 ## A blocked Cholesky factor
 
 A `LinearMixedModel` object contains two blocked matrices; a symmetric matrix `A` (only the lower triangle is stored) and a lower-triangular `L` which is the lower Cholesky factor of the updated and inflated `A`.
@@ -220,18 +229,30 @@ To modify the optimization process the input fields can be changed after constru
 
 Suppose, for example, that the user wishes to try a [Nelder-Mead](https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method) optimization method instead of the default [`BOBYQA`](https://en.wikipedia.org/wiki/BOBYQA) (Bounded Optimization BY Quadratic Approximation) method.
 ```@example Main
-fm2 = LinearMixedModel(@formula(reaction ~ 1+days+(1+days|subj)), sleepstudy);
-fm2.optsum.optimizer = :LN_NELDERMEAD;
-fit!(fm2)
-fm2.optsum
+fm2nm = LinearMixedModel(@formula(reaction ~ 1+days+(1+days|subj)), sleepstudy);
+fm2nm.optsum.optimizer = :LN_NELDERMEAD;
+fit!(fm2nm; thin=1)
+fm2nm.optsum
 DisplayAs.Text(ans) # hide
 ```
 
 The parameter estimates are quite similar to those using `:LN_BOBYQA` but at the expense of 140 functions evaluations for `:LN_NELDERMEAD` versus 57 for `:LN_BOBYQA`.
+When plotting the progress of the individual fits, it becomes obvious that `:LN_BOBYQA` has fully converged by the time `:LN_NELDERMEAD` begins to approach the optimum.
+
+```@example Main
+using Gadfly
+nm = fm2nm.optsum.fitlog
+bob = fm2.optsum.fitlog
+convdf = DataFrame(algorithm=[repeat(["NelderMead"], length(nm));
+                           repeat(["BOBYQA"], length(bob))],
+                   objective=[last.(nm); last.(bob)],
+                   step=[1:length(nm); 1:length(bob)])
+plot(convdf, x=:step, y=:objective, color=:algorithm, Geom.line)
+```
 
 Run time can be constrained with  `maxfeval` and `maxtime`.
 
-See the documentation for the [`NLopt`](https://github.com/JuliaOpt/NLopt.jl) package for details about the various settings. 
+See the documentation for the [`NLopt`](https://github.com/JuliaOpt/NLopt.jl) package for details about the various settings.
 
 ### Convergence to singular covariance matrices
 

src/generalizedlinearmixedmodel.jl

@@ -174,20 +174,22 @@ function fit(
     fast::Bool=false,
     nAGQ::Integer=1,
     progress::Bool=true,
+    thin::Int=typemax(Int),
 )
     return fit(
         GeneralizedLinearMixedModel,
         f,
         columntable(tbl),
         d,
         l;
-        wts=wts,
-        offset=offset,
-        contrasts=contrasts,
-        verbose=verbose,
-        fast=fast,
-        nAGQ=nAGQ,
-        progress=progress,
+        wts,
+        offset,
+        contrasts,
+        verbose,
+        fast,
+        nAGQ,
+        progress,
+        thin,
     )
 end
 
@@ -204,15 +206,17 @@ function fit(
     fast::Bool=false,
     nAGQ::Integer=1,
     progress::Bool=true,
+    thin::Int=typemax(Int),
 )
     return fit!(
         GeneralizedLinearMixedModel(
             f, tbl, d, l; wts=wts, offset=offset, contrasts=contrasts
         );
-        verbose=verbose,
-        fast=fast,
-        nAGQ=nAGQ,
-        progress=progress,
+        verbose,
+        fast,
+        nAGQ,
+        progress,
+        thin,
     )
 end
 
@@ -230,25 +234,29 @@ function fit(
     fast::Bool=false,
     nAGQ::Integer=1,
     progress::Bool=true,
+    thin::Int=typemax(Int),
 )
     return fit(
         GeneralizedLinearMixedModel,
         f,
         tbl,
         d,
         l;
-        wts=wts,
-        contrasts=contrasts,
-        offset=offset,
-        verbose=verbose,
-        fast=fast,
-        nAGQ=nAGQ,
-        progress=progress,
+        wts,
+        contrasts,
+        offset,
+        verbose,
+        fast,
+        nAGQ,
+        progress,
+        thin,
     )
 end
 
 """
-    fit!(m::GeneralizedLinearMixedModel[, verbose=false, fast=false, nAGQ=1, progress=true])
+    fit!(m::GeneralizedLinearMixedModel; fast=false, nAGQ=1,
+                                         verbose=false, progress=true,
+                                         thin::Int=1)
 
 Optimize the objective function for `m`.
 
@@ -261,13 +269,16 @@ during the iterations to minimize the deviance.  There is a delay before this di
 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.
+
+At every `thin`th iteration  is recorded in `fitlog`, optimization progress is saved in `m.optsum.fitlog`.
 """
 function fit!(
     m::GeneralizedLinearMixedModel{T};
     verbose::Bool=false,
     fast::Bool=false,
     nAGQ::Integer=1,
     progress::Bool=true,
+    thin::Int=typemax(Int),
 ) where {T}
     β = m.β
     lm = m.LMM
@@ -284,16 +295,23 @@ function fit!(
     end
     setpar! = fast ? setθ! : setβθ!
     prog = ProgressUnknown("Minimizing"; showspeed=true)
+    # start from zero for the initial call to obj before optimization
+    iter = 0
+    fitlog = optsum.fitlog
     function obj(x, g)
         isempty(g) || throw(ArgumentError("g should be empty for this objective"))
         val = deviance(pirls!(setpar!(m, x), fast, verbose), nAGQ)
+        iszero(rem(iter, thin)) && push!(fitlog, (copy(x), val))
         verbose && println(round(val; digits=5), " ", x)
         progress && ProgressMeter.next!(prog; showvalues=[(:objective, val)])
+        iter += 1
         return val
     end
     opt = Opt(optsum)
     NLopt.min_objective!(opt, obj)
     optsum.finitial = obj(optsum.initial, T[])
+    empty!(fitlog)
+    push!(fitlog, (copy(optsum.initial), optsum.finitial))
     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
@@ -303,10 +321,14 @@ function fit!(
             xmin_[i] = zero(T)
         end
     end
+    loglength = length(fitlog)
     if xmin ≠ xmin_
         if (zeroobj = obj(xmin_, T[])) ≤ (fmin + 1.e-5)
             fmin = zeroobj
             copyto!(xmin, xmin_)
+        elseif length(fitlog) > loglength
+            # remove unused extra log entry
+            pop!(fitlog)
         end
     end
     ## ensure that the parameter values saved in m are xmin

src/linearmixedmodel.jl

@@ -183,9 +183,18 @@ function fit(
     progress::Bool=true,
     REML::Bool=false,
     σ=nothing,
+    thin=typemax(Int),
 )
     return fit(
-        LinearMixedModel, f, Tables.columntable(tbl); wts, contrasts, progress, REML, σ
+        LinearMixedModel,
+        f,
+        Tables.columntable(tbl);
+        wts,
+        contrasts,
+        progress,
+        REML,
+        σ,
+        thin,
     )
 end
 
@@ -198,8 +207,9 @@ function fit(
     progress=true,
     REML=false,
     σ=nothing,
+    thin=typemax(Int),
 )
-    return fit!(LinearMixedModel(f, tbl; contrasts, wts, σ); progress, REML)
+    return fit!(LinearMixedModel(f, tbl; contrasts, wts, σ); progress, REML, thin)
 end
 
 function _offseterr()
@@ -220,11 +230,12 @@ function fit(
     REML::Bool=false,
     offset=[],
     σ=nothing,
+    thin=typemax(Int),
 )
     return if !isempty(offset)
         _offseterr()
     else
-        fit(LinearMixedModel, f, tbl; wts, contrasts, progress, REML, σ)
+        fit(LinearMixedModel, f, tbl; wts, contrasts, progress, REML, σ, thin)
     end
 end
 
@@ -240,13 +251,17 @@ function fit(
     REML::Bool=false,
     offset=[],
     σ=nothing,
-    fast::Bool=false,
-    nAGQ::Integer=1,
+    thin=typemax(Int),
+    fast=nothing,
+    nAGQ=nothing,
 )
     return if !isempty(offset)
         _offseterr()
     else
-        fit(LinearMixedModel, f, tbl; wts, contrasts, progress, REML, σ)
+        if !isnothing(fast) || !isnothing(nAGQ)
+            @warn "fast and nAGQ arguments are ignored when fitting a LinearMixedModel"
+        end
+        fit(LinearMixedModel, f, tbl; wts, contrasts, progress, REML, σ, thin)
     end
 end
 
@@ -398,13 +413,24 @@ function feL(m::LinearMixedModel)
 end
 
 """
-    fit!(m::LinearMixedModel[; progress::Bool=true, REML::Bool=false])
+    fit!(m::LinearMixedModel; progress::Bool=true, REML::Bool=false,
+                              σ::Union{Real, Nothing}=nothing,
+                              thin::Int=typemax(Int))
 
 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.
+
+At every `thin`th iteration  is recorded in `fitlog`, optimization progress is
+saved in `m.optsum.fitlog`.
 """
-function fit!(m::LinearMixedModel{T}; progress::Bool=true, REML::Bool=false) where {T}
+function fit!(
+    m::LinearMixedModel{T};
+    progress::Bool=true,
+    REML::Bool=false,
+    σ::Union{Real,Nothing}=nothing,
+    thin::Int=typemax(Int),
+) 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
@@ -413,14 +439,21 @@ function fit!(m::LinearMixedModel{T}; progress::Bool=true, REML::Bool=false) whe
     opt = Opt(optsum)
     optsum.REML = REML
     prog = ProgressUnknown("Minimizing"; showspeed=true)
+    # start from zero for the initial call to obj before optimization
+    iter = 0
+    fitlog = optsum.fitlog
     function obj(x, g)
         isempty(g) || throw(ArgumentError("g should be empty for this objective"))
         val = objective(updateL!(setθ!(m, x)))
+        iszero(rem(iter, thin)) && push!(fitlog, (copy(x), val))
         progress && ProgressMeter.next!(prog; showvalues=[(:objective, val)])
+        iter += 1
         return val
     end
     NLopt.min_objective!(opt, obj)
     optsum.finitial = obj(optsum.initial, T[])
+    empty!(fitlog)
+    push!(fitlog, (copy(optsum.initial), optsum.finitial))
     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
@@ -431,10 +464,14 @@ function fit!(m::LinearMixedModel{T}; progress::Bool=true, REML::Bool=false) whe
             xmin_[i] = zero(T)
         end
     end
-    if xmin_ ≠ xmin
+    loglength = length(fitlog)
+    if xmin ≠ xmin_
         if (zeroobj = obj(xmin_, T[])) ≤ (fmin + 1.e-5)
             fmin = zeroobj
             copyto!(xmin, xmin_)
+        elseif length(fitlog) > loglength
+            # remove unused extra log entry
+            pop!(fitlog)
         end
     end
     ## ensure that the parameter values saved in m are xmin
@@ -826,7 +863,7 @@ end
 
 Read, check, and restore the `optsum` field from a JSON stream or filename.
 """
-function restoreoptsum!(m::LinearMixedModel, io::IO)
+function restoreoptsum!(m::LinearMixedModel{T}, io::IO) where {T}
     dict = JSON3.read(io)
     ops = m.optsum
     okay =
@@ -852,6 +889,13 @@ function restoreoptsum!(m::LinearMixedModel, io::IO)
     ops.returnvalue = Symbol(dict.returnvalue)
     # provides compatibility with fits saved before the introduction of fixed sigma
     ops.sigma = get(dict, :sigma, nothing)
+    fitlog = get(dict, :fitlog, nothing)
+    ops.fitlog = if isnothing(fitlog)
+        # compat with fits saved before fitlog
+        [(ops.initial, ops.finitial, ops.final, ops.fmin)]
+    else
+        [(convert(Vector{T}, first(entry)), T(last(entry))) for entry in fitlog]
+    end
     return m
 end