[R-sig-ME] Random vs. fixed effects

[Gabor Grothendieck]

Just one other comment so that we are not unfair to lmer.  I was using
the CRAN version of lmer.  Using the development version in package
lme4a there are no numerical problems for the 4 level case. Thus it
appears that lme, lme4a and ADMB all handle this case.

> library(lme4a)
Loading required package: Matrix
Loading required package: lattice

Attaching package: 'Matrix'


        The following object(s) are masked from package:base :

         det

> set.seed(1)
> f <- function(n, k) {
+   set.seed(1)
+   x <- 1:n
+   fac <- gl(k, 1, n)
+   fac.eff <- rnorm(k, 0, 4)[fac]
+   e <- rnorm(n)
+   y <- 1 + 2 * x + fac.eff + e
+   # lme(y ~ x, random = ~ 1 | fac)
+   lmer(y ~ x + (1 | fac))
+ }
> n <- 10000
> f(n, 4)
Linear mixed model fit by REML
Formula: y ~ x + (1 | fac)
 REML
28684

Random effects:
 Groups   Name        Variance Std.Dev.
 fac      (Intercept) 19.5486  4.4214
 Residual              1.0247  1.0123
Number of obs: 10000, groups: fac, 4

Fixed effects:
             Estimate Std. Error t value
(Intercept) 1.313e+00  2.211e+00       1
x           2.000e+00  3.507e-06  570338

Correlation of Fixed Effects:
  (Intr)
x -0.008
> packageDescription("lme4a")$Version
[1] "0.999375-45"


On Fri, Apr 23, 2010 at 12:41 PM, Gabor Grothendieck
<ggrothendieck at gmail.com> wrote:
> Here is a simulation of 10k cases using 4 and 50 level factors for the
> random effect.  With 4 levels there are numerical problems and the
> accuracy of the random effect is terrible.  With 50 levels there are
> no numerical problems and the accuracy is much better.
>
>> library(lme4)
>> set.seed(1)
>> n <- 10000
>> k <- 4
>> f <- function(n, k) {
> + set.seed(1)
> + x <- 1:n
> + fac <- gl(k, 1, n)
> + fac.eff <- rnorm(k, 0, 4)[fac]
> + e <- rnorm(n)
> + y <- 1 + 2 * x + fac.eff + e
> + lmer(y ~ x + (1|fac))
> + }
>
>> # simulation with 4 level random effect
>> f(n, 4)
> Linear mixed model fit by REML
> Formula: y ~ x + (1 | fac)
>   AIC   BIC logLik deviance REMLdev
>  28733 28762 -14363    28702   28725
> Random effects:
>  Groups   Name        Variance Std.Dev.
>  fac      (Intercept) 1.1162   1.0565
>  Residual             1.0298   1.0148
> Number of obs: 10000, groups: fac, 4
>
> Fixed effects:
>             Estimate Std. Error t value
> (Intercept) 1.313e+00  5.286e-01       2
> x           2.000e+00  3.515e-06  568923
>
> Correlation of Fixed Effects:
>  (Intr)
> x -0.033
> Warning message:
> In mer_finalize(ans) : false convergence (8)
>
>> # simulation with 50 level random effect
>> f(n, 50)
> Linear mixed model fit by REML
> Formula: y ~ x + (1 | fac)
>   AIC   BIC logLik deviance REMLdev
>  29040 29069 -14516    29009   29032
> Random effects:
>  Groups   Name        Variance Std.Dev.
>  fac      (Intercept) 11.2016  3.3469
>  Residual              1.0251  1.0125
> Number of obs: 10000, groups: fac, 50
>
> Fixed effects:
>             Estimate Std. Error t value
> (Intercept) 1.396e+00  4.738e-01       3
> x           2.000e+00  3.507e-06  570242
>
> Correlation of Fixed Effects:
>  (Intr)
> x -0.037
>
>
>
>
> On Fri, Apr 23, 2010 at 9:38 AM, Schultz, Mark R. <Mark.Schultz2 at va.gov> wrote:
>> I just read a post by Andrew Dolman suggesting that a factor with only 3
>> levels should be treated as a fixed effect. This seems to be a perennial
>> question with mixed models. I'd really like to hear opinions from
>> several experts as to whether there is a consensus on the topic. It
>> really makes me uncomfortable that such an important modeling decision
>> is made with an "ad hoc" heuristic.
>>
>>
>>
>> Thanks,
>>
>> Mark Schultz, Ph.D.
>>
>> Bedford VA Hospital
>>
>> Bedford, Ma.
>>
>>
>>        [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>