14
$\begingroup$

Given three variables, y and x, which are positive continuous, and z, which is categorical, I have two candidate models given by:

fit.me <- lmer( y ~ 1 + x + ( 1 + x | factor(z) ) )

and

fit.fe <- lm( y ~ 1 + x )

I hope to compare these models to determine which model is more appropriate. It seems to me that in some sense fit.fe is nested within fit.me. Typically, when this general scenario holds, a chi-squared test can be performed. In R, we can perform this test with the following command,

anova(fit.fe,fit.me)

When both models contain random-effects (generated by lmer from the lme4 package), the anova() command works fine. Owing to boundary parameters, it is normally advisable to test the resulting Chi-Square statistic via simulation, nonetheless, we can still use the statistic in the simulation procedure.

When both models contain only fixed-effects, this approach---and, the associated anova() command---work fine.

However, when one model contains random effects and the reduced model contains only fixed-effects, as in the above scenario, the anova() command doesn't work.

More specifically, I get the following error:

 > anova(fit.fe, fit.me)
 Error: $ operator not defined for this S4 class

Is there anything wrong with using the Chi-Square approach from above (with simulation)? Or is this simply a problem of anova() not knowing how to deal with linear models generated by different functions?

In other words, would it be appropriate to manually generate the Chi-Square statistic derived from the models? If so, what are the appropriate degrees of freedom for comparing these models? By my reckoning:

$$ F = \frac{\left((SSE_{reduced}-SSE_{full})/(p-k)\right)}{\left((SSE_{full})/(n-p-1)\right)} \sim F_{p-k,n-p-1} $$

We are estimating two parameters in the fixed effects model (slope and intercept) and two more parameters (variance parameters for the random slope and random intercept) in the mixed-effects model. Typically, the intercept parameter isn't counted in the degrees of freedom computation, so that implies that $k=1$ and $p=k+2=3$; having said that I'm not sure if the variance parameters for the random-effects parameters should be included in the degrees of freedom computation; the variance estimates for fixed-effect parameters are not considered, but I believe that to be because the parameter estimates for fixed effects are assumed to be unknown constants whilst they are considered to be unknowable random variables for mixed effects. I would appreciate some assistance on this issue.

Finally, does anybody have a more appropriate (R-based) solution to comparing these models?

Improve this question
$\endgroup$
4
  • 5
    $\begingroup$ If you replace lm() with gls() from the nlme package, and lmer() with lme() (again from the nlme package), everything will work fine. But note that you will get a conservative test (too large p-values), as the parameters for the simpler model is on the boundary of the parameter space. And really the choice of whether to include the random effects should be based on theory (e.g., the sampling plan), not on a statistical test. $\endgroup$
    – Karl Ove Hufthammer
    Mar 14, 2015 at 19:58
  • 1
    $\begingroup$ What do you want to do with the models? One model may be better for some purposes and the other model better for other purposes. All models are wrong, so the question is not which model is right, but which is more useful for your particular problem. $\endgroup$
    – Kodiologist
    Mar 14, 2015 at 21:40
  • 1
    $\begingroup$ @Kodiologist Basically, I want to ensure that the parameter estimates for the fixed effects are reliable. The standard errors of which may be unreliable if the observations are assumed to be independent. In addition, it would be nice to make some statement about how the variable the random effect is, but I guess this isn't quite as essential. $\endgroup$
    – user9171
    Mar 14, 2015 at 21:46
  • 2
    $\begingroup$ @user9171 A good way to check stability (reliability) in a model's parameter estimates is to use bootstrapping. Graph bootstrap distributions for each parameter the two models share, with one graph per parameter and model. Tighter distributions imply higher stability. You will likely find that the simpler model yields more stable estimates, because fewer parameters allow for more precise estimation of each parameter. $\endgroup$
    – Kodiologist
    Mar 14, 2015 at 21:54

1 Answer 1

Reset to default
9
$\begingroup$

Technically, you can get it working by just switching the order of the parameters:

> anova(fit.me, fit.fe)

Will work just fine. If you pass an object generated by lmer first, the anova.merMod will be called instead of anova.lm (which does not know how to handle lmer objects). See:

?anova.merMod

Although, choosing a mixed model or a fixed model is a modeling choice which needs to take into account the experimental design, not a model selection problem. See @BenBolker's https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#testing-significance-of-random-effects for more details:

Consider not testing the significance of random effects.

Improve this answer
$\endgroup$
1
  • $\begingroup$ +1. I took the liberty to insert a link to @BenBolker's FAQ that contains some further discussion and references. $\endgroup$
    – amoeba
    Mar 9, 2018 at 11:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.