Mismatching results for singular fit with different R/lme4 versions

Question

I am trying to match the estimate of random effects from R version 3.5.3 (lme4 1.1-18-1) to R version 4.1.1 (lme4 1.1-27.1). However, there is a small difference of random effects between these two versions when there is singular fit. I'm fine with singularity warnings, but it is puzzling that different versions of R/lme4 produce slightly different results.

The following scripts are from R version 3.5.3 (lme4 1.1-18-1) and R version 4.1.1 (lme4 1.1-27.1) with the dataset Arabidopsis from lme4.

> sessionInfo()
R version 3.5.3 (2019-03-11)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252 
[2] LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
 [1] minqa_1.2.4     MASS_7.3-51.1   compiler_3.5.3  Matrix_1.2-15  
 [5] tools_3.5.3     Rcpp_1.0.1      splines_3.5.3   nlme_3.1-137   
 [9] grid_3.5.3      nloptr_1.2.1    lme4_1.1-18-1   lattice_0.20-38
> library(lme4)
Loading required package: Matrix
> options(digits = 15)
> ##########
> #Example1#
> ##########
> fit1 <- lmer(total.fruits~(1|reg)+(1|reg:popu),data=Arabidopsis,control=lmerControl(optimizer="bobyqa"))
> VarCorr(fit1)
 Groups   Name        Std.Dev.       
 reg:popu (Intercept)  7.744768797534
 reg      (Intercept) 10.629179104291
 Residual             39.028818969641
> ##########
> #Example2#
> ##########
> fit2 <- lmer(total.fruits~(1|reg)+(1|reg:popu)+(1|reg:popu:amd)+(1|reg:popu:amd:status),data=Arabidopsis,control=lmerControl(optimizer="bobyqa"))
> fit2@theta
[1] 0.150979711638631 0.000000000000000 0.189968995915902
[4] 0.260818869156072
> VarCorr(fit2)
 Groups              Name        Std.Dev.       
 reg:popu:amd:status (Intercept)  5.841181759473
 reg:popu:amd        (Intercept)  0.000000000000
 reg:popu            (Intercept)  7.349619506926
 reg                 (Intercept) 10.090696322743
 Residual                        38.688521100461
> ##########
> #Example3#
> ##########
> devfun353 <- lmer(total.fruits~(1|reg)+(1|reg:popu)+(1|reg:popu:amd)+(1|reg:popu:amd:status),data=Arabidopsis,control=lmerControl(optimizer="bobyqa"),devFunOnly = T)
> save.image('myEnvironment353.Rdata')

