[R-sig-ME] blmer(), minimum amount of prior to get a model to converge

[Ben Bolker]

   I got away from singularity for df>=3 (very large values of df, i.e. 
 >20, are giving crazy answers, not sure why -- please disregard ...).
   You could try values between 2 and 3 if you wanted to be very precise 
(in a little bit of playing I got singularity with df=2.9, not with df=3 
...)

   cheers
     Ben Bolker

## fit initial model
m1 <- lmer(math ~ ses*sector + (ses | sch.id), data = hsb)
## explore singularity: not being reported as such on my machine
## because min theta is slightly greater than 1e-4 threshold ...
## but probably singular, in reality
isSingular(m1)
lwr <- getME(m1, "lower")
theta <- getME(m1, "theta")
theta[lwr==0]
## 0.0004824945
dd <- getME(m1,"devfun")
dd(c(theta[1:2],0)) - dd(theta)
## theta=0 is indeed a better fit ...

## fit baseline wishart
m2 <- blmer(math ~ ses*sector + (ses | sch.id), data = hsb,
             cov.prior = wishart(df=2))

## experiment with a variety of df values
wvec <- c(2:20,(3:10)*10)
res <- matrix(NA,nrow=length(wvec),ncol=3)
for (i in seq_along(wvec)) {
     m <- update(m2,cov.prior=wishart(df=wvec[i]))
     res[i,] <- getME(m,"theta")
}
par(las=1,bty="l")
matplot(wvec,res,type="b",pch=1,log="xy",col=c(1,2,4))


On 10/4/20 3:34 PM, Simon Harmel wrote:
> Dear Vincent,
> 
> Unfortunately, I couldn't improve the singularity problem using the prior
> specification solution you suggested (In fact, I played around with
> `wishart(df = level.dim + 2.5)` more than your suggested range without
> success).
> 
> library(lme4)
> library(blme)
> 
> hsb <- read.csv('
> https://raw.githubusercontent.com/rnorouzian/e/master/hsb.csv')
> <https://raw.githubusercontent.com/rnorouzian/e/master/hsb.csv')m4>
> 
> m1 <- lmer(math ~ ses*sector + (ses | sch.id), data = hsb)
> 
> m2 <- blmer(math ~ ses*sector + (ses | sch.id), data = hsb, cov.prior = ???)
> 
> On Sun, Oct 4, 2020 at 11:10 AM Vincent Dorie <vdorie using gmail.com> wrote:
> 
>> I agree, that that is a meaningful distinction. You can use a prior to
>> nudge the estimate away from the boundary of the space, which
>> addresses singularity. You can also use a prior to add information to
>> the likelihood, which addresses convergence. However, in the latter
>> scenario that information would modify your estimate in a subjective
>> manner, and it would be impossible to say that it was better than
>> simply living with an optimizer warning unless you actually had prior
>> information.
>>
>> On Sat, Oct 3, 2020 at 9:36 PM Ben Bolker <bbolker using gmail.com> wrote:
>>>
>>>     Thanks Vincent.
>>>
>>>     FWIW it would make me really happy if people distinguished clearly
>>> between
>>>
>>> * "singular/nonsingular" - an issue with the 'true' best estimate, i.e.
>>> whether the MLE for the variance-covariance matrix of the REs is
>>> positive definite vs. being only positive *semi*definite
>>>
>>> * "converged/nonconverged" - a question of whether we think the
>>> numerical optimization has worked correctly or not
>>>
>>>     cheers
>>>      Ben Bolker
>>>
>>>
>>> On 10/3/20 9:22 PM, Vincent Dorie wrote:
>>>> There's no single minimum amount, but you can decrease the relative
>>>> impact of the prior by fitting a sequence of models until convergence
>>>> becomes a problem again.
>>>>
>>>> # default
>>>> m2 <- blmer(math ~ ses*sector + (ses | sch.id), data = hsb, cov.prior
>>>> = wishart(df = level.dim + 2.5))
>>>> # point at which blme model is same as lme4
>>>> m3 <- blmer(math ~ ses*sector + (ses | sch.id), data = hsb, cov.prior
>>>> = wishart(df = level.dim + 1))
>>>> # fit models in sequence with df from level.dim + 2.5 to level.dim + 1
>>>>
>>>> Technically, any prior which goes to zero when the determinant of the
>>>> covariance of the random effects go to zero should have the desired
>>>> effect (df > level.dim + 1), but there may be limitations introduced
>>>> by the optimizer.
>>>>
>>>> Vince
>>>>
>>>>
>>>> On Sat, Oct 3, 2020 at 1:17 AM Simon Harmel <sim.harmel using gmail.com>
>> wrote:
>>>>>
>>>>> Hello all,
>>>>>
>>>>> This may be a simple/naive question, but I have a non-converging
>> lmer()
>>>>> model due to singularity.
>>>>>
>>>>> I was wondering what is the minimum prior specification in `blmer()`
>> to get
>>>>> this singular model to converge?
>>>>>
>>>>> library(lme4)
>>>>> library(blme)
>>>>> hsb <- read.csv('
>>>>> https://raw.githubusercontent.com/rnorouzian/e/master/hsb.csv')m4 <-
>> m1 <-
>>>>> lmer(math ~ ses*sector + (ses | sch.id), data = hsb)
>>>>>
>>>>> m2 <- blmer(math ~ ses*sector + (ses | sch.id), data = hsb,
>> cov.prior = ???)
>>>>>
>>>>>           [[alternative HTML version deleted]]
>>>>>
>>>>> _______________________________________________
>>>>> R-sig-mixed-models using r-project.org mailing list
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>>
>>>> _______________________________________________
>>>> R-sig-mixed-models using r-project.org mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>>
>>>
>>> _______________________________________________
>>> R-sig-mixed-models using r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>> _______________________________________________
>> R-sig-mixed-models using r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>