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

[Gabor Grothendieck]

Here it is redone with lme.  lme seems not to exhibit the numerical
problems for 4 levels that we saw with lmer.

> library(nlme)
> 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)
+ }
> n <- 10000
> f(n, 4) # 4 levels
Linear mixed-effects model fit by REML
  Data: NULL
  Log-restricted-likelihood: -14342.06
  Fixed: y ~ x
(Intercept)           x
   1.313495    1.999999

Random effects:
 Formula: ~1 | fac
        (Intercept) Residual
StdDev:    4.421380 1.012295

Number of Observations: 10000
Number of Groups: 4

> ######################
> f(n, 50) # 50 levels
Linear mixed-effects model fit by REML
  Data: NULL
  Log-restricted-likelihood: -14515.87
  Fixed: y ~ x
(Intercept)           x
   1.396288    2.000000

Random effects:
 Formula: ~1 | fac
        (Intercept) Residual
StdDev:    3.322084  1.01249

Number of Observations: 10000
Number of Groups: 50



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
>>
>