By the way, the reason for printing out the number of levels in the factors
as part of the random-effects summary is so that the user will notice when
these are small, which usually means a mistake in the model formula.  Of
course, it requires that the user look at that part of the summary and
anyone who has used SAS is accustomed to ignoring large parts of the output.

Anyway, the issue is now recorded as https://github.com/lme4/lme4/issues/24


On Wed, Dec 12, 2012 at 5:46 PM, Douglas Bates <bates@stat.wisc.edu> wrote:

> Stéphane,
>
> You are correct that there should not be an attempt to fit such a model.
>  In theory the deviance function is degenerate (a singular Hessian matrix
> at every parameter value).  In practice it will not be exactly degenerate
> and the optimizer will attempt to locate a minimum determined by round-off
> error.
>
> Again, it is probably good to throw an error in this case and to warn when
> the number of levels in a grouping factor for random effects term(s) is
> small (say, < 5).  I'll add an issue on the github repository,
> https://github.com/lme4/lme4/issues, regarding this.
>
>
>
> On Wed, Dec 12, 2012 at 3:30 AM, laurent stephane <laurent_step@yahoo.fr>wrote:
>
>> Dear Ben,
>>
>>  In the company where I work we often use mixed models for ANOVA models
>> with random effects, and the parameters of interest are the variances. All
>> my colleagues are SAS users. In such a situation, SAS 9.2 returns a zero
>> variance for the random effect with only one level. My colleagues have
>> often not been paying attention, thus they often reported a zero variance
>> in some statistical reports. Hence the zero variance is not really a good
>> solution.
>>
>>
>> Generally R does not return a negligible variance in such a situation.
>> Here is an example :
>>
>> > # simulate a one-way random effect anova model
>> > set.seed(666)
>> > sigma.b <- 1
>> > sigma.w <- 2
>> > J <- 6
>> > I <- 4
>> > Mu <- 0
>> > y <- c( sapply( rnorm(I, Mu, sigma.b), function(mu) rnorm(J, mu,
>> sigma.w) ) )
>> > dat <- data.frame(y=y, F1=gl(I,J))
>> >
>> > # model 1 :
>> > lmer( y ~ (1|F1), data=dat )
>> Linear mixed model fit by REML
>> Formula: y ~ (1 | F1)
>>    Data: dat
>>    AIC   BIC logLik deviance REMLdev
>>  115.8 119.3 -54.88    110.8   109.8
>> Random effects:
>>  Groups   Name        Variance Std.Dev.
>>  F1       (Intercept) 0.84857  0.92118
>>  Residual             5.53586  2.35284
>> Number of obs: 24, groups: F1, 4
>>
>> Fixed effects:
>>             Estimate Std. Error t value
>> (Intercept)   0.3041     0.6654   0.457
>> >
>> > # add a one-level factor
>> > dat <- cbind(dat, F2="A")
>> >
>> > # model 2 :
>> > # here the F1 variance and the residual variance are identical to model
>> 1
>> > # but the F2 variance has no sense
>> > lmer( y ~ (1|F1) + (1|F2), data=dat )
>> Linear mixed model fit by REML
>> Formula: y ~ (1 | F1) + (1 | F2)
>>    Data: dat
>>    AIC   BIC logLik deviance REMLdev
>>  117.8 122.5 -54.88    111.6   109.8
>> Random effects:
>>  Groups   Name        Variance Std.Dev.
>>  F1       (Intercept) 0.84857  0.92118
>>  F2       (Intercept) 0.61510  0.78428
>>  Residual             5.53586  2.35284
>> Number of obs: 24, groups: F1, 4; F2, 1
>>
>> Fixed effects:
>>             Estimate Std. Error t value
>> (Intercept)   0.3041     1.0285   0.296
>> >
>> > # model 3:
>> > # the F1 variance of model 1 is divided into the F1 variance and the
>> F1:F2 variance
>> > lmer( y ~ (1|F1) + (1|F2) + (1|F1:F2), data=dat )
>> Linear mixed model fit by REML
>> Formula: y ~ (1 | F1) + (1 | F2) + (1 | F1:F2)
>>    Data: dat
>>    AIC   BIC logLik deviance REMLdev
>>  119.8 125.7 -54.88    111.6   109.8
>> Random effects:
>>  Groups   Name        Variance Std.Dev.
>>  F1:F2    (Intercept) 0.42429  0.65137
>>  F1       (Intercept) 0.42429  0.65137
>>  F2       (Intercept) 0.61510  0.78428
>>  Residual             5.53586  2.35284
>> Number of obs: 24, groups: F1:F2, 4; F1, 4; F2, 1
>>
>> Fixed effects:
>>             Estimate Std. Error t value
>> (Intercept)   0.3041     1.0281   0.296
>>
>>
>>
>> All the best,
>> Stéphane
>>
>>
>>
>> ________________________________
>>  De : Ben Bolker <bbolker@gmail.com>
>> À : r-sig-mixed-models@r-project.org
>> Envoyé le : Mercredi 12 décembre 2012 3h58
>> Objet : Re: [R-sig-ME] random effect with only one level
>>
>> laurent stephane <laurent_step@...> writes:
>>
>> >
>> > Hi,
>> >
>> >  The lmer() function seems to work well for a mixed model having a
>> random
>> effect with only one level: it
>> > doesn't crash and there's not even a warning. I'm not a specialist
>> about the
>> numerical mathematics behind
>> > lmer(), but isn't it possible to see somewhere in the output there's a
>> problem
>> such as a convergence
>> > problem for instance ? (in such a situation SAS 9.2 returns a "note"
>> about a
>> problem with some Hessian matrix).
>> >
>> > Stéphane
>> >
>>
>>   I tried out a trivial example.  I agree that it's a bit odd that
>> this doesn't break down anywhere, but the results are mostly sensible --
>> or at least not dangerous.  The random-effects variance reported
>> is meaningless, but the coefficients and residual variances are
>> correct (they match the results of lm()), and the BLUP/conditional
>> mode of the single random-effects level is essentially zero.
>>
>>   Does that agree with your results?
>>
>>   So it might be nice to have a warning, but it doesn't strike
>> me as really dangerous.
>>
>>
>> set.seed(101)
>> d <- data.frame(x=runif(200),f="A")
>> d <- within(d, y <- rnorm(nrow(d),mean=2*x+3))
>>
>> library(lme4)
>> m2 <- lmer(y~x+(1|f),data=d)
>> m3 <- lm(y~x,data=d)
>> all.equal(fixef(m2),coef(m3)) ## TRUE
>> all.equal(attr(VarCorr(m2),"sc"),summary(m3)$sigma)  ## TRUE
>> ranef(m2)  ## 2e-12
>>
>> _______________________________________________
>> R-sig-mixed-models@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>         [[alternative HTML version deleted]]
>>
>>
>> _______________________________________________
>> R-sig-mixed-models@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>>
>

	[[alternative HTML version deleted]]

