[R-sig-ME] Point estimate outside of its confidence interval with lmer

[Viechtbauer, Wolfgang (NP)]

After simulating a bunch of datasets, I was able to find one where this issue arises.

##############################

dat <- structure(list(id = c(1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,
4L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 8L, 9L, 9L, 9L, 9L, 10L,
10L, 10L), yij = c(1.3342, -1.3854, -2.3733, -0.2884, -0.567, -0.6045, 0.494,
0.7738, 0.7991, -0.2539, -0.3956, 1.8508, 2.0849, 0.8782, 1.1726, 0.0818,
0.4099, 0.2279, 0.3244, -0.4167, -0.998, -0.9063, 0.9197, 0.806, -0.1283,
4.1764, -0.0946, -0.2065, 0.0359, 1.2685, -0.0226, 0.6875, 0.1928), obs =
c(1L, 2L, 3L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 1L, 1L, 2L, 3L, 4L, 5L, 1L,
2L, 3L, 4L, 5L, 1L, 2L, 1L, 1L, 2L, 3L, 4L, 1L, 2L, 3L )), row.names = c(NA,
-33L), class = "data.frame")

library(lme4)
library(nlme)
library(metafor)

res1 <- lmer(yij ~ 1 + (1 | id), data=dat)
res2 <- lme(yij ~ 1, random = ~ 1 | id, data=dat)
res3 <- rma.mv(yij ~ 1, V=0, random = ~ 1 | id/obs, data=dat)

summary(res1)
summary(res2)
summary(res3)

# so we get the same results from all fits

# profile likelihood CI for sigma based on ML estimation
confint(res1, parm=".sigma")

# the CI for sigma does not include the REML estimate of sigma
sigma(res1)

# get the profile likelihood CI for sigma based on REML estimation
confint(res3, sigma2=2)

# examine the REML profile for sigma^2
ci <- confint(res3, sigma2=2)
profile(res3, sigma2=2, cline=TRUE, steps=100, xlim=c(0.4,2.2))
abline(v=ci$random[1,2:3], lty="dotted")

# note: for this small dataset, the ML fit is really rather different
summary(update(res1, REML=FALSE))

# double-check the profile likelihood CI based on ML estimation
res3 <- update(res3, method="ML")
res3
confint(res3, sigma2=2)

# the lower bound is rather different compared to confint(res1, parm=".sigma")

# examine the ML profile for sigma^2
ci <- confint(res3, sigma2=2)
profile(res3, sigma2=2, cline=TRUE, steps=100, xlim=c(0.4,2.2))
abline(v=ci$random[1,2:3], lty="dotted")

##############################

So, as far as I can tell, there are really two issues here:

1) As noted already, confint.merMod() uses the ML likelihood to profile all parameters including the variance components. This can in principle lead to the CI not including a particular REML estimate of a parameter. However, the switch to ML makes sense in so far as confint.merMod() also profiles fixed effects for which one must use ML estimation. It would be *really* confusing if the CIs for the fixed effects are based on the regular while the CIs for variance components are based on restricted likelihood. Might be nice though to emphasize under help(confint.merMod) that this is happening.

2) However, in the example above, the profiling itself goes wrong, leading to an incorrect lower bound of the profile likelihood CI for sigma. After playing around with some of the settings to profile.merMod(), I got the same CI as metafor with:

confint(res1, parm=".sigma", delta.cutoff=1/64)

So, @Emmanuel, you might want to check if the issue you are encountering is not one of type 1) but rather of type 2).

Best,
Wolfgang

