Fit models with PIRLS and refit without PIRLS yields similar estimates to "slow" model but is faster

[MaximilianNuber]

Dear all,

I am trying to fit many GLMM in parallel, but still needed some speed-up, so I decided to use PIRLS in MixedModels.jl.
However, I discovered that fitting a model with PIRLS (fast = true) and then refitting the model yields parameters which are not quite the same, but much closer to the estimates of the non-PIRLS optimisation. And fitting first with PIRLS, then without, took less time than only fitting without PIRLS.
I wondered, if the estimates from the run with PIRLS remain in the model and are a good starting point for the estimation with NLopt?
Is this comparable to init_from_lmm?
So my question is, if using this would be considered good practice.

Unfortunately, I cannot share my own data, so I worked the following MWE.
The time difference between fitting the full model as compared to fast-fitting, then slow-fitting is more pronounced in my own data.
I would just like to know if I can reasonably take advantage of this?

Thanks for any help.

verbagg = MixedModels.dataset(:verbagg)
const vaform = @formula(r2 ~ 1 + anger + gender + btype + situ + (1|subj) + (1|item));

@time mod_fast = fit(GeneralizedLinearMixedModel, vaform, verbagg, Bernoulli(), fast = true)

Output:

2.348491 seconds (71.50 k allocations: 9.973 MiB)
Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 1)
  r2 ~ 1 + anger + gender + btype + situ + (1 | subj) + (1 | item)
  Distribution: Bernoulli{Float64}
  Link: LogitLink()

   logLik    deviance     AIC       AICc        BIC    
 -4075.7917  8151.5833  8167.5833  8167.6024  8223.0537

Variance components:
        Column   Variance Std.Dev. 
subj (Intercept)  1.794431 1.339564
item (Intercept)  0.246843 0.496833

 Number of obs: 7584; levels of grouping factors: 316, 24

Fixed-effects parameters:
─────────────────────────────────────────────────────
                   Coef.  Std. Error      z  Pr(>|z|)
─────────────────────────────────────────────────────
(Intercept)    0.208273    0.405425    0.51    0.6075
anger          0.0543791   0.0167533   3.25    0.0012
gender: M      0.304089    0.191223    1.59    0.1118
btype: scold  -1.0165      0.257531   -3.95    <1e-04
btype: shout  -2.0218      0.259235   -7.80    <1e-14
situ: self    -1.01344     0.210888   -4.81    <1e-05
─────────────────────────────────────────────────────

Refit with "fast = false":

@btime refit!(mod_fast, fast = false)

Output:

6.589997 seconds (106.80 k allocations: 3.902 MiB)
Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 1)
  r2 ~ 1 + anger + gender + btype + situ + (1 | subj) + (1 | item)
  Distribution: Bernoulli{Float64}
  Link: LogitLink()

   logLik    deviance     AIC       AICc        BIC    
 -4075.6999  8151.3997  8167.3997  8167.4188  8222.8701

Variance components:
        Column   Variance Std.Dev. 
subj (Intercept)  1.794848 1.339719
item (Intercept)  0.245308 0.495286

 Number of obs: 7584; levels of grouping factors: 316, 24

Fixed-effects parameters:
─────────────────────────────────────────────────────
                   Coef.  Std. Error      z  Pr(>|z|)
─────────────────────────────────────────────────────
(Intercept)    0.197208    0.405177    0.49    0.6265
anger          0.0574887   0.0167574   3.43    0.0006
gender: M      0.321003    0.19126     1.68    0.0933
btype: scold  -1.05835     0.256796   -4.12    <1e-04
btype: shout  -2.10498     0.25852    -8.14    <1e-15
situ: self    -1.05525     0.210295   -5.02    <1e-06
─────────────────────────────────────────────────────

Now the full model for comparison:

@time mod_full = fit(GeneralizedLinearMixedModel, vaform, verbagg, Bernoulli(), fast = false)

Output:

10.163416 seconds (221.99 k allocations: 15.436 MiB, 0.94% gc time)
Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 1)
  r2 ~ 1 + anger + gender + btype + situ + (1 | subj) + (1 | item)
  Distribution: Bernoulli{Float64}
  Link: LogitLink()

   logLik    deviance     AIC       AICc        BIC    
 -4075.6999  8151.3997  8167.3997  8167.4187  8222.8701

Variance components:
        Column   Variance Std.Dev. 
subj (Intercept)  1.794811 1.339706
item (Intercept)  0.245319 0.495296

 Number of obs: 7584; levels of grouping factors: 316, 24

Fixed-effects parameters:
─────────────────────────────────────────────────────
                   Coef.  Std. Error      z  Pr(>|z|)
─────────────────────────────────────────────────────
(Intercept)    0.199036    0.405176    0.49    0.6233
anger          0.0574321   0.0167572   3.43    0.0006
gender: M      0.320691    0.191258    1.68    0.0936
btype: scold  -1.05882     0.256801   -4.12    <1e-04
btype: shout  -2.10545     0.258525   -8.14    <1e-15
situ: self    -1.05544     0.210299   -5.02    <1e-06
─────────────────────────────────────────────────────

Here my loaded packages:

using Base.Threads
using GLM
using MixedModels
using StatsBase
using StatsModels
using CategoricalArrays
using DataFrames
using CSV

using Base.Iterators
using SplittablesBase: SplittablesBase, amount, halve
using BenchmarkTools

I am on Julia 1.8.1 on an Ubuntu server.