[R-sig-ME] Bug?: var estimates on exp scale for gaussian(link=log)

[Gorjanc Gregor]

Hi!

I would like to fit log-normal fixed model with lmer. Fixed effects estimates
are returned on log scale, while variance components are returned on exp scale.
Why is this the case?

Here is my example. I fit several models to show the behaviour.

Thanks, Gregor

> ## --- GENERATE DATA ---
>
> set.seed(54763278)
> z <- rnorm(n=1000, mean=10, sd=1)
>
> ai <- rep(c(1, 2), each=500)
> a <- c(-2, 2)[ai]
>
> bi <- gl(n=10, k=1, length=1000)
> b <- seq(-5, 5, length.out=10)[bi]
>
> ci <- factor(round(runif(n=1000, min=1, max=100)))
> c <- rnorm(n=100)[ci]
>
> di <- factor(round(runif(n=1000, min=1, max=50)))
> d <- rnorm(n=50)[di]
>
> z <- z + a + b + c + d
> y <- exp(z)
>
> pod <- data.frame(y, ai, bi, ci, di)
>

> ## --- FIT NORMAL MIXED MODEL ---
>
> ## z ~ Normal
> tmp <- lmer(z ~ ai + bi + (1 | ci) + (1 | di), data=pod)
> tmp
...
Random effects:
 Groups   Name        Variance Std.Dev.
 ci       (Intercept) 1.17     1.080
 di       (Intercept) 0.99     0.995
 Residual             1.09     1.042
number of obs: 1000, groups: ci, 100; di, 50

Fixed effects:
            Estimate Std. Error t value
(Intercept)  -1.3399     0.2354    -5.7
ai            4.0293     0.0706    57.1
...
> ## y ~ Log-Normal
> tmp <- lmer(y ~ ai + bi + (1 | ci) + (1 | di), data=pod)
> tmp
...
Random effects:
 Groups   Name        Variance Std.Dev.
 ci       (Intercept) 1.33e+14 11531505
 di       (Intercept) 6.10e+13  7809115
 Residual             7.35e+15 85757196
number of obs: 1000, groups: ci, 100; di, 50

Fixed effects:
             Estimate Std. Error t value
(Intercept) -30997885   12065854   -2.57
ai           20613508    5480798    3.76
...

> ## --- FIT LOG-NORMAL MIXED MODEL ---
>
> ## y ~ Log-Normal
> tmp <- lmer(y ~ ai + bi + (1 | ci) + (1 | di),
+             data=pod, family=gaussian(link="log"))
CHOLMOD warning: w"
·L_ ¿<_ ¿z
Error in devLaplace(PQLpars) :
  Cholmod error `matrix not positive definite' at
file:../Supernodal/t_cholmod_super_numeric.c, line 613

btw, why is this happening, there is no such problem with identity
link?

> ## y ~ Log-Normal
> tmp <- lmer(y ~ ai + bi + (1 | ci),
+             data=pod, family=gaussian(link="log"))
...
Random effects:
 Groups   Name        Variance Std.Dev.
 ci       (Intercept) 4.39e+13  6622360
 Residual             1.64e+14 12824146
number of obs: 1000, groups: ci, 100

Fixed effects:
             Estimate Std. Error  t value
(Intercept) -2.16e+00   2.83e+05 -7.6e-06
ai           4.91e+00   3.99e-01    12.33
...

Here fixed effect estimates are on log scale, while variance component
estimates are huge --> is this on exp scale?