[R-sig-ME] lmer() producing wildly inaccurate values for some GLM models
Nathaniel Smith
njs at pobox.com
Sat Oct 20 01:54:12 CEST 2007
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
More information about the R-sig-mixed-models
mailing list