`confint(::MixedModelProfile)` swaps `lower` and `upper` when profiled model is fitted with `REML`

[yjunechoe]

All works fine for models fitted with ML:

m = fit(MixedModel, @formula(reaction ~ 1 + days + (1|subj)), MixedModels.dataset(:sleepstudy));
confint(profile(m))

DictTable with 3 columns and 4 rows:
 par   estimate  lower    upper
 ────┬───────────────────────────
 β1  │ 251.405   231.992  270.818
 β2  │ 10.4673   8.88655  12.048
 σ   │ 30.8954   27.8139  34.5911
 σ1  │ 36.0121   26.0079  52.9358

But with REML, lower and upper values are seemingly swapped for the βs:

refit!(m, REML=true);
confint(profile(m))

DictTable with 3 columns and 4 rows:
 par   estimate  lower    upper
 ────┬───────────────────────────
 β1  │ 251.405   265.901  236.909
 β2  │ 10.4673   11.6974  9.23715
 σ   │ 30.9912   27.8913  34.7109
 σ1  │ 37.1238   26.6233  55.1929

[palday]

Thanks for including an MWE!

I'm looking into what causes this -- both REML and ML models use the same confint method, so this is a little bit surprising.