Random vs. fixed effects

library(nlme)
set.seed(1)
f <- function(n, k) {

n <- 10000
f(n, 4) # 4 levels

######################
f(n, 50) # 50 levels

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