[R-sig-ME] lmer Problem

[Ben Bolker]

  I'm taking the liberty of cc'ing this back to r-sig-mixed-models,
**without** the data, because I think the answer will be of general
interest.  You're not doing anything wrong, and I don't think there's
anything pathological about your data -- there's just not quite enough
of it, or it's a little too noisy.  This collapse of the correlation to
+/- 1 is analogous to an estimate of zero variance for a random effect,
just more confusing.  Some of your choices:

  * leave the model as is but treat with caution
  * force slope-intercept correlation to zero by using (1|pid) +
(0+Actual|pid) -- **HOWEVER** as pointed out previously on this list
this is really just a reparameterization of the model, and not
necessarily the best idea
  * drop the slope:group random effect (i.e. just use (1|pid))
  * use blme::blmer to impose a weak prior on the variance-covariance matrix
  * get more data (ha!)


Here's what I did to explore the issue:

####
library(Hmisc)
library(lme4)

load("thorpe_example.RData")

## devel version gives warnings about max|grad| >0.001, but
## this is a false positive: the fit is singular, and we shouldn't
## be checking the gradient under those conditions, or at least
## should be checking it more carefully with respect to the boundary
## (the gradient _parallel_ to the boundary should be zero ...)

summary(dth.fit1 <- lmer(Estimate~Actual+(Actual|pid),data=dth))

pp <- profile(dth.fit1)
library(lattice)
xyplot(pp)

library(ggplot2)
theme_set(theme_bw())
## ggplot complains about 'labelled' objects
dth2 <- transform(subset(dth,select=c(Actual,Estimate,pid)),
                  pid=factor(pid),
                  Actual=as.numeric(Actual),
                  Estimate=as.numeric(Estimate))
ggplot(dth2,aes(Actual,Estimate,colour=pid))+geom_point()+
    geom_smooth(se=FALSE,method="lm")+theme(legend.position="none")

## switching optimizers doesn't help (because we got the right
## answer in the first place)
dth.fit2 <- update(dth.fit1,control=lmerControl(optimizer="Nelder_Mead"))

## centering the predictor doesn't help
dth.fit3 <- update(dth.fit1,
                   data=transform(dth,Actual=Actual-mean(Actual)))


th <- getME(dth.fit1,"theta")
## should th[3] really be exactly zero?
lme4:::checkConv(dth.fit1 at optinfo$derivs,getME(dth.fit1,"theta"),
                 lbound=getME(dth.fit1,"lower"))
## right now check.conv.singular="ignore", but we shouldn't
## be testing gradients and Hessians when we have a singular fit
## (this would be classified as a singular fit since the default
## tolerance for singularity is 1e-4)

## compute and display slice
library(bbmle)
dd <- update(dth.fit1,devFunOnly=TRUE)
ss <- slice2D(getME(dth.fit1,"theta"),dd)
splom(ss)

## plot random effects
dotplot(ranef(dth.fit1,condVar=TRUE),
              scales = list(x = list(relation = 'free')))[["pid"]]

## fit with blme

library(blme)
summary(dth.fit4 <- blmer(Estimate~Actual+(Actual|pid),data=dth))
VarCorr(dth.fit4)


On 14-01-16 10:49 AM, Kevin E. Thorpe wrote:
> Dear Ben.
> 
> My apologies for sending this directly to you off-list, but given that I
> am sending you data I felt it was safest.
> 
> My problem is that in the attached analysis I am getting a correlation
> between the intercept and slope random effects of exactly 1.00.  I
> thought that this might just be rounding in the output, but I think not,
> given the attached caterpillar plot.
> 
> Either I am doing something wrong, or I have yet another pathological
> data set.
> 
> I have included sessionInfo in my output (the Rout file).
> 
> If you could take a look and offer any suggestions, I would be most
> appreciative.
> 
> Thank you,
> 
> Kevin
>