[R-sig-ME] Random vs. fixed effects
Gabor Grothendieck
ggrothendieck at gmail.com
Sat Apr 24 19:34:25 CEST 2010
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
>>
>
More information about the R-sig-mixed-models
mailing list