reset optsum.final in unfit! (#828)

* reset optsum.final in unfit!

* NEWS, version bump

* tweak test tolerance

* fix somewhat misguided test

* update CI

Author: palday

.github/workflows/minimum.yml

@@ -21,8 +21,7 @@ jobs:
     strategy:
       matrix:
         julia-version: [min]
-        julia-arch: [x64]
-        os: [ubuntu-22.04, macos-13, windows-2019]
+        os: [ubuntu-22.04, macos-14, windows-2019]
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2

NEWS.md

@@ -1,3 +1,7 @@
+MixedModels v4.35.1 Release Notes
+==============================
+- The final parameter vector `optsum.final` is now reset in calls to `unfit!`. This has the secondary effect of correctly starting the fit for `prfit!` at the initial parameter vector instead of the final parameter vector of any previous optimization. [#828]
+
 MixedModels v4.35.0 Release Notes
 ==============================
 - `StatsAPI.cooksdistance(::LinearMixedModel)` is now defined and exported. [#825]
@@ -632,3 +636,4 @@ Package dependencies
 [#815]: https://github.com/JuliaStats/MixedModels.jl/issues/815
 [#823]: https://github.com/JuliaStats/MixedModels.jl/issues/823
 [#825]: https://github.com/JuliaStats/MixedModels.jl/issues/825
+[#828]: https://github.com/JuliaStats/MixedModels.jl/issues/828

Project.toml

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

README.md

@@ -50,13 +50,13 @@ Typical distribution forms are _Bernoulli_ for binary data or _Poisson_ for coun
 
 ## Currently Tested Platforms
 
-|OS      | OS Version    |Arch |Julia           |
-|:------:|:-------------:|:---:|:--------------:|
-|Linux   | Ubuntu 20.04  | x64 |v1.8            |
-|Linux   | Ubuntu 20.04  | x64 |current release |
-|Linux   | Ubuntu 20.04  | x64 |nightly         |
-|macOS   | Monterey 12   | x64 |v1.8            |
-|Windows | Server 2019   | x64 |v1.8            |
+|OS      | OS Version    |Arch    |Julia           |
+|:------:|:-------------:|:------:|:--------------:|
+|Linux   | Ubuntu 22.04  | x64    |v1.10           |
+|Linux   | Ubuntu 22.04  | x64    |current release |
+|Linux   | Ubuntu 22.04  | x64    |nightly         |
+|macOS   | Sonoma 14     | aarm64 |v1.10           |
+|Windows | Server 2019   | x64    |v1.10           |
 
 Note that previous releases still support older Julia versions.
 

src/linearmixedmodel.jl

@@ -1253,8 +1253,13 @@ end
 Mark a model as unfitted.
 """
 function unfit!(model::LinearMixedModel{T}) where {T}
-    model.optsum.feval = -1
-    model.optsum.initial_step = T[]
+    optsum = model.optsum
+    optsum.feval = -1
+    optsum.initial_step = T[]
+    # for variances (bounded at zero), we have ones, while
+    # for everything else (bounded at -Inf), we have zeros
+    map!(T ∘ iszero, optsum.initial, optsum.lowerbd)
+    copyto!(optsum.final, optsum.initial)
     reevaluateAend!(model)
 
     return model

test/pirls.jl

@@ -110,9 +110,9 @@ end
     cbpp = dataset(:cbpp)
     gm2 = fit(MixedModel, first(gfms[:cbpp]), cbpp, Binomial(); wts=float(cbpp.hsz), progress=false, init_from_lmm=[:β, :θ])
     @test weights(gm2) == cbpp.hsz
-    @test deviance(gm2,true) ≈ 100.09585619892968 atol=0.0001
-    @test sum(abs2, gm2.u[1]) ≈ 9.723054788538546 atol=0.0001
-    @test logdet(gm2) ≈ 16.90105378801136 atol=0.001
+    @test deviance(gm2, true) ≈ 100.09585619892968 rtol=0.0001
+    @test sum(abs2, gm2.u[1]) ≈ 9.723054788538546 rtol=0.0001
+    @test logdet(gm2) ≈ 16.90105378801136 rtol=0.0001
     @test isapprox(sum(gm2.resp.devresid), 73.47174762237978, atol=0.001)
     @test isapprox(loglikelihood(gm2), -92.02628186840045, atol=0.001)
     @test !dispersion_parameter(gm2)
@@ -281,7 +281,5 @@ end
     form = @formula(recalled ~ serialpos + zerocorr(serialpos | subject) + (1 | subject & session))
     glmm = @test_logs((:warn, r"Evaluation at default initial parameter vector failed"),
                       GeneralizedLinearMixedModel(form, df, Bernoulli()));
-    glmm.optsum.ftol_rel = 1e-7
-    fit!(glmm; init_from_lmm=[:β, :θ], fast=true, progress=false)
-    @test deviance(glmm) ≈ 996.0402 atol=0.01
+    @test all(==(1e-8), glmm.optsum.initial)
 end

test/prima.jl

@@ -16,14 +16,19 @@ prmodel = LinearMixedModel(formula(model), dataset(:sleepstudy))
 prmodel.optsum.backend = :prima
 
 @testset "$optimizer" for optimizer in (:cobyla, :lincoa)
-    MixedModels.unfit!(prmodel)
+    unfit!(prmodel)
     prmodel.optsum.optimizer = optimizer
     fit!(prmodel; progress=false, fitlog=false)
     @test isapprox(loglikelihood(model), loglikelihood(prmodel))
 end
 
+@testset "refit!" begin
+    refit!(prmodel; progress=false, fitlog=true)
+    @test prmodel.optsum.fitlog[1][1] == [1.0]
+end
+
 @testset "failure" begin
-    MixedModels.unfit!(prmodel)
+    unfit!(prmodel)
     prmodel.optsum.optimizer = :bobyqa
     prmodel.optsum.maxfeval = 5
     @test_logs((:warn, r"PRIMA optimization failure"),