>-----Original Message-----
>From: Emmanuel Curis [mailto:emmanuel.curis using parisdescartes.fr]
>Sent: Friday, 13 May, 2022 17:19
>To: Viechtbauer, Wolfgang (NP)
>Cc: Ben Bolker; r-sig-mixed-models using r-project.org
>Subject: Re: [R-sig-ME] Point estimate outside of its confidence interval with
>lmer
>
>Thanks Ben & Wolfgang for your explainations.
>
>I'll try to prepare an anonymised and selected dataset to check that
>it is the only reason of this discrepancy.
>
>Wolfgang, by bootstrap as a check, I meant that since bootstrapped
>dataset are refited using REML, the sampled sigma are from REML hence
>IC based on them should include the observed one — and this is the
>case, in favor of a REML/ML discrepancy.  But that could also be that
>the khi-square approximation is just false, indeed.  IC are not much
>larger by bootsrap, however, instead shifted when comparing boostrap
>IC for REML and ML fits.
>
>Best regards,
>
>Le Thu, May 12, 2022 at 04:01:04PM +0000, Viechtbauer, Wolfgang (NP) a écrit :
>> Hi Ben,
>>
>> There is nothing unusual about profiling the restricted likelihood and in fact,
>LRTs and profile likelihood CIs for variance components tend to perform slightly
>better. See, for example:
>>
>> Morrell, C. H. (1998). Likelihood ratio testing of variance components in the
>linear mixed-effects model using restricted maximum likelihood. Biometrics,
>54(4), 1560-1568.
>>
>> This is also discussed in various books on mixed-effects models. For example,
>Verbeke & Molenberghs' "Linear Mixed Models for Longitudinal Data" (see section
>6.3.2: "Further, in contrast to the LR test for fixed effects (see Section
>6.2.5), valid LR tests are also obtained under REML estimation. [...]"). So, if
>one can conduct a LRT for variance components under REML estimation, one can also
>profile the restricted likelihood and use this to construct CIs.
>>
>> Of course, the usual caveats apply when the parameter tested falls on the
>boundary of the parameter space.
>>
>> Best,
>> Wolfgang
>>
>> >-----Original Message-----
>> >From: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces using r-project.org] On
>> >Behalf Of Ben Bolker
>> >Sent: Thursday, 12 May, 2022 15:55
>> >To: r-sig-mixed-models using r-project.org
>> >Subject: Re: [R-sig-ME] Point estimate outside of its confidence interval with
>> >lmer
>> >
>> >   For what it's worth, there is a (relatively ancient) Github issue
>> >discussing this:
>> >
>> >https://github.com/lme4/lme4/issues/325
>> >
>> >   Further evidence that we ought to add a warning, or at least a
>> >message, saying that the profiles are based on the ML fit.
>> >
>> >   Can anyone point me to some theory that supports/describes the
>> >process of constructing CIs based on the profile of the REML statistic ... ?
>> >
>> >   cheers
>> >    Ben Bolker
>> >
>> >
>> >On 2022-05-12 7:01 a.m., Viechtbauer, Wolfgang (NP) wrote:
>> >> Dear Emmanuel,
>> >>
>> >> Please see below for my responses.
>> >>
>> >> Best,
>> >> Wolfgang
>> >>
>> >>> -----Original Message-----
>> >>> From: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces using r-project.org]
>On
>> >>> Behalf Of Emmanuel Curis
>> >>> Sent: Thursday, 12 May, 2022 9:30
>> >>> To: r-sig-mixed-models using r-project.org
>> >>> Subject: [R-sig-ME] Point estimate outside of its confidence interval with
>> >lmer
>> >>>
>> >>> Hello,
>> >>>
>> >>> I have encountered an unexpected result for some datasets when using
>confint
>> >>> after fitting a model with lmer: the confidence intervals for the standard
>> >>> deviations in the model did not included the point estimate, given by
>summary
>> >>> for instance.
>> >>>
>> >>> I think the problem is a mix of small sample size and REML vs ML, but I
>would
>> >>> be happy to have confirmation that my interpretation is correct, since I'm
>> >>> not very familiar with profiling and REML vs ML... and wonder if something
>> >>> more problematic is occuring.
>> >>>
>> >>> My interpretation is that lmer fits using REML by default,
>> >>
>> >> Correct.
>> >>
>> >>> hence the variances are estimated "unbiased" (the equivalent of dividing by
>n
>> >>> - k for the residual variance of a linear model, but I'm not sure exactly
>> >>> what is n and k for random effects variance estimation).
>> >>
>> >> Only for some very special cases can we show that the variance components
>are
>> >> estimated unbiasedly. However, yes, REML estimation tends to provide
>estimates
>> >> of variance components that are approximately unbiased in many cases.
>> >>
>> >> Also, for more complex models, changing a ML estimate into a REML one isn't
>> >> simply a multiplication of the ML estimate by some factor that depends on
>> >> things like the sample size (and in more complex models, number of clusters)
>> >> and the number of parameters.
>> >>
>> >>> But when confint is called, using profile, the ML profile is used, hence
>> >>> variance confidence intervals are estimated "biased" (the equivalent
>> >>> of dividing by n).  However, when n is small, dividing by n - k and n
>> >>> may give very different results, hence in some cases the profiled
>> >>> confidence interval for ML estimate does not include the REML
>> >>> estimate, for the standard deviations.  I guess in practice the
>> >>> difference between REML and ML is more complex than just using n or
>> >>> n-k, but would it be idea?
>> >>
>> >> Correct - see above. So I wouldn't just heuristically apply some
>'correction'
>> >> that only comes from simple regression models.
>> >>
>> >> Actually, when constructing a profile likelihood CI for a variance component
>> >> (which is what confint() does by default for lmer model objects), one can
>also
>> >> profile the restricted log likelihood function, which would guarantee that
>the
>> >> estimate falls into the CI. I actually thought that this is what confint()
>> >> does for lmer objects (fitted with REML), but I checked and you are correct
>-
>> >> even when the model was fitted with REML, all profile likelihood CIs
>> >> (including those for variance components) are based on the regular
>likelihood.
>> >>
>> >> Without a reproducible example, I cannot tell you whether this is really the
>> >> reason for the estimate falling outside of the CI, but it certainly could
>be.
>> >>
>> >>> Support for this is that
>> >>> 1) when using bootstrap, there seems to be no such discrepancy
>> >>> 2) when fitting using REML = FALSE, the profiled IC is the same and
>> >>>     does include the (different) point estimate
>> >>>
>> >>> Does it sound correct?
>> >>
>> >> Not sure how 1) helps to diagnose this, but 2) certainly provides support
>for
>> >> the hypothesis that the discrepancy arises from different likelihoods being
>> >> used for estimation and profiling.
>> >>
>> >>> Additionnal questions, if this interpretation is correct:
>> >>> - would it make sense to make confidence intervals based on REML
>> >>>    profiles, and not ML profiles? if so, how?
>> >>
>> >> As noted above, this can be done for variance components, but I don't see a
>> >> way of doing this with confint() for lmer model objects.
>> >>
>> >>> - wouldn't a warning be a good idea when point estimates are outside
>> >>>    CI, with the explaination if it is indeed REML vs ML?
>> >>
>> >> That's up to the developers of lme4 to decide.
>> >>
>> >>> - if this is indeed a "small sample size" problem, I guess in such
>> >>>    cases any asymptotic result is difficult to trust, right?  Does it
>> >>>    mean profiled interval cannot be trusted also, neither nested
>> >>>    models tests, and that only bootstrap may be used?
>> >>
>> >> Profiling the likelihood relies on the asymptotic behavior of the likelihood
>> >> ratio test. So yes, profiling likelihood CIs may not have a nominal coverage
>> >> rate in small samples (for some definition of 'small', which is context
>> >> dependent).
>> >>
>> >>> - in such cases, is there any argument to prefer REML over ML or
>> >>>    vice-versa?
>> >>
>> >> One can make arguments in favor of both ML and REML. REML tends to provide
>> >> approximately unbiased estimates, but with larger mean-squared error, while
>ML
>> >> tends to provide negative biased estimates, but with smaller MSE.
>> >>
>> >>> Thanks in advance for your help,
>> >>>
>> >>> --
>> >>>                                 Emmanuel CURIS
>> >>>                                 emmanuel.curis using parisdescartes.fr
>> >>>
>> >>> Page WWW: http://emmanuel.curis.online.fr/index.html