[R-sig-ME] lmer() producing wildly inaccurate values for some GLM models

[Nathaniel Smith]

A few more notes on this, in response to a followup off-list:

On Fri, Oct 19, 2007 at 01:44:15AM -0700, Nathaniel Smith wrote:
> for (s in 1:subjects) {
>   subj <- c(subj, rep(s, samples))
>   x.subj <- rgamma(samples, 2, scale=2)
>   x <- c(x, x.subj)
>   i.subj <- rnorm(1, i, sqrt(i.var))
>   b.subj <- rnorm(1, b, sqrt(b.var))
>   true.scale <- c(true.scale, (i.subj + b*x.subj)/true.shape)

The 'b' in the last line line above should, of course, instead be
b.subj.  (This is from the second test case.)

With that fix, the output changes a bit, but not much:

---------------------------->8----------------------------
              (Intercept)        x        logLik
gaussian.lm      192.1596 2.474115 -59410.907309
gaussian.lmer    192.1226 2.483921 -59226.456758
gamma.glm        192.3719 2.420455 -58651.253027
gamma.lmer       192.4162 2.398446  -1039.161105
ig.glm           192.4670 2.395605 -58895.259401
ig.lmer          192.4670 2.395606     -6.385455
> VarCorr(gaussian.lmer)$subj
2 x 2 Matrix of class "dpoMatrix"
            (Intercept)         x
(Intercept)   252.30861 12.853766
x              12.85377  0.698928
> VarCorr(gamma.lmer)$subj
2 x 2 Matrix of class "dpoMatrix"
            (Intercept)         x
(Intercept)   39.109410 2.1247739
x              2.124774 0.1177906
> VarCorr(ig.lmer)$subj
2 x 2 Matrix of class "dpoMatrix"
              (Intercept)             x
(Intercept)  1.529936e-06 -6.192435e-20
x           -6.192435e-20  5.148065e-13
---------------------------->8----------------------------

The logLik is still low, the estimated effects are still off, etc.  

It was also pointed out that of course I'm fitting a model that allows
the intercept and slope to covary, which is not true in my original
data.  One would hope that lmer would discover this and set the
effect correlations to 0, but in any case, if I fit the more
restricted model, then I get a better gaussian fit but
still-anomalous results from the others:

gaussian.restricted <- lmer(y ~ x + (1 | subj) + (0 + x | subj))
gamma.restricted <- lmer(y ~ x + (1 | subj) + (0 + x | subj),
                         family=Gamma(link="identity"))
ig.restricted <- lmer(y ~ x + (1 | subj) + (0 + x | subj),
                      family=inverse.gaussian(link="identity"))

Gaussian model:
Random effects:
 Groups   Name        Variance  Std.Dev.
 subj     (Intercept)  299.7609 17.3136 
 subj     x              1.0649  1.0320 
 Residual             8136.4169 90.2021 
number of obs: 10000, groups: subj, 10; subj, 10

Gamma model:
Random effects:
 Groups   Name        Variance   Std.Dev.  
 subj     (Intercept) 5.8623e+01 7.6565e+00
 subj     x           9.9887e-11 9.9943e-06
 Residual             1.9977e-01 4.4696e-01
number of obs: 10000, groups: subj, 10; subj, 10

Inverse gaussian model:
Random effects:
 Groups   Name        Variance   Std.Dev.  
 subj     (Intercept) 1.5280e-06 1.2361e-03
 subj     x           5.1481e-13 7.1750e-07
 Residual             1.0296e-03 3.2088e-02
number of obs: 10000, groups: subj, 10; subj, 10

-- Nathaniel

-- 
Eternity is very long, especially towards the end.
  -- Woody Allen