> sessionInfo()
R version 4.1.1 (2021-08-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
 [1] minqa_1.2.4        MASS_7.3-54        compiler_4.1.1     minque_2.0.0       Matrix_1.3-4      
 [6] tools_4.1.1        Rcpp_1.0.7         tinytex_0.34       splines_4.1.1      nlme_3.1-152      
[11] grid_4.1.1         xfun_0.27          nloptr_1.2.2.2     boot_1.3-28        lme4_1.1-27.1     
[16] ADDutil_2.2.1.9005 lattice_0.20-44   
> library(lme4)
Loading required package: Matrix
Warning message:
package ‘lme4’ was built under R version 4.1.2 
> options(digits = 15)
> ##########
> #Example1#
> ##########
> fit1 <- lmer(total.fruits~(1|reg)+(1|reg:popu),data=Arabidopsis,control=lmerControl(optimizer="bobyqa"))
> VarCorr(fit1)
 Groups   Name        Std.Dev.       
 reg:popu (Intercept)  7.744768797534
 reg      (Intercept) 10.629179104291
 Residual             39.028818969641
> ##########
> #Example2#
> ##########
> fit2 <- lmer(total.fruits~(1|reg)+(1|reg:popu)+(1|reg:popu:amd)+(1|reg:popu:amd:status),data=Arabidopsis,control=lmerControl(optimizer="bobyqa"))
boundary (singular) fit: see ?isSingular
> fit2@theta
[1] 0.150979743348540 0.000000000000000 0.189969036985684 0.260818797487214
> VarCorr(fit2)
 Groups              Name        Std.Dev.       
 reg:popu:amd:status (Intercept)  5.841182965248
 reg:popu:amd        (Intercept)  0.000000000000
 reg:popu            (Intercept)  7.349621069388
 reg                 (Intercept) 10.090693513643
 Residual                        38.688520961140
> ##########
> #Example3#
> ##########
> devfun411 <- lmer(total.fruits~(1|reg)+(1|reg:popu)+(1|reg:popu:amd)+(1|reg:popu:amd:status),data=Arabidopsis,control=lmerControl(optimizer="bobyqa"),devFunOnly = T)
> load('myEnvironment353.Rdata')
> devfun353 <- lme4:::mkdevfun(environment(devfun353))
> minqa::bobyqa(c(1,1,1,1),devfun353,0,control = list(iprint=2))
npt = 6 , n =  4 
rhobeg =  0.2 , rhoend =  2e-07 
start par. =  1 1 1 1 fn =  6443.44054431489 
rho:    0.020 eval:  11 fn:      6393.61 par: 0.00000 0.621363 0.744867 0.823498 
rho:   0.0020 eval:  38 fn:      6361.97 par:0.156855  0.00000 0.190090 0.234676 
rho:  0.00020 eval:  49 fn:      6361.94 par:0.150719  0.00000 0.190593 0.249106 
rho:  2.0e-05 eval:  67 fn:      6361.94 par:0.150988  0.00000 0.189943 0.260821 
rho:  2.0e-06 eval:  74 fn:      6361.94 par:0.150980  0.00000 0.189965 0.260811 
rho:  2.0e-07 eval:  82 fn:      6361.94 par:0.150980  0.00000 0.189969 0.260819 
At return
eval:  90 fn:      6361.9381 par: 0.150980  0.00000 0.189969 0.260819
parameter estimates: 0.150979722854965, 0, 0.189968942342717, 0.260818725554898 
objective: 6361.93810274656 
number of function evaluations: 90 
> minqa::bobyqa(c(1,1,1,1),devfun411,0,control = list(iprint=2))
npt = 6 , n =  4 
rhobeg =  0.2 , rhoend =  2e-07 
start par. =  1 1 1 1 fn =  6443.44054431489 
rho:    0.020 eval:  11 fn:      6393.61 par: 0.00000 0.621363 0.744867 0.823498 
rho:   0.0020 eval:  38 fn:      6361.97 par:0.156855  0.00000 0.190090 0.234676 
rho:  0.00020 eval:  49 fn:      6361.94 par:0.150719  0.00000 0.190593 0.249106 
rho:  2.0e-05 eval:  67 fn:      6361.94 par:0.150988  0.00000 0.189943 0.260821 
rho:  2.0e-06 eval:  74 fn:      6361.94 par:0.150980  0.00000 0.189965 0.260811 
rho:  2.0e-07 eval:  82 fn:      6361.94 par:0.150980  0.00000 0.189969 0.260819 
At return
eval:  90 fn:      6361.9381 par: 0.150980  0.00000 0.189969 0.260819
parameter estimates: 0.150979722854965, 0, 0.189968942342717, 0.260818725554898 
objective: 6361.93810274656 
number of function evaluations: 90 

When the model is simpler, there is no singularity warning and the results match. (See example 1 in both scripts) When model is relatively complex, there is singularity warning and the results are slightly off (See example 2 in both scripts). The difference is <1e-5 in this case but I have observed <1e-4 before. Can anyone shed some lights on why the results are slightly different? and is it even possible to match the results to at least 1e-8?

Not sure if this is useful but I also extract devfun from 3.5.3 and run it in 4.1.1. The results match. (see example 3) In addition, when I read iteration history from BOBYQA, the $\theta$ of the term that leads to singularity warning oscillates between 0 and small numbers (around 1e-7 to 1e-9).

This post discusses similar topics. It also shows the singularity warning leads to slightly different estimate. There is no obvious change in LME4 NEWS that cause the difference. This FAQ and ?isSingular give great explanation on singularity warning but does not address the issue of mismatching directly.

TL;DR: Sometimes when there is singularity warning (I am ok with), the random effects are slightly different under different R/lme4 versions. Why is this happening and how to address